\section{Introduction}

Models of mass transfer of substances (solutes) through semi-perme\-able membranes appear in various contexts, such as biomedical and chemical engineering applications~\cite{Friedman}. Examples include the modelling of electrokinetic flows, solute dynamics across arterial walls, and cellular signal transduction (see, e.g., \cite{Brera2010,Zunino,Cangiani_Natalini_2010} and the references therein).

This work is concerned with the development and analysis of numerical methods for a class of continuum models for mass transfer based on initial/boundary value multi-compart\-ment partial differential equation (PDE) problems, closed by nonlinear interface conditions. The interface conditions considered are the Kedem-Katchalsky (KK) equations, which represent an established model for the mass transfer mechanisms~\cite{KK58,Katchalsky}. More specifically, we consider a generic model problem consisting of a system of semilinear advection-diffusion-reaction parabolic PDE problems in multi-compartment configurations, coupled with nonlinear interface KK-type conditions. The focus is to address some challenges in the numerical solution of these models, such as the treatment of nonlinearities due to both the interface modelling and the nonlinear reactions, the discontinuity of the state variables across the interface, as well as the development of stable numerical methods in the advection-dominated regime.

Numerical methods for mass transfer problems based on conforming finite elements have been developed for the solution of solute dynamics across arterial walls; see~\cite{Zunino,Quarte02,Quarte05} and the references therein for more details. Some existence results for the purely diffusing interface problem without forcing, coupled with KK-type interface conditions, along with some numerical experiments are given in~\cite{Calabro06}. Further, numerical approaches to the treatment of interface conditions for PDE problems, resulting to globally continuous solutions can be found, e.g., in~\cite{babuska,AkrivisDougalis91,MR1622502,MR2013126,Melenk_interface10}.

Discontinuous Galerkin (dG) methods (see, e.g., \cite{nodal_dg_book,riviere_book,dipietro_ern_book} and the references therein)
and mortar methods (see, e.g., \cite{mortar_survey} for a survey) have also been proposed for the treatment of coupled systems via interface conditions in various contexts~\cite{Stenberg98,girault_riviere,electroporation,ern_moz,MR2208566,kan_riv,MR1975269,Hansbo05}. Also, the advantages of dG methods for interfacing different numerical methods (numerical interfaces) have been identified~\cite{per_scho,MR1956028}, as well as their use on transmission-type/high-contrast problems, yielding continuous solutions across the transmission interface, has been investigated~\cite{MR2002258,MR2257119,MR2491426,MR2837483}.

Here, we consider a dG method of interior penalty type for the solution of the semilinear parabolic advection-diffusion-reaction PDE system coupled with nonlinear interface conditions of KK-type across the subdomains. The use of dG is motivated partly by the observation that the interface conditions, yielding \emph{discontinuous} solutions across the interface, can be imposed by modifying the interior penalty dG numerical fluxes. Another important factor for employing a dG method is the desired stability property of the numerical method in the advection-dominated regime.  A priori bounds for the  proposed spatially discrete dG method in both the $L^\infty(L^2)$- and $L^2(H^1)$-type norms are presented for a range of reaction fields, under the simplifying assumption that the finite element mesh is aligned with the subdomain interfaces.

A priori error bounds for interior penalty dG methods for parabolic problems have been considered in various settings (see, e.g.,~\cite{riviere_book} for an exposition and the more recent~\cite{Dolej08}).  DG methods for semilinear parabolic spatially self-adjoint problems with locally Lipschitz continuous nonlinearity have been analysed in~\cite{Lasis07}.
In the present analysis, advection terms are included and systems of equations are considered. In the presence of advection, the analysis of the symmetric interior penalty dG method in~\cite{Lasis07} would require the assumption of quasi-uniformity of the mesh. To avoid this assumption, a different continuation argument is employed in the derivation of the a priori bounds presented here, at the expense of a stricter growth condition on the nonlinearity of the forcing term. This continuation argument is inspired by the derivation of a posteriori bounds for semilinear parabolic phase-field models~\cite{nochetto,bartels}.  The nonlinear interface terms are tackled using a non-standard elliptic projection which is inspired by a classical construction of Douglas and Dupont~\cite{douglas_dupont} for the treatment of nonlinear boundary conditions.

The paper is organised as follows. In Section~\ref{prelim}, the PDE model along with a short derivation of the nonlinear interface conditions is presented. Section~\ref{dg_heat} is devoted to the description of the dG method proposed for the advection-diffusion part of the spatial operator incorporating the nonlinear interface conditions. Section~\ref{err_ell_problem} contains error estimates for the new elliptic projection, which is, in turn, utilised in the subsequent a priori error analysis presented in Section~\ref{err_par_prob}. Section~\ref{sec:numex} contains numerical experiments highlighting the stability and the optimal rate of convergence of the proposed method in practice. Finally, some conclusions are offered in Section~\ref{conclusions}.

\section{Interface modelling and governing PDEs}\label{prelim}

We shall consider systems of parabolic semilinear PDEs describing the flux of solutes around and through a semi-permeable membrane. The membrane is modelled as an internal boundary equipped with nonlinear interface conditions which are described in the following section.

\subsection{Interface modelling}\label{sec:tras}

We outline the Kedem-Katchalsky (KK) equations modelling solutes flow across semi-permeable membranes. The KK equations have been introduced in~\cite{KK58}; we refer to~\cite{Teorell53} for earlier works and to~\cite{Friedman} for a thorough exposition.

It is assumed that the membrane separates two compartments $\Omega^1$ and $\Omega^2$ filled with a free fluid which is called the solvent and that the membrane characteristics are uniform in space and time. The KK equations specify the dependence of the solutes and solvent fluxes across the membrane in terms of two driving forces, namely the hydrostatic and osmotic pressure jumps. In the case of a single solute, the solvent and solute fluxes from $\Omega^1$ to $\Omega^2$ normal to the membrane walls are given by, respectively,
\begin{eqnarray}\label{eq:KK1}
    \vspace{2 mm}
    J_v&=& L_P (\delta {\rm p}-\sigma_d \, \delta\pi),\\
    \label{eq:KK2}
    J_s&=&\omega \, \delta\pi+J_v(1-\sigma_f)\bar{u},
\end{eqnarray}
in terms of the hydrostatic and osmotic pressure jumps $\delta {\rm p}={\rm p}^1-{\rm p}^2$ and $\delta \pi=RT\delta u$, with $R$ being the ideal gas constant, $T$ denoting temperature, and $\delta u=u^1-u^2$ the solute concentration jump across the membrane.  Here, $\bar{u}$ represents the average concentration of the solute across the membrane.  The above constitutive laws are characterised by the phenomenological coefficients of filtration $L_P$, reflection $\sigma_d$ and $\sigma_f$, and permeation $\omega$. These coefficients may depend on the concentration while they are assumed to be constant with respect to both the position along the membrane and time.

Equation~\eqref{eq:KK1}, which is known as Starling's law of filtration, shows that the solvent flow is affected by the osmotic flow of the solute. This observation introduces a nonlinearity in the transport term of the solute flux in equation~\eqref{eq:KK2}. Indeed, substituting $J_v$ into~\eqref{eq:KK2} we get the final model for the solute flux:
\[
J_s=\omega \, R \, T \, \delta u + L_P ( \delta {\rm p} - \sigma_d R T \delta u)(1-\sigma_f) \bar{u}
=p(u^1,u^2) \, \delta u- r    \bar{u}(\mbf{b}\cdot \mbf{n})|_{\Omega^2}.
\]
In this last expression, we have collected the diffusive part of the flux by introducing the nonlinear {\em permeability} function $p(u^1,u^2)$, and written the advective transport in terms of the friction coefficient $r\in [0,1]$ and the normal component of the transport field yielded by the hydrostatic pressure, denoted by $\mbf{b}\cdot \mbf{n}$.

In presence of $n$ solutes, it is usually assumed that there is no direct coupling between the solutes' fluxes~\cite{Friedman}. Under this simplifying assumption the KK equations describing the flux of the solvent and $n$ solutes read
\begin{eqnarray}\label{eq:KK1_many}
    J_v&=&L_P \big(\delta {\rm p}-\sum_{j=1}^n\sigma_{d,j} \, \delta\pi_j\big),\\
    \label{eq:KK2_many}
    J_{s,i} &=&\omega \, \delta\pi_i+J_v(1-\sigma_f)\bar{u}_i,\qquad i=1,\dots ,n,
\end{eqnarray}
with $J_{s,i}$, $\delta \pi_i$, and $\bar{u}_i$ representing, respectively, the flux, osmotic pressure, and average concentration of the $i$-th solute.
Proceeding as in the case of one solute we get the following expression for the solute fluxes:
\begin{eqnarray}
\nonumber
J_{s,i}&=&
p_i(u_i^1,u_i^2)
 \delta u_i-\sum_{j=1}^n \tilde{p}_{i,j}(\mbf{u}^1,\mbf{u}^2) \, \delta u_j
-r_i \bar{u}_i (\mbf{b}\cdot \mbf{n})|_{\Omega^2}
\\
\label{eq:KK2_final}
&=&
\mbf{p}_i(\mbf{u}^1,\mbf{u}^2)\cdot \delta\mbf{u}
-r_i \bar{u}_i (\mbf{b}\cdot \mbf{n})|_{\Omega^2},
\qquad i=1,\dots,n
\end{eqnarray}
with $\delta\mbf{u}=\mbf{u}^1-\mbf{u}^2$, $\mbf{u}^j=(u_1^j,\dots,u_n^j)^{{\mathsf T}}$, $j=1,2$, where $u_i^j$ denotes the concentration of the $i$-th solute in the $j$-th compartment,  $p_i$ the corresponding permeability, and $\tilde{p}_{i,j}$ the cross-coefficients expressing the crowding effect inside the membrane. For the second equality, we have collected in the vector {\em permeability} function $\mbf{p}_i(\mbf{u}^1,\mbf{u}^2)$ all diffusion terms.

It remains to fix a model for the average concentrations inside the membrane.  In the case of relatively thick membranes~\cite{Zunino,Quarte05}, it is appropriate to consider a mechanical approach describing solutes flow within the membrane through an advection-diffusion model. This approach yields a weighted arithmetic average of the concentration at the membrane's faces. Thus, denoting by $u^1_i$ and $u^2_i$, $i=1,\dots, n$, such concentrations, we get
\begin{equation} \label{eq:tiltedaverage}
\bar{u}_i=\delta_w u_i:=v^1_iu^1_i+v^2_iu^2_i,
\end{equation}
for some given weights $v_i=(v^1_i,v^2_i)$ with $v^1_i+v^2_i=1$. These can be expressed in terms of the internal P\`eclet number of the membrane advection-diffusion model and are so that the upwind value dominates~\cite{Zunino}, {\em cf.} the conditions given
below in~\eqref{eq:weights_conditions}.

In what follows, the fluxes given by~\eqref{eq:KK2_final} together with the model \eqref{eq:tiltedaverage} for the average concentration inside the membrane are used to close the PDE problem with appropriate interface conditions.

\subsection{Notation}\label{not}

We denote by $L^p(\omega)$, $1\le p\le +\infty$, the standard Lebesgue spaces, $\omega\subset\mathbb{R}^d$, $d \in \{ 2,3\}$, with corresponding norms $\|\cdot\|_{p,\omega}$; if $\omega=\Omega$, we shall write instead $\|\cdot\|_{p}$. Also the norm of $L^2(\omega)$ will be denoted by $\ltwo{\cdot}{\omega}$ and if $\omega=\Omega$ by $\ltwo{\cdot}{}$ for brevity; by $\langle \cdot,\cdot\rangle$ we write the standard $L^2$-inner product on $\Omega$; when the arguments are vectors of $L^2$-functions, the $L^2$-inner product is modified in the standard fashion.  We denote by $H^s(\omega)$ the standard Hilbertian Sobolev space of index $s\in\mathbb{R}$ of real-valued functions defined on $\omega\subset\mathbb{R}^d$; in particular $H^1_0(\omega)$ signifies the space of functions in $H^1(\omega)$ whose traces onto the boundary $\partial\omega$ vanish.  For $1\le p\le +\infty$, we denote the standard Bochner spaces $L^p(0,T;X)$, with $X$ being a Banach space with norm $\|\cdot\|_X$. Finally, we denote by $C(0,T;X)$ the space of continuous functions $v: [0,T]\to X$ with norm $\|v\|_{C(0,T;X)}:=\max_{0\le t\le T}\|v(t)\|_{X}<+\infty$.

Let $\Omega$ be a bounded open domain with Lipschitz boundary in $\mathbb{R}^d$, and let $\partial\Omega$ be the boundary of $\Omega$. The domain $\Omega$ is subdivided into two subdomains $\Omega^1$ and $\Omega^2$, such that $\Omega=\Omega^1\cup\Omega^2\cup\Gamma_{\inter}$, where $\Gamma_{\inter}:=(\partial {\Omega}^1\cap\partial{\Omega}^2)\backslash \partial\Omega$, see Figure~\ref{fig:two_compartments}. For $i=1,2$, we assume that $\Omega^i$ has Lipschitz boundary and that $\partial \Omega^i\cap\partial\Omega$ has positive $(d-1)$-dimensional (Hausdorff) measure. We define $\mathbb{H}^s:=[H^s(\Omega^1\cup\Omega^2)]^{n}$, $s \in \mathbb{R}$.
\begin{figure}
\centering
\setlength{\unitlength}{2.5cm}
\begin{picture}(0,1.2)
\put(-1,0){\line(0,1){1}}
\put(-1,0){\line(1,0){2}}
\put(-1,1){\line(1,0){2}}
\put(1,0){\line(0,1){1}}
\put(0,0){\line(0,1){1}}
\put(-.5,.5){\makebox(0,0){$\Omega^1$}}
\put(.5,.5){\makebox(0,0){$\Omega^2$}}
\put(.1,.5){\makebox(0,0){{{$\Gamma_{\inter}$}}}}
\end{picture}
\caption{The domain of solution $\Omega$ is subdivided into two subdomains $\Omega^1$, $\Omega^2$. The interface boundary is defined as
$\Gamma_{\inter}:=(\partial {\Omega}^1\cap\partial{\Omega}^2)\backslash \partial\Omega$.}
\label{fig:two_compartments}
\end{figure}

We shall employ the following notational convention: vectors are indicated with lower case bold symbols, $n\times n$ diagonal matrices with upper case (non-bold) symbols, and $n\times d$ tensors with upper case bold symbols.

The gradient $\nabla {\mbf{v}}$ of a vector function ${\mbf{v}}:\Omega^1\cup\Omega^2\to\mathbb{R}^n$ in $\mathbb{H}^1$ is a mapping $\Omega^1\cup\Omega^2\to\mathbb{R}^{n\times d}$ gained from componentwise application of the gradient operation: $\nabla {\mbf{v}}:=(\nabla v_1,\dots,\nabla v_n)^{{\mathsf T}}$. Similarly the divergence $\nabla \cdot\mbf{Q}$ of the tensor-valued function $\mbf{Q}:\Omega^1\cup\Omega^2\to\mathbb{R}^{n\times d}$ is $\nabla \cdot\mbf{Q}:=(\nabla \cdot Q_1,\dots,\nabla \cdot Q_n)^{{\mathsf T}}$ where the $Q_i$ are rows of $\mbf Q$.

\subsection{Model Problem}\label{mod}

For a time interval $[0,T]$, $T>0$, and for
\[
{\mbf{u}}\in L^2(0,T;\mathbb{H}^1), \qquad {\mbf{u}}_t\in L^2(0,T;\mathbb{H}^{-1}),
\]
with ${\mbf{u}}:=(u_1,\dots,u_n)^{{\mathsf T}}$, we consider the system of semilinear parabolic equations
\begin{equation}\label{pde}
 {\mbf{u}}_t - \nabla\cdot (A\nabla {\mbf{u}} - U\mbf{B})+\mbf{F}({\mbf{u}})=\mbf{0}\quad\text{in } (0,T]\times(\Omega^1\cup\Omega^2),
\end{equation}
with $U$ denoting the diagonal matrix $U=\diag(u_1,\dots,u_n)$.  Here, $\mbf{B}$ is an $n\times d$ tensor field with rows $B_i\in C^1(0,T;W^{1,\infty}({\Omega}\backslash\Gamma_{\inter})^d\cap W^{\infty}({\rm div},\Omega))$, $i=1,\dots, n$, and $A\in [C([0,T] \times \Omega^1\cup\Omega^2)]^{n\times n}$ diagonal, with $A=\diag(a_1,a_2,\dots,a_n)$, where $a_i:[0,T]\times \Omega^1\cup\Omega^2\to\mathbb{R}$, $i=1,\dots, n$.  We assume that there exists a constant $\alpha_{\min}>0$ of uniform parabolicity such that $a_i(t,x)\ge\alpha_{\min}$ for all $i=1,\dots,n$ and $(t,x) \in [0,T] \times \Omega$. For simplicity, we also require that the matrix $1/2\diag(\nabla\cdot \mbf{B})$ is positive semi-definite. Finally, $\mbf{F}: \mathbb{R}^n\to\mathbb{R}^n$ is a vector field satisfying the growth condition
\begin{equation}\label{growth}
|\mbf{F}({\mbf{w}})-\mbf{F}({\mbf{v}})|\le C(1+|{\mbf{w}}|+|{\mbf{v}}|)^{\gamma}|{\mbf{w}}-{\mbf{v}}|,
\end{equation}
for ${\mbf{w}}, {\mbf{v}}\in\mathbb{R}^n$ and $\gamma \ge 0$ constant, where $|\cdot|$ denotes the Euclidean distance on $\mathbb{R}^n$. The admissible values of the constant $\gamma$ will be discussed in detail at various instances in the text.

We impose the initial condition
\begin{equation}\label{modelbc}
  {\mbf{u}}(0,x)={\mbf{u}}_0(x)\text{ on }\{0\}\times\Omega,
\end{equation}
for ${\mbf{u}}_0\in [L^2(\Omega)]^n$. On $[0,T]\times\partial\Omega$, we impose mixed Dirichlet and Neumann boundary conditions as follows. Let $i$ be an index running over the set $\{1,\ldots,n\}$. For each $i$ the boundary $\partial\Omega$ is split into $\partial\Omega=\Gamma^{\text{\rm D}}_i\cup\Gamma^{\text{\rm N}}_i$, with $\Gamma^{\text{\rm D}}_i$ being of positive $(d-1)$-dimensional (Hausdorff) measure. Further, we subdivide $\partial\Omega=\partial\Omega_i^-\cup\partial\Omega_i^+$, where $\partial\Omega_i^-:=\{\mbf{x}\in\partial\Omega: (B_i\mbf{n})(\mbf{x})<0\}$ and $\partial\Omega_i^+=\partial\Omega\backslash\partial\Omega_i^-$ are the \emph{inflow} and \emph{outflow} parts of the boundary $\partial\Omega$ for the $i$-th equation. Finally, we assign
\begin{equation}\label{bc}
\begin{aligned}
u_i &= g^{\text{\rm D}}_i \quad\text{on}\ \Gamma^{\text{\rm D}}_i, \\
a_i\nabla u_i\cdot \mbf{n} &= g^{\text{\rm N}}_i \quad\text{on}\ \Gamma^{\text{\rm N}}_i\cap\partial\Omega_i^+,\\
\big(a_i\nabla u_i - B_i^{\mathsf T}u_i\big)\cdot \mbf{n} &= g^{\text{\rm N}}_i \quad\text{on}\ \Gamma^{\text{\rm N}}_i\cap\partial\Omega_i^-,
\end{aligned}
\end{equation}
for Dirichlet and Neumann data $g^{\text{\rm D}}_i\in H^{1/2}(\Gamma^{\text{\rm D}}_i)$, $g^{\text{\rm N}}_i\in L^2(\Gamma^{\text{\rm N}}_i)$, respectively; here and in what follows $\mbf{n}$ denotes the unit outward normal vector to $\partial\Omega$.  We denote by $\chi_i^-:\partial\Omega_i^-\to\mathbb{R}$ the characteristic function of $\partial\Omega_i^-$.  Upon defining $\mathcal{X}^-:=\diag(\chi_1^{-},\dots,\chi_n^{-})$,  and  $\mathcal{X}^+:=I-\mathcal{X}^-$, we can write (\ref{bc}) collectively as
\begin{equation}\label{bcvec}
{\mbf{u}} = \mbf{g}_{\text{\rm D}} \quad\text{on}\ \Gamma_{\text{\rm D}}\quad\text{and}\quad \big(A\nabla {\mbf{u}}-\mathcal{X}^-U\mbf{B}\big)\mbf{n} = \mbf{g}_{\text{\rm N}} \quad\text{on}\ \Gamma_{\text{\rm N}},
\end{equation}
where
$\mbf{g}_{\text{\rm D}}:=(g^{\text{\rm D}}_1,\dots,g^{\text{\rm D}}_n)^{{\mathsf T}}$ and $\mbf{g}_{\text{\rm N}}:=(g^{\text{\rm N}}_1,\dots,g^{\text{\rm N}}_n)^{{\mathsf T}}$.

The model problem is completed imposing, across $\Gamma_{\inter}$, the fluxes described in Section~\ref{sec:tras}. In view of~\eqref{eq:KK2_final} and \eqref{eq:tiltedaverage}, we define the friction coefficients and weights $r_i,\upsilon_i^{1,2}:\Gamma_{\inter}\rightarrow [0,1]$ and the {\em permeabilities} $\mbf{p}_i:\mathbb{R}^{2n}\to\mathbb{R}^n$, $i=1,\dots,n$, as functions of the traces of ${\mbf{u}}$ from both sides of the interface.  The interface conditions read
\[
\begin{aligned}
(a_i\nabla u_i-u_iB_i^{\mathsf T})\cdot \mbf{n}|_{\Omega^1} &= \mbf{p}_i({\mbf{u}}^1,{\mbf{u}}^2) \cdot
({\mbf{u}}^2-{\mbf{u}}^1)
-r_i(\upsilon_i^1 u_i^1+\upsilon_i^2u_i^2)(B_i\mbf{n})|_{\Omega^1},\quad\text{on}\ \Gamma_{\inter},\\
(a_i\nabla u_i-u_iB_i^{\mathsf T})\cdot \mbf{n}|_{\Omega^2} &= \mbf{p}_i({\mbf{u}}^1,{\mbf{u}}^2)\cdot
({\mbf{u}}^1-{\mbf{u}}^2)
-r_i(\upsilon_i^1 u_i^1+\upsilon_i^2u_i^2)(B_i
\mbf{n})|_{\Omega^2},\quad\text{on}\ \Gamma_{\inter},
\end{aligned}
\]
where ${\mbf{u}}^j:={\mbf{u}}|_{\bar{\Omega}_j\cap\Gamma_{\inter}}$, $j=1,2$. Introducing,
$
\Upsilon^{j} =\textrm{diag}(\upsilon^{j}_1,\dots,\upsilon^{j}_n), \quad  j \in \{1,2\}$,
$
{\rm R} =\textrm{diag}(r_1,\dots,r_n)$, and
$\mbf{P}({\mbf{u}}) =(\mbf{p}_1({\mbf{u}}^1,{\mbf{u}}^2),\dots,\mbf{p}_n({\mbf{u}}^1,{\mbf{u}}^2))^{{\mathsf T}}$,
the interface conditions can be written in vector notation as
\begin{equation}\label{transmission}
\begin{aligned}
(A\nabla {\mbf{u}}-U\mbf{B}) \, \mbf{n}|_{\Omega^1} &= \mbf{P}({\mbf{u}}) ({\mbf{u}}^2-{\mbf{u}}^1)
-{\rm R}(\Upsilon^1 U^1+\Upsilon^2U^2)(\mbf{B}\mbf{n})|_{\Omega^1},\quad\text{on}\ \Gamma_{\inter},\\
(A\nabla {\mbf{u}}-U\mbf{B}) \, \mbf{n}|_{\Omega^2} &= \mbf{P}({\mbf{u}}) ({\mbf{u}}^1-{\mbf{u}}^2)
-{\rm R}(\Upsilon^1 U^1+\Upsilon^2U^2)(\mbf{B}\mbf{n})|_{\Omega^2},\quad\text{on}\ \Gamma_{\inter}.
\end{aligned}
\end{equation}

We make the following assumptions on the weights and permeabilities in accordance with Section~\ref{sec:tras}. For every $i=1,\dots,n$, the weights $\upsilon_i^{1,2}$ satisfy, for any  $\mbf{x}\in\Gamma_{\inter}$,
\begin{equation}\label{eq:weights_conditions}
\upsilon_i^1(\mbf{x})+\upsilon_i^2(\mbf{x})=1,\qquad
\left\{\begin{array}{l}
\vspace{1mm}
\upsilon_i^1(\mbf{x})\ge \upsilon_i^2(\mbf{x})\quad\textrm{if }(B_i\mbf{n}|_{\partial\Omega^1})(\mbf{x})\ge 0,\\
\upsilon_i^1(\mbf{x})< \upsilon_i^2(\mbf{x})\quad\textrm{otherwise.}
\end{array}\right.
\end{equation}

We let $\tilde{\mbf{p}}:\mathbb{R}^{2n}\rightarrow \mathbb{R}^n$ denote the function describing the diffusive flux across the interface $\Gamma_{\inter}$, that is
\begin{equation}\label{eq:pjump}
\tilde{\mbf{p}}({\mbf{x}}^1,{\mbf{x}}^2)
:=\mbf{P}({\mbf{x}}^1,{\mbf{x}}^2)
({\mbf{x}}^1-{\mbf{x}}^2)\quad\forall\, {\mbf{x}}^1,{\mbf{x}}^2\in\mathbb{R}^n,
\end{equation}
and assume that $\tilde{\mbf{p}}\in C^{1,1}(\mathbb{R}^{2n})$ and that its Jacobian $\tilde{\mbf{p}}'$ is bounded.

Throughout this work, we assume that the problem given by~\eqref{pde}, \eqref{modelbc}, \eqref{bcvec}, and \eqref{transmission} has a unique solution that remains bounded up to, and including, the final time $T$.

\section{Space discretisation by the discontinuous Galerkin method}\label{dg_heat}

\subsection{Finite element spaces}
Let $\mathcal{T}$ be a shape-regular and locally quasi-uniform subdivision of $\Omega$ into disjoint open elements $ \kappa \in \mathcal{T}$, such that $\Gamma_{\inter}\subset \cup_{ \kappa \in \mathcal{T}}\partial\kappa=:\Gamma$, the \emph{skeleton}. Further we decompose $\Gamma$ into three disjoint subsets $\Gamma=\partial\Omega\cup\Gamma_{\text{\rm int}}\cup\Gamma_{\inter}$, where $\Gamma_{\text{\rm int}}:=\Gamma\backslash(\partial\Omega\cup\Gamma_{\inter})$. We assume that the subdivision $\mathcal{T}$ is constructed via mappings $F_{\kappa}$, where $F_{\kappa}:\hat{\kappa}\to\kappa$ are smooth maps with non-singular Jacobian, and $\hat{\kappa}$ is the reference $d$-dimensional simplex or the reference $d$-dimensional (hyper)cube. It is assumed that the union of the closures of the elements $ \kappa \in \mathcal{T}$ forms a covering of the closure of $\Omega$; i.e., $\bar{\Omega}=\cup_{ \kappa \in \mathcal{T}}\bar{\kappa}$.

For $m \in \mathbb N$ we denote by $\mathbb{P}_m(\hat{\kappa})$ the set of polynomials of total degree at most $m$ if $\hat{\kappa}$ is the reference simplex, and the set of all tensor-product polynomials on $\hat{\kappa}$ of degree $k$ in each variable, if $\hat{\kappa}$ is the reference hypercube.  Let $m_{\kappa}\in\mathbb{N}$ be given for each $ \kappa \in \mathcal{T}$.  We consider the $hp$-discontinuous finite element space
\begin{equation}\label{eq:FEM-spc}
{\rm V}_h:=\{v\in L^2(\Omega):v|_{\kappa}\circ F_{\kappa}
   \in \mathbb{P}_{m_{\kappa}}(\hat{\kappa}),\, \kappa \in \mathcal{T}\},
\end{equation}
and set $\mathbb{V}_h:=[{\rm V}_h]^n$.

Next, we introduce relevant trace operators. Let $\kappa^+$, $\kappa^-$ be two elements sharing an edge $e:=\partial\kappa^+\cap\partial\kappa^-\subset\Gamma_{\text{\rm int}}\cup \Gamma_{\inter}$. Denote the outward normal unit vectors on $e$ of $\partial\kappa^+$ and $\partial\kappa^-$ by $\mbf{n}^+$ and $\mbf{n}^-$, respectively. For functions $\mbf{q}:\Omega\to\mathbb{R}^n$ and $\mbf{Q}:\Omega\to\mathbb{R}^{n\times d}$ that may be discontinuous across $\Gamma$, we define the following quantities: for $\mbf{q}^+:=\mbf{q}|_{\kappa^+}$, $\mbf{q}^-:=\mbf{q}|_{\kappa^-}$ and $\mbf{Q}^+:=\mbf{Q}|_{\kappa^+}$, $\mbf{Q}^-:=\mbf{Q}|_{\kappa^-}$ on the restriction to $e$, we set
$ \mean{\mbf{q}}:=\frac{1}{2}(\mbf{q}^+ + \mbf{q}^-)$, $ \mean{\mbf{Q}}:=\frac{1}{2}(\mbf{Q}^+ + \mbf{Q}^-)$,
and
$ \jumptwo{\mbf{q}}:=\mbf{q}^+\otimes \mbf{n}^++\mbf{q}^-\otimes \mbf{n}^-$,
$ \jump{\mbf{Q}}:=\mbf{Q}^+ \mbf{n}^++\mbf{Q}^- \mbf{n}^-$,
where $\otimes$ denotes the standard tensor product operator, with $\mbf{q}\otimes\mbf{w}=\mbf{q}\mbf{w}^{{\mathsf T}}$. If $e\in \partial\kappa\cap\Gamma_{\partial}$, these definitions are modified as follows: $\mean{\mbf{q}}:=\mbf{q}^+ ,\ \mean{\mbf{Q}}:=\mbf{Q}^+$ and $\jumptwo{\mbf{q}}:=q^+\otimes\mbf{n},\ \jump{\mbf{Q}}:=\mbf{Q}^+\mbf{n}$.

Further, we introduce the mesh quantities $\mbf{\tt h}:\Omega\to \mathbb{R}$, ${\tt m}:\Omega\to \mathbb{R}$ with $\mbf{\tt h}(x)= \diam{\kappa}$, ${\tt m}(x)= m_{\kappa}$, if $x\in \kappa$, and the averaged values $\mbf{\tt h}(x)= \mean{\mbf{\tt h}}$, ${\tt m}(x)= \mean{{\tt m}}$, if $x\in \Gamma$.  Finally, we define $h_{\max}:=\max_{x\in{\Omega}}\mbf{\tt h}$ and $h_{{\rm min}}:=\min_{x\in{\Omega}}\mbf{\tt h}$.

We shall assume the existence of a constant $C_A\ge1$ independent of $\mathcal{T}$ such that, on any face  that is \emph{not} contained in $\Gamma_{\inter}$, given the two elements $\kappa$, $\kappa'$ sharing that face, the diffusion matrix $A$ satisfies
\begin{equation}\label{blv_A}
C_A^{-1} \le \big\| A\|_{\infty,\kappa}\|A^{-1}\big\|_{\infty,\kappa'}\le C_A.
\end{equation}
This assumption can be removed using the ideas from \cite{MR2491426}, but we refrain from doing so here for simplicity of the presentation.

The next result is a modification of the classical trace estimate for functions in $H^1(\Omega^1\cup\Omega^2)+{\rm V}_h$; see~\cite{BuffaOrtner} for similar results.

\begin{lemma}\label{broken_trace
Assume that the mesh $\mathcal{T}$ is both shape-regular and locally quasi-uniform. Then for $v\in H^1(\Omega^1\cup\Omega^2)+{\rm V}_h$, the following trace estimate holds:
\begin{equation}\label{trace_estimate}
\sum_{j=1}^2\ltwo{v|_{\Omega^j}}{\Gamma_{\inter}}^2\le c_1 \epsilon \Big(\sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d}\ltwo{\nabla v}{\kappa}^2 + \ltwo{\mbf{\tt h}^{-1/2}\jump{v}}{\Gamma_{\dint}}^2\Big)+ c_2 \epsilon^{-1}\ltwo{v}{}^2,
\end{equation}
for any $\epsilon>h_{\max}$ and for some constants $c_1>0$ and $c_2>0$, depending only on the shape-regularity of the mesh and on the domain $\Omega$.
\end{lemma}
\begin{proof}
We use the decomposition of $v\in H^1(\Omega^1\cup\Omega^2)+{\rm V}_h$ into a conforming part $v_c$ and a non-conforming part $v_d:=v-v_c\in {\rm V}_h$. This decomposition is described in \cite{karakashian_pascal,karakashian_pascal_convergence} for functions in ${\rm V}_h$; the extension to $H^1(\Omega^1\cup\Omega^2)+{\rm V}_h$ follows by taking $v_c\in H^1(\Omega^1\cup\Omega^2)$. Using Theorem 2.1(iii) from \cite{karakashian_pascal_convergence}, there exists a $v_c^i\in  H^1(\Omega^1\cup\Omega^2)$, $i=1,2$, such that
\begin{equation}\label{recovery}
\sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d}\ltwo{D^{\alpha}(v-v_c^i)}{\kappa\cap\Omega^i}^2
\le
C \ltwo{\mbf{\tt h}^{1/2-|\alpha |}\jump{v}}{\Gamma_{\dint}\cap\Omega^i}^2,
\end{equation}
where $D^{\alpha}$ is the differentiation operator for a multi-index $\alpha$ with $|\alpha|=0,1$. Hence, (\ref{recovery}) implies
\begin{equation}\label{recovery_transmission}
\sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d}\ltwo{D^{\alpha}(v-v_c)}{\kappa}^2
\le
 C \ltwo{\mbf{\tt h}^{1/2-|\alpha |}\jump{v}}{\Gamma_{\dint}}^2,
\end{equation}
for $(v_c)|_{\Omega^j}:=v_c^j$, $j=1,2$.

The triangle inequality implies
\begin{equation}\label{triangle_in}
\sum_{j=1}^2\ltwo{v|_{\Omega^j}}{\Gamma_{\inter}}^2
\le
2\sum_{j=1}^2\big(\ltwo{v_c|_{\Omega^j}}{\Gamma_{\inter}}^2+\ltwo{(v-v_c)|_{\Omega^j}}{\Gamma_{\inter}}^2\big).
\end{equation}
To bound the first term on the right-hand side of (\ref{triangle_in}), we note that $v_c\in H^1(\Omega^1\cup\Omega^2)$, giving
\begin{equation}\label{conforming_part}
\begin{aligned}
\sum_{j=1}^2\ltwo{v_c|_{\Omega^j}}{\Gamma_{\inter}}^2
&\le
C \sum_{j=1}^2 \big( \ltwo{v_c}{\Omega^j}^2 + \ltwo{v_c}{\Omega^j}\ltwo{\nabla v_c}{\Omega^j}\big)^{1/2}
\le
C \big( \ltwo{v_c}{}^2 + \ltwo{v_c}{}\ltwo{\nabla v_c}{}\big)^{1/2}\\
&\le
C\big(  \epsilon^{-1}\ltwo{v_c}{}^2 + \epsilon\ltwo{\nabla v_c}{}^2\big)^{1/2},
\end{aligned}
\end{equation}
for any $\epsilon>0$ sufficiently small.  To further bound the right-hand side of (\ref{conforming_part}), we use the triangle inequality for each term, viz.,
\begin{equation}\label{triangle_two}
\ltwo{v_c}{}\le \ltwo{v-v_c}{} +\ltwo{v}{}
\quad\text{and}\quad
\ltwo{\nabla v_c}{}^2\le 2\sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d}\big(\ltwo{\nabla(v-v_c)}{\kappa}^2 +\ltwo{\nabla v}{\kappa}^2\big),
\end{equation}
in conjunction with (\ref{recovery_transmission}), to arrive at
\begin{equation}\label{conforming_part_two}
\begin{aligned}
\sum_{j=1}^2\ltwo{v_c|_{\Omega^j}}{\Gamma_{\inter}}^2
&\le
C\Big(  \epsilon^{-1} \ltwo{\mbf{\tt h}^{1/2}\jump{v}}{\Gamma_{\dint}}^2 +  \epsilon^{-1} \ltwo{v}{}^2  +  \epsilon \ltwo{\mbf{\tt h}^{-1/2}\jump{v}}{\Gamma_{\dint}}^2 +\epsilon \sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d}\ltwo{\nabla v}{\kappa}^2\Big)^{1/2}\\
&\le
C\Big(   \epsilon^{-1} \ltwo{v}{}^2  +  \epsilon
\ltwo{\mbf{\tt h}^{-1/2}\jump{v}}{\Gamma_{\dint}}^2 +\epsilon \sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d}\ltwo{\nabla v}{\kappa}^2\Big)^{1/2},
\end{aligned}
\end{equation}
using the assumption that $\epsilon>h_{\max}$.

To bound the second term on the right-hand side of (\ref{triangle_in}), we use the trace estimate on each element in conjunction with (\ref{recovery_transmission}), viz.,
\begin{equation}\label{nonconforming_part}
\begin{aligned}
\sum_{j=1}^2\ltwo{(v-v_c)|_{\Omega^j}}{\Gamma_{\inter}}^2
&\le
C\!\!\!\!\!\!\!\!\sum_{ \kappa \in \mathcal{T}: \bar{\kappa}\cap\Gamma_{\inter}\ne \emptyset}\!\!\!\!\!\!\!\!\big(     \ltwo{\mbf{\tt h}^{-1/2}(v-v_c)}{\kappa}^2  +   \ltwo{\mbf{\tt h}^{1/2}\nabla(v-v_c)}{\kappa}^2\big)
\le
C  \ltwo{\jump{v}}{\Gamma_{\dint}}^2.
\end{aligned}
\end{equation}
Noting the assumption $\epsilon>h_{\max}$, the use of (\ref{conforming_part_two}) and (\ref{nonconforming_part}) in (\ref{triangle_in}) concludes the proof.
\end{proof}

\subsection{Space discretization}

The discretization of the space variables will be based on a discontinuous Galerkin method of interior penalty type for the diffusion part and of upwind type for the advection. Special care has to be given to the incorporation of the interface conditions.

More specifically, we shall introduce a \emph{interior penalty discontinuous Galerkin} (IPDG, for short)  discretisation of the advection-diffusion operator
\begin{equation}\label{eq:sec-system}
-\nabla\cdot (A \nabla{\mbf{w}}-W\mbf{B})
,
\end{equation}
where $W:=\diag(w_1,w_2,\dots, w_n)$, for ${\mbf{w}}:=(w_1,w_2,\dots, w_n)^{{\mathsf T}}$.

To this end, we define $B({\mbf{w}}_h,{\mbf{v}}_h)$ to be the IPDG bilinear form
\begin{equation}\label{eq:bilinear}
\begin{aligned}
&\sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d}\int_{\kappa} (A\nabla{\mbf{w}}_h-W_h\mbf{B}): \nabla {\mbf{v}}_h
+\int_{\Gamma_{\inter}}\Big(
\mean{W_h\mbf{B}}+\mathcal{B}_{{\mathcal {I}}}\jumptwo{{\mbf{w}}_h}
\Big):\jumptwo{{\mbf{v}}_h}\\
&-\int_{\Gamma_{\dint}}\Big(\mean{A\nabla {\mbf{w}}_h-W_h\mbf{B}}:\jumptwo{{\mbf{v}}_h}+\mean{A\nabla {\mbf{v}}_h}:\jumptwo{{\mbf{w}}_h}
-(\Sigma+\mathcal{B})\jumptwo{{\mbf{w}}_h}:\jumptwo{{\mbf{v}}_h}\Big)\\
&-\int_{\Gamma_{\text{\rm D}}}\Big((A\nabla {\mbf{w}}_h-\mathcal{X}^+W_h\mbf{B}):({\mbf{v}}_h\otimes \mbf{n})+(A\nabla {\mbf{v}}_h):({\mbf{w}}_h\otimes \mbf{n})
-\Sigma{\mbf{w}}_h\cdot{\mbf{v}}_h\Big)\\
&+\int_{\Gamma_{\text{\rm N}}}(\mathcal{X}^+W_h\mbf{B}):({\mbf{v}}_h\otimes \mbf{n}),
\end{aligned}
\end{equation}
and define
\begin{equation}\label{eq:nonlinear}
N({\mbf{w}}_h,{\mbf{v}}_h):=
\int_{\Gamma_{\inter}}\left(\mbf{P}({\mbf{w}}_h)\jumptwo{{\mbf{w}}_h}
-({\rm I}-{\rm R})
\left(
\mean{W_h\mbf{B}}+\mathcal{B}_{{\mathcal {I}}}\jumptwo{{\mbf{w}}_h}
\right)
\right):\jumptwo{{\mbf{v}}_h}.
\end{equation}
Here, $\Sigma:=C_{\sigma}A{\tt m}^2\mbf{\tt h}^{-1}$ denotes the \emph{discontinuity-penalization parameter} matrix with $C_{\sigma}>1$ constant.  Furthermore, $\mathcal{B}:=\frac{1}{2}\diag(|B_1\cdot \mbf{n}|,\dots,|B_n\cdot \mbf{n}|)$ and $\mathcal{B}_{{\mathcal {I}}}:=(\Upsilon^1-\frac{1}{2}{\rm I})\mbf{B}\mbf{n}|_{\Omega^1}= (\Upsilon^2-\frac{1}{2}{\rm I})\mbf{B}\mbf{n}|_{\Omega^2}$ is diagonal with non negative entries.

\begin{remark} We comment on the interface terms. The diffusion term, appearing in $N$, is simply given by $\int_{\Gamma_{\inter}}\mbf{P}({\mbf{w}})\jumptwo{{\mbf{w}}}:\jumptwo{{\mbf{v}}_h}$. This resembles the typical jump stabilisation term with the permeability coefficient replacing the discontinuity-penalization parameter, rendering its implementation within a dG computer code straightforward.

To ensure the coercivity of $B$, the advective interface term has been split as ${\rm R} ={\rm I}-({\rm I}-{\rm R})$, resulting into contributions in both $B$ and $N$. Indeed, the advective interface contribution in $B$ can be recast using the weighted mean $\mean{W_h\mbf{B}}^\upsilon:=\Upsilon^1W_h\mbf{B}|_{\Omega^1} + \Upsilon^2W_h\mbf{B}|_{\Omega^2}$, so that
\begin{equation}\label{eq:weighted_average}
\mean{W_h\mbf{B}}^\upsilon:\jumptwo{{\mbf{v}}_h}=\left(
\mean{W_h\mbf{B}}+\mathcal{B}_{{\mathcal {I}}}\jumptwo{{\mbf{w}}_h}
\right):\jumptwo{{\mbf{v}}_h},
\end{equation}
thereby resembling the typical dG upwinding for linear advection problems and, hence, ensuring the coercivity of $B$.
\end{remark}
\begin{remark}
In this setting, $\Gamma_{\text{\rm N}}$ can have non-trivial intersection with both $\partial\Omega^-_i$ and $\partial\Omega^+_i$, thereby extending the discontinuous Galerkin method proposed in \cite{houston-etal:01}.
\end{remark}

\section{Elliptic projection error}\label{err_ell_problem}
The a priori error analysis is based on a (non-standard) elliptic projection inspired by a construction of Douglas and Dupont for the treatment of nonlinear boundary conditions in the context of conforming finite element methods \cite{douglas_dupont}.

\begin{definition}
For each $t\in [0,T]$ we define the \emph{elliptic projection} ${\mbf{w}}_h\in \mathbb{V}_h$ to be the solution of the problem: find ${\mbf{w}}_h\equiv{\mbf{w}}_h(t)\in \mathbb{V}_h$, such that
\begin{equation}\label{ellipitc_projection}
{B}({\mbf{u}}-{\mbf{w}}_h,{\mbf{v}}_h)+\lambda \ltwoin{{\mbf{u}}-{\mbf{w}}_h}{{\mbf{v}}_h}
+N({\mbf{u}},{\mbf{v}}_h)-N({\mbf{w}}_h,{\mbf{v}}_h)
=
0\quad \forall {\mbf{v}}_h\in \mathbb{V}_h,
\end{equation}
for some fixed $\lambda>0$, and ${\mbf{u}}$ denoting the exact solution.
\end{definition}

The constant $\lambda>0$ in the definition above will be chosen large enough to ensure the uniqueness of the projection ${\mbf{w}}_h$, cf. Lemma \ref{ell_proj_lemma} below.

Next, denoting by $\mathbb{S}^s:=\mathbb{H}^s+\mathbb{V}_h$, $s\in\mathbb{R}$,
we define the dG-norm on $\mathbb{S}^1$
\begin{equation}
\label{eq:dGnorm}
\begin{aligned}
\ndg{{\mbf{w}}}:=\bigg(&\sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d}\big(\ltwo{\sqrt{A}\nabla{\mbf{w}}}{\kappa}^2+\frac{1}{2}\ltwo{\sqrt{\diag(\nabla\cdot \mbf{B})}{\mbf{w}}}{\kappa}^2\big)
+\ltwo{\sqrt{\Sigma}\jumptwo{{\mbf{w}}}}{\Gamma_{\text{\rm D}}\cup\Gamma_{\text{\rm int}}}^2
\\
&+\ltwo{\sqrt{\mathcal{B}}\jumptwo{{\mbf{w}}}}{\Gamma\backslash\Gamma_{\inter}}^2
+\ltwo{
\sqrt{\mathcal{B}_{{\mathcal {I}}}}
\jumptwo{{\mbf{w}}}}{\Gamma_{\inter}}^2
\bigg)^{1/2},
\end{aligned}
\end{equation}
where $
\ltwo{\mbf{Q}}{\kappa}^2 := \int_\kappa \sum_{i=1}^n|Q_i(x)|^2 \, {\rm d}x$,
denotes the Frobenius norm whenever $\mbf{Q}$ is a $n\times d$ tensor. We assume that~\eqref{eq:dGnorm} is a norm. This is satisfied when standard assumptions on the solution in conjunction with the boundary conditions hold on each subdomain, e.g., $\Gamma_{\text{\rm D}}\cap\partial\Omega^j$ has positive $(d-1)$-dimensional (Hausdorff) measure for $j=1,2$. If the interface manifold $\Gamma_{\inter}$ is not characteristic to the advection field, such hypotheses can be further relaxed.

For the remainder of this work, we shall make the simplifying assumption that $\mbf{B}$ is such that:
\begin{equation}\label{wind_assumption}
 B_i\cdot\nabla (v_h)_i\in {\rm V}_h,\quad \text{ for }\ i=1,\dots,n,
\end{equation}
for any function ${\mbf{v}}_h:=((v_h)_1,\dots,(v_h)_n)^{{\mathsf T}}\in \mathbb{V}_h$. We note, however, that this appears not to be a genuine limitation in the arguments presented below: ideas on how to circumvent this assumption have been presented, e.g., in \cite{houston-etal:01, ayuso_marini}, for the case of scalar linear advection-diffusion problems.

The next two results show the coercivity and the continuity of the bilinear form $B(\cdot,\cdot)$. Their proofs follow straightforward variations of well-known arguments (see, e.g., \cite{arnold,houston-etal:01}) and are, therefore, omitted for brevity.

\begin{lemma}\label{coercivity}
There exists a positive constant $C_{coer}\in\mathbb{R}$, such that, for ${\mbf{v}}_h\in \mathbb{V}_h$,
\[
B({\mbf{v}}_h,{\mbf{v}}_h)\ge C_{coer}\ndg{{\mbf{v}}_h}^2.
\]
\end{lemma}
\begin{lemma}\label{continuity}
Let $\Pi:[L^2(\Omega)]^n\to \mathbb{V}_h$ denote the $L^2$-orthogonal projection onto $\mathbb{V}_h$.
For any ${\mbf{w}}\in \mathbb{H}^{s}$, $s>3/2$ and ${\mbf{v}}_h\in \mathbb{V}_h$ we have
\[
|{B}(\mbf{\eta},{\mbf{v}}_h)|
\le   C_{cont} \ndg{\mbf{\eta}}_{\mathcal{B}}
\ndg{{\mbf{v}}_h},
\]
where $\mbf{\eta}:={\mbf{w}}-\Pi {\mbf{w}}$ and
\begin{equation}\label{eq:Aeta}
\ndg{\mbf{\eta}}^2_{\mathcal{B}}:=
\ndg{\mbf{\eta}}^2
+\ltwo{{\Sigma}^{-1/2}\mean{A\nabla \mbf{\eta}}}{\Gamma_{\text{\rm D}}\cup\Gamma_{\text{\rm int}}}^2
+\ltwo{\sqrt{\mathcal{B}}\mean{\mbf{\eta}}}{\Gamma}^2
+\ltwo{\sqrt{\mathcal{B}_{{\mathcal {I}}}}\mean{\mbf{\eta}}}{\Gamma_{\inter}}^2.
\end{equation}
\end{lemma}

The next result establishes the well-posedness of the
problem~\eqref{ellipitc_projection} and relevant approximation properties.

\begin{lemma}\label{ell_proj_lemma}
Assume that ${\mbf{u}}\in \mathbb{H}^s$, $s>3/2$
for all $t\in(0,T]$.
For $\lambda>0$ sufficiently large and for $h_{\max}$ sufficiently small, the variational problem (\ref{ellipitc_projection}) has a unique solution ${\mbf{w}}_h\in \mathbb{V}_h$ for each $t\in(0,T]$. Moreover, the following bound holds:
\begin{equation}
\label{ell_proj_bound}
C_{coer}\ndg{\mbf{\rho}}^2+ \lambda \ltwo{\mbf{\rho}}{}^2
\le \ndg{\mbf{\eta}}_{\mathcal{B},\lambda}^2,
\end{equation}
and, if also ${\mbf{u}}_t\in \mathbb{H}^s$, then
\begin{equation}\label{ell_proj_bound_t}
C_{coer}\ndg{\mbf{\rho}_t}^2+ \lambda \ltwo{\mbf{\rho}_t}{}^2 \le \ndg{\mbf{\eta}_t}_{\mathcal{B},\lambda}^2+\ndg{\mbf{\eta}}_{\mathcal{B},\lambda}^2
\end{equation}
where $\mbf{\rho}:={\mbf{u}}-{\mbf{w}}_h$, $\mbf{\eta}:={\mbf{u}}-\Pi{\mbf{u}}$, and
\[
\ndg{\mbf{\eta}}_{\mathcal{B},\lambda}^2:=
C_{c}\ndg{\mbf{\eta}}^2_{\mathcal{B}}+7\lambda \ltwo{\mbf{\eta}}{}^2,
\]
with $C_{c}:=(4C_{cont}^2+3C_{coer}^2)/C_{coer}$.
\end{lemma}

{\bf Proof.}
Well-posedness of (\ref{ellipitc_projection}) is established by proving that the associated mapping is strongly monotone on $\mathbb{V}_h$.  Using the assumption that $\tilde{\mbf{p}}\in C^{0,1}(\mathbb{R}^{2n})$
, we  get
\begin{equation}\label{eq:N_bound}
\begin{aligned}
|N({\mbf{v}},{\mbf{z}})-N({\mbf{w}},{\mbf{z}})|&\le \int_{\Gamma_{\inter}} (|\tilde{\mbf{p}} ({\mbf{v}})-\tilde{\mbf{p}} ({\mbf{w}})|+C_\mathcal{B} |{\mbf{v}}-{\mbf{w}})|)\,|\jumptwo{{\mbf{z}}}|\\
&\le \int_{\Gamma_{\inter}} \left(C^p|{\mbf{v}}-{\mbf{w}}|+C_\mathcal{B}|{\mbf{v}}-{\mbf{w}}|\right) \,|\jumptwo{\mbf{z}}|\\
&\le C^{p}_{\mathcal{B}} \sum_{j=1}^2(\ltwo{({\mbf{v}}-{\mbf{w}})|_{\Omega^j}}{\Gamma_{\inter}}^2 +\ltwo{{\mbf{z}}|_{\Omega^j}}{\Gamma_{\inter}}^2),
\end{aligned}
\end{equation}
where $C^p$ is a Lipschitz constant for the function $\tilde{\mbf{p}}$,  $C_{\mathcal{B}}>0$ is a  constant proportional to $\displaystyle{\max_{i=1,\dots,n}\ltwo{B_i\mbf{n}}{\infty,\Gamma_{\inter}}}$, and $C^{p}_{\mathcal{B}}>0$ is a constant depending on both $C^p$ and $C_{\mathcal{B}}$.

This, in conjunction with the coercivity of the bilinear form $B$ provided by Lemma~\ref{coercivity}, gives
\[
\begin{aligned}
&B({\mbf{v}}_h-{\mbf{w}}_h,{\mbf{v}}_h-{\mbf{w}}_h)+\lambda \ltwoin{{\mbf{v}}_h-{\mbf{w}}_h}{{\mbf{v}}_h-{\mbf{w}}_h}
+N({\mbf{v}}_h,{\mbf{v}}_h-{\mbf{w}}_h)-N({\mbf{w}}_h,{\mbf{v}}_h-{\mbf{w}}_h)\\
&\ge C_{coer}\ndg{{\mbf{v}}_h-{\mbf{w}}_h}^2
+\lambda \ltwo{{\mbf{v}}_h-{\mbf{w}}_h}{}^2-C^p_{\mathcal{B}}\sum_{j=1}^2\ltwo{({\mbf{v}}_h-{\mbf{w}}_h)|_{\Omega^j}}{\Gamma_{\inter}}^2\\
&\ge \frac{C_{coer}}{2}\ndg{{\mbf{v}}_h-{\mbf{w}}_h}^2
+\big(\lambda
-\frac{2 C^p_\mathcal{B} c_1c_2}{C_{coer}\alpha_{\min}}\big)\ltwo{{\mbf{v}}_h-{\mbf{w}}_h}{}^2,
\end{aligned}
\]
the last inequality owning to the trace estimate~\eqref{trace_estimate} with $\epsilon= C_{coer}\alpha_{\min}/(2C^p_{\mathcal{B}}c_1)$. Hence strong monotonicity is ensured as soon as $\lambda>2 C^p_\mathcal{B}c_1c_2/C_{coer}\alpha_{\min}$.

To show~\eqref{ell_proj_bound} and~\eqref{ell_proj_bound_t}, we split the error ${\mbf{u}}-{\mbf{w}}_h=\mbf{\eta}+\mbf{\xi}$, with $\mbf{\eta}:={\mbf{u}}-\Pi{\mbf{u}}$ and $\mbf{\xi}:=\Pi{\mbf{u}}-{\mbf{w}}_h$.  Then, setting ${\mbf{v}}_h=\mbf{\xi}$ in (\ref{ellipitc_projection}), we deduce
\begin{equation}\label{error_eq_ell_proj}
B(\mbf{\xi},\mbf{\xi})+\lambda \ltwoin{\mbf{\xi}}{\mbf{\xi}}= -{B}(\mbf{\eta},\mbf{\xi})-\lambda \ltwoin{\mbf{\eta}}{\mbf{\xi}}-N({\mbf{u}},\mbf{\xi})+N({\mbf{w}}_h,\mbf{\xi}).
\end{equation}
Using the coercivity and the continuity of $B$ along with the bound~\eqref{eq:N_bound} in (\ref{error_eq_ell_proj}) gives
\begin{equation}\label{error_in_elliptic}
 \frac{3}{4}C_{coer}\ndg{\mbf{\xi}}^2+ \frac{3}{4}\lambda \ltwo{\mbf{\xi}}{}^2\le\frac{C_{cont}^2}{C_{coer}}\ndg{\mbf{\eta}}^2_{\mathcal{B}}+\lambda \ltwo{\mbf{\eta}}{}^2+C^p_\mathcal{B}\sum_{j=1}^2\big(\ltwo{\mbf{\eta}|_{\Omega^j}}{\Gamma_{\inter}}^2+\ltwo{\mbf{\xi}|_{\Omega^j}}{\Gamma_{\inter}}^2\big).
\end{equation}
Using (\ref{trace_estimate}) with
$\epsilon= C_{coer}\alpha_{\min}/(4C^p_\mathcal{B}c_1)$ on the last term on the right-hand side of (\ref{error_in_elliptic}), we arrive at
\[
 \frac{1}{2}C_{coer}\ndg{\mbf{\xi}}^2+ \frac{3}{4}\lambda \ltwo{\mbf{\xi}}{}^2\le
C_{c}\ndg{\mbf{\eta}}^2_{\mathcal{B}}+\lambda \ltwo{\mbf{\eta}}{}^2
  +C^p_\mathcal{B}c_2 \epsilon^{-1}(\ltwo{\mbf{\eta}}{}^2
+\ltwo{\mbf{\xi}}{}^2).
\]
Choosing $\lambda> 16(C^p_\mathcal{B})^2c_1c_2/(C_{coer}\alpha_{\min})$, we deduce
\begin{equation}\label{error_in_elliptic_three}
 \frac{1}{2}C_{coer}\ndg{\mbf{\xi}}^2+ \frac{1}{2}\lambda \ltwo{\mbf{\xi}}{}^2\le
C_{c}\ndg{\mbf{\eta}}^2_{\mathcal{B}}+\frac{5}{4}\lambda \ltwo{\mbf{\eta}}{}^2.
\end{equation}
A triangle inequality already gives~\eqref{ell_proj_bound}.

In view of obtaining~\eqref{ell_proj_bound_t}, we differentiate~\eqref{ellipitc_projection} with respect to $t$ and then test with ${\mbf{v}}_h=\mbf{\xi}_t$:
\begin{equation}\label{error_in_elliptic_four}
\begin{aligned}[l]
&B\left(\mbf{\xi}_t,\mbf{\xi}_t\right)+\lambda \ltwoin{\mbf{\xi}_t}{\mbf{\xi}_t}=
-{B}\left(\mbf{\eta}_t,\mbf{\xi}_t\right)-\lambda \ltwoin{\mbf{\eta}_t}{\mbf{\xi}_t}\\
&\quad -\int_{\Gamma_{\inter}}
\frac{d}{d t}
\big(
\tilde{\mbf{p}} ({\mbf{u}})-\tilde{\mbf{p}} ({\mbf{w}}_h)
+({\rm I}-{\rm R})
\left(
\mean{(U-W_h)\mbf{B}}+\mathcal{B}_{{\mathcal {I}}}\jumptwo{\mbf{\rho}}
\right)
\big):\jumptwo{\mbf{\xi}_t}.
\end{aligned}
\end{equation}
Using the assumption that the Jacobian $\tilde{\mbf{p}}'\in C^{0,1}(\mathbb{R}^{2n})$ and is bounded we obtain
\begin{equation}\label{eq:p_t_bound}
\begin{aligned}
\left|\int_{\Gamma_{\inter}}\right.&\left.  \frac{d}{d t} \big( \tilde{\mbf{p}} ({\mbf{u}})-\tilde{\mbf{p}} ({\mbf{w}}_h) +({\rm I}-{\rm R}) \left( \mean{(U-W_h)\mbf{B}}+\mathcal{B}_{{\mathcal {I}}}\jumptwo{\mbf{\rho}} \right) \big):\jumptwo{\mbf{\xi}_t}\right|\\
&\le \int_{\Gamma_{\inter}}\Big(|(\tilde{\mbf{p}}'({\mbf{u}})-\tilde{\mbf{p}}'({\mbf{w}}_h))[{\mbf{u}}_t]| +|\tilde{\mbf{p}}'({\mbf{w}}_h)[\mbf{\rho}_t]|+C_{\mathcal{B}}\sum_{j=1}^2(|\mbf{\rho}_t|_{\Omega^j}|+|\mbf{\rho}|_{\Omega^j}|)\Big) |\jumptwo{\mbf{\xi}_t}| \\
&\le \Big(\sum_{j=1}^2(\tilde{C}^{p}_{\mathcal{B},{\mbf{u}}} \ltwo{\mbf{\rho}|_{\Omega^j}}{\Gamma_{\inter}} +C^{p}_{\mathcal{B}}     \ltwo{\mbf{\rho}_t|_{\Omega^j}}{\Gamma_{\inter}})\Big) \sum_{j=1}^2\ltwo{\mbf{\xi}_t|_{\Omega^j}}{\Gamma_{\inter}}\\
&\le \sum_{j=1}^2\left(\tilde{C}^{p}_{\mathcal{B},{\mbf{u}}} \ltwo{\mbf{\rho}|_{\Omega^j}}{\Gamma_{\inter}}^2 +2C^{p}_{\mathcal{B}} \ltwo{\mbf{\eta}_t|_{\Omega^j}}{\Gamma_{\inter}}^2 +(\tilde{C}^{p}_{\mathcal{B},{\mbf{u}}} +3C^{p}_{\mathcal{B}})\ltwo{\mbf{\xi}_t|_{\Omega^j}}{\Gamma_{\inter}}^2\right) \\
&\le  C^{p}_{\mathcal{B},{\mbf{u}}}\sum_{j=1}^2\left( \ltwo{\mbf{\rho}|_{\Omega^j}}{\Gamma_{\inter}}^2 +\ltwo{\mbf{\eta}_t|_{\Omega^j}}{\Gamma_{\inter}}^2 +\ltwo{\mbf{\xi}_t|_{\Omega^j}}{\Gamma_{\inter}}^2\right),
\end{aligned}
\end{equation}
with $\tilde{C}^{p}_{\mathcal{B},{\mbf{u}}}$, $C^{p}_{\mathcal{B},{\mbf{u}}}>0$ constants depending on $C^{p}_{\mathcal{B}}$ and on $\ltwo{{\mbf{u}}_t}{L^{\infty}([0,T]\times \Omega)}$.  Here and in what follows, square brackets are used to denote the argument of a linear operator.

Applying~\eqref{eq:p_t_bound} to the right-hand side of~\eqref{error_in_elliptic_four} and using the coercivity and continuity of $B$ yields
\begin{equation}\label{error_in_elliptic_five}
\begin{aligned}
 \frac{3}{4}C_{coer}\ndg{\mbf{\xi}_t}^2+ \frac{3}{4}\lambda \ltwo{\mbf{\xi}_t}{}^2&\le
\frac{C_{cont}^2}{C_{coer}} \ndg{\mbf{\eta}_t}^2_{\mathcal{B}}
+\lambda \ltwo{\mbf{\eta}_t}{}^2\\
&+C^{p}_{\mathcal{B},{\mbf{u}}}\sum_{j=1}^2\left(
\ltwo{\mbf{\rho}|_{\Omega^j}}{\Gamma_{\inter}}^2
+\ltwo{\mbf{\eta}_t|_{\Omega^j}}{\Gamma_{\inter}}^2
+\ltwo{\mbf{\xi}_t|_{\Omega^j}}{\Gamma_{\inter}}^2\right).
\end{aligned}
\end{equation}

As before, using the trace estimate (\ref{trace_estimate}) on the last term on the right-hand side of (\ref{error_in_elliptic_five}) with
$\epsilon= C_{coer}\alpha_{\min}/(4C^{p}_{\mathcal{B},{\mbf{u}}}c_1)$, we arrive to
\[
\begin{aligned}
 \frac{1}{2}C_{coer}\ndg{\mbf{\xi}_t}^2+ \frac{1}{2}\lambda \ltwo{\mbf{\xi}_t}{}^2&\le
\frac{4C_{cont}^2+3C_{coer}^2}{C_{coer}} \ndg{\mbf{\eta}_t}^2_{\mathcal{B}}
        +\frac{5}{4}\lambda \ltwo{\mbf{\eta}_t}{}^2
+\frac{1}{4}\ndg{\mbf{\eta}}_{\mathcal{B},\lambda}^2,
\end{aligned}
\]
for any $\lambda> 16(\max\{C^{p}_{\mathcal{B},{\mbf{u}}},C^{p}_{\mathcal{B}}\})^2c_1c_2/(C_{coer}\alpha_{\min})$.  Now~\eqref{ell_proj_bound_t} easily follows by the triangular inequality.
\vspace{-.2cm} \\* \mbox{} \hfill $\Box$

We conclude this section with an $L^2$-error bound of the elliptic projection~(\ref{ellipitc_projection}). This is obtained by an Aubin-Nitsche duality-type argument, inspired by a construction of Douglas and Dupont~\cite{douglas_dupont} for nonlinear boundary conditions.

The interface operator $N$ given in~\eqref{eq:nonlinear} consists of a nonlinear component driven by the function $\tilde{\mbf{p}}({\mbf{w}})=\mbf{P}({\mbf{w}})\jumptwo{{\mbf{w}}}$ and a linear component which we can characterise by  introducing the linear operator $L[{\mbf{w}}]:=-({\rm I}-{\rm R}) \mean{W\mbf{B}}^\upsilon +\mathcal{B}_{{\mathcal {I}}}\jumptwo{{\mbf{w}}}$. We abbreviate $\mathbb{S} := \mathbb{S}^1$ and let $\mathbb{S}^*$ be the dual space of $\mathbb{S}$. It is convenient to momentarily view $N$ as an operator from $\mathbb{S} \to \mathbb{S}^*$, indicated with a calligraphic font:
\begin{equation*}
\mathcal N : \; \mathbb{S} \to \mathbb{S}^*, \; {\mbf{w}} \mapsto \Bigl( {\mbf{v}} \mapsto \int_{\Gamma_{\inter}}\left(\tilde{\mbf{p}}({\mbf{w}})
+L[{\mbf{w}}]
\right):\jumptwo{{\mbf{v}}} \Bigr),
\end{equation*}
where the dependence on ${\mbf{v}}$ represents a linear mapping $\mathbb{S} \to \mathbb R$ in $\mathbb{S}^*$. Thus the derivative $\mathcal N'$ is a mapping $\mathbb{S} \to L(\mathbb{S},\mathbb{S}^*)$, where $L(\mathbb{S},\mathbb{S}^*)$ denotes the linear mappings from $\mathbb{S}$ to $\mathbb{S}^*$. Therefore the integral
\[
\mathcal{P}(t,{\mbf{v}}) := \int_0^1 \mathcal N'({\mbf{w}}^\theta(t,\cdot))({\mbf{v}})\, d\theta,
\]
where ${\mbf{w}}^\theta := \theta {\mbf{u}} + (1-\theta) {\mbf{w}}_h$,
belongs to $\mathbb{S}^*$ for each $t \in (0,T)$, ${\mbf{v}} \in \mathbb{S}$. In particular $\mathcal{P}(t, {\mbf{u}}(t,\cdot)-{\mbf{w}}_h(t,\cdot)) \in \mathbb{S}^*$ and
\begin{align}\nonumber
\mathcal{P}(t, {\mbf{u}}(t,\cdot)-{\mbf{w}}_h(t,\cdot)) \;
= & \int_0^1 \mathcal N'({\mbf{w}}^\theta(t,\cdot))({\mbf{u}}(t,\cdot)-{\mbf{w}}_h(t,\cdot)) \, d\theta\\
= & \int_0^1 \partial_\theta (\mathcal N({\mbf{w}}^\theta(t,\cdot))) \, d\theta= \, \mathcal N({\mbf{u}}(t,\cdot)) - \mathcal N({\mbf{w}}_h(t,\cdot)),\label{eq:PN}
\end{align}
using that $[0,1] \to \mathbb{S}^*, \; \theta \mapsto \mathcal N({\mbf{w}}^\theta(t,\cdot))$ is continuously differentiable as $\tilde{\mbf{p}} \in C^{1,1}(\mathbb{R}^{2n})$. We shall frequently abbreviate $\mathcal{P}(t,{\mbf{z}}(t,\cdot))$ by $\mathcal{P}{\mbf{z}}$ below.

We assume that there is an $s \in (3/2,2]$ such that for all $\mbf{\alpha} \in [L^2(\Omega)]^n$ and $\mbf{\beta} \in [H^{1/2}(\Gamma_{\inter})]^{2n}$ there exists a solution $\mbf{\zeta} \in \mathbb{H}^s$ of the linear dual equation:
\begin{equation}\label{eq:weak_dual_prob}
B({\mbf{v}},\mbf{\zeta})+\lambda \ltwoin{{\mbf{v}}}{\mbf{\zeta}}
+\langle{\mathcal{P}}{\mbf{v}},{\mbf{\zeta}}\rangle
=({\mbf{v}},\mbf{\alpha}) + \langle{\mbf{v}}, \mbf{\beta}\rangle_{\Gamma_{\inter}} \quad \forall \, {\mbf{v}}\in\mathbb{H}^1.
\end{equation}
Further, we assume that the dual solution $\mbf{\zeta}$ satisfies the elliptic regularity bound
\begin{equation}\label{eq:dual_elliptic_regularity}
  \sum_{j=1}^2\ltwo{\mbf{\zeta}}{H^s(\Omega^j)}\lesssim\ltwo{\mbf \alpha}{} + \ltwo{\mbf \beta}{H^{1/2}(\Gamma_{\inter})}.
\end{equation}
The bound is motivated by \cite{mitreamitrea} and \cite{savare}: Suppose $\Omega^1$ and $\Omega^2$ are smooth or creased domains. Assuming the existence of $\mbf{\zeta}$ allows to decouple the problem into the subdomains, using $\mbf{\zeta}$ for the boundary data on the interface. With $\mathbb{H}^1$ control available, lower-order terms may be moved to the right-hand side before applying elliptic regularity bounds such as those mentioned above, leading typically to control in fractional Sobolev norms.

\begin{lemma}\label{ell_proj_l2_lemma}
Assume that the hypothesis of Lemma \ref{ell_proj_lemma} and~\eqref{eq:weak_dual_prob} with~\eqref{eq:dual_elliptic_regularity} hold true.
For $\lambda>0$ sufficiently large, for $h_{\max}$ sufficiently small, the following error bounds holds:
\begin{align}
\label{ell_proj_l2_bound}
\ltwo{\mbf{\rho}}{}
&\le  C  (1+h^2_{\max}\lambda)^{1/2} h^{s-1}_{\max}\,
	{\ndg{\mbf{\eta}}}_{\mathcal{B},\lambda},
\end{align}
If in addition the function
$\tilde{\mbf{p}}$
defined in \eqref{eq:pjump} is twice differentiable with bounded second partials and ${\mbf{u}}, {\mbf{u}}_t \in W^{1,\infty}([0,T] \times \Gamma_{\inter})$ then
\begin{align}
\label{ell_proj_l2_bound_t}
\ltwo{\mbf{\rho}_t}{}^2
&\le  C  (1+h^2_{\max}\lambda)^{1/2} h^{s-1}_{\max}\,({\ndg{\mbf{\eta}_t}}_{\mathcal{B},\lambda}+{\ndg{\mbf{\eta}}}_{\mathcal{B},\lambda}).
\end{align}
The constant $C$ depends only on $C_A$ and the shape-regularity of the mesh.
\end{lemma}

{\bf Proof.} Let ${\mbf{z}}$ solve \eqref{eq:weak_dual_prob} with $\mbf{\alpha} = \mbf{\rho}$ and $\mbf{\beta} = \mbf{0}$.
Owing to the adjoint consistency of the symmetric interior penality bilinear form $B$ as well as of $\mathcal P$, the
dual solution ${\mbf{z}}$ also satisfies \eqref{eq:weak_dual_prob} for ${\mbf{v}} \in \mathbb{V}_h$: Because ${\mbf{z}}$ belongs to $\mathbb{H}^1$, the integral terms over $\Gamma_{\dint}$ in $B$ vanish, so that upon integration by parts $B({\mbf{v}},{\mbf{z}})$ is equal to the $L^2(\Omega)$ scalar product of ${\mbf{v}}$ and the adjoint differential operator applied to ${\mbf{z}}$, as well as $L^2$ scalar products over $\Gamma_{\inter}$ and $\partial \Omega$. In all scalar products, also including those arising from $\mathcal P$, derivatives are then only acting on ${\mbf{z}}$ and not on ${\mbf{v}}$, implying by density of smooth functions in $L^2$ that \eqref{eq:weak_dual_prob} also holds for dG test functions ${\mbf{v}}$. This allows us to test in~\eqref{eq:weak_dual_prob} with ${\mbf{v}}=\mbf{\rho}={\mbf{u}}-{\mbf{w}}_h$ and get:
\begin{equation}\label{eq:dual_l2}
\begin{aligned}
\ltwo{\mbf{\rho}}{}^2&={B}(\mbf{\rho},{\mbf{z}})+\lambda \ltwoin{\mbf{\rho}}{{\mbf{z}}}
+\langle{\mathcal{P}}\mbf{\rho},{\mbf{z}}\rangle\\
&={B}(\mbf{\rho},{\mbf{z}})+\lambda \ltwoin{\mbf{\rho}}{{\mbf{z}}}
+N({\mbf{u}},{\mbf{z}})-N({\mbf{w}}_h,{\mbf{z}})\\
&=
B(\mbf{\rho},\mbf{\eta}^z)
+\lambda \ltwoin{\mbf{\rho}}{\mbf{\eta}^z}
+N({\mbf{u}},\mbf{\eta}^z)-N({\mbf{w}}_h,\mbf{\eta}^z),
\end{aligned}
\end{equation}
with $\mbf{\eta}^z={\mbf{z}}-\Pi{\mbf{z}}$.  The second equality follows from \eqref{eq:PN}; the last equality  follows from the definition of the elliptic projection~(\ref{ellipitc_projection}).

We now bound each term on the right-hand side of~\eqref{eq:dual_l2}.
Observe that
\begin{equation}\label{eq:dual_B_bound}
\begin{aligned}
\left
B(\mbf{\rho},\mbf{\eta}^z)\right|&\le C \ndg{\mbf{\rho}}(
\ndg{\mbf{\eta}^z}_{\mathcal{B}}^2
+
\sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d} (\max_{i=1,\dots,n}\ltwo{a_i^{-1/2}B_i}{\infty,\kappa}^2)\ltwo{\mbf{\eta}^z}{\kappa}^2)^{1/2}\\
\vspace{2mm}
&\le C \ndg{\mbf{\rho}}(
\ndg{\mbf{\eta}^z}_{\mathcal{B}}^2
+\lambda\ltwo{\mbf{\eta}^z}{}^2)^{1/2}
\le C \ndg{\mbf{\rho}} \,{\ndg{\mbf{\eta}^z}}_{\mathcal{B},\lambda},
\end{aligned}
\end{equation}
for $\lambda$ big enough.
The nonlinear interface term in~(\ref{eq:dual_l2}) is bounded
as in the proof of Lemma~\ref{ell_proj_lemma}, yielding
\begin{equation}\label{eq:dual_inter_bound}
\begin{aligned}
&\left| N({\mbf{u}},\mbf{\eta}^z)\right. \left. \hspace{-1mm}-N({\mbf{w}}_h,\mbf{\eta}^z)\right|
\le 4C^p_{\mathcal{B}}\big(\sum_{j=1}^2\ltwo{\mbf{\rho}|_{\Omega^j}}{\Gamma_{\inter}}^2\big)^{1/2}
\big(\sum_{j=1}^2\ltwo{\mbf{\eta}^z|_{\Omega^j}}{\Gamma_{\inter}}^2\big)^{1/2}\\
\vspace{2mm}
&\quad\le \big(C_{coer}\ndg{\mbf{\rho}}^2+4C^p_{\mathcal{B}}c_2\epsilon^{-1}\ltwo{\mbf{\rho}}{}^2\big)^{1/2}
\big(C_{coer}\ndg{\mbf{\eta}^z}^2+4C^p_{\mathcal{B}}c_2\epsilon^{-1}\ltwo{\mbf{\eta}^z}{}^2\big)^{1/2}\\
&\quad\le \big(C_{coer}\ndg{\mbf{\rho}}^2+\lambda\ltwo{\mbf{\rho}}{}^2\big)^{1/2}
\big(C_{coer}\ndg{\mbf{\eta}^z}^2+\lambda\ltwo{\mbf{\eta}^z}{}^2\big)^{1/2},
\end{aligned}
\end{equation}
using once again (\ref{trace_estimate}) with
$\epsilon= C_{coer}\alpha_{\min}/(4C^p_\mathcal{B}c_1)$,
for $\lambda> 16(C^p_\mathcal{B})^2c_1c_2/(C_{coer}\alpha_{\min})$.

Using the above inequalities
and applying the bounds on the elliptic projection error given in
Lemma~\ref{ell_proj_lemma}, we obtain
\begin{equation}\label{rho_l2}
\ltwo{\mbf{\rho}}{}^2\le C
{\ndg{\mbf{\eta}}}_{\mathcal{B},\lambda}
{\ndg{\mbf{\eta}^z}}_{\mathcal{B},\lambda}
\end{equation}
for $\lambda$ big enough.

The term in $\mbf{\eta}^z$ can be bounded using standard approximation estimates for the error of the (orthogonal) $L^2$-projection (see, e.g., \cite{schwab}) yielding
\begin{equation}\label{eta_z}
{\ndg{\mbf{\eta}^z}}_{\mathcal{B},\lambda}^2
\le C\sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d} h_\kappa^{2(s-1)}(1+h_\kappa^2\lambda) \, \ltwo{{\mbf{z}}}{H^s(\kappa)}^2.
\end{equation}

Further, inserting~(\ref{eta_z}) into~(\ref{rho_l2}) and then using \eqref{eq:dual_elliptic_regularity}
we get
\begin{equation}\label{rho_l2_1}
\begin{aligned}
	\ltwo{\mbf{\rho}}{}^2 &\le C(1+h^2_{\max}\lambda)^{1/2} h^{s-1}_{\max}\,
	{\ndg{\mbf{\eta}}}_{\mathcal{B},\lambda}\sum_{j=1}^2\ltwo{{\mbf{z}}}{H^s(\Omega^j)}\\
	&\le C  (1+h^2_{\max}\lambda)^{1/2} h^{s-1}_{\max}\,
	{\ndg{\mbf{\eta}}}_{\mathcal{B},\lambda}
	\ltwo{\mbf{\rho}}{},\end{aligned}
\end{equation}
thus yielding~\eqref{ell_proj_l2_bound}.

We now prove~\eqref{ell_proj_l2_bound_t}. Let $\tilde{{\mbf{z}}}$ be the solution \eqref{eq:weak_dual_prob} with $\mbf{\alpha} = \mbf{\rho}_t$ and $\mbf{\beta} = \mbf{0}$.
Using the linearity of $\mathcal{P}$ in the second argument, we have
\begin{align*}
\textstyle d_t \mathcal{P}(t,\mbf{\rho}(t)) = \, & ( \partial_{\smbf{\rho}} \mathcal{P} ) \bigr|_{(t,\smbf{\rho}(t))} \mbf{\rho}_t + (\partial_t \mathcal{P})\bigr|_{(t,\smbf{\rho}(t))}
= \mathcal{P} \bigr|_{(t,\smbf{\rho}(t))} \mbf{\rho}_t + (\partial_t \mathcal{P})\bigr|_{(t,\smbf{\rho}(t))}.
\end{align*}
Testing in~\eqref{eq:weak_dual_prob} with ${\mbf{v}}=\mbf{\rho}_t$ gives, for each $t$,
\begin{equation}\label{eq:dual_l2_t}
\begin{aligned}
\ltwo{\mbf{\rho}_t(t)}{}^2&={B}(\mbf{\rho}_t,\tilde{{\mbf{z}}})+\lambda \ltwoin{\mbf{\rho}_t}{\tilde{{\mbf{z}}}}
+\langle{\mathcal{P}}\mbf{\rho}_t,{\tilde{{\mbf{z}}}}\rangle\\
&
B(\mbf{\rho}_t,\tilde{{\mbf{z}}})+\lambda \ltwoin{\mbf{\rho}_t}{\tilde{{\mbf{z}}}}
+ \langle d_t  {\mathcal{P}}(t,\mbf{\rho}(t)),{\tilde{{\mbf{z}}}}\rangle- \langle \partial_t \mathcal{P}, {\tilde{{\mbf{z}}}} \rangle\\
&
B(\mbf{\rho}_t,\tilde{{\mbf{z}}})+\lambda \ltwoin{\mbf{\rho}_t}{\tilde{{\mbf{z}}}}
+ \bigl\langle d_t \bigl( N({\mbf{w}},\cdot)-N({\mbf{w}}_h,\cdot) \bigr), \tilde{{\mbf{z}}} \bigl\rangle  - \langle \partial_t \mathcal{P}, {\tilde{{\mbf{z}}}} \rangle\\
&
B(\mbf{\rho}_t,\mbf{\eta}^{\tilde{z}})+\lambda \ltwoin{\mbf{\rho}_t}{\mbf{\eta}^{\tilde{z}}}
+ \bigl\langle d_t \bigl( N({\mbf{w}},\cdot)-N({\mbf{w}}_h,\cdot) \bigr), \mbf{\eta}^{\tilde{z}} \bigl\rangle - \langle \partial_t \mathcal{P}, {\tilde{{\mbf{z}}}} \rangle.
\end{aligned}
\end{equation}
with $\mbf{\eta}^{\tilde{z}}=\tilde{{\mbf{z}}}-\Pi\tilde{{\mbf{z}}}$.  The last equality follows from differentiating with respect to time the definition of the elliptic projection~(\ref{ellipitc_projection}).

The first three terms on the right-hand side of~\eqref{eq:dual_l2_t} can be bounded by the same argument used above. In particular, following the argument in~\eqref{eq:p_t_bound} and then applying the bounds on the elliptic projection error given in Lemma~\ref{ell_proj_lemma} yields
\begin{equation}\label{eq:Ntbound}
\begin{aligned}
& |d_t \left(N({\mbf{u}},\mbf{\eta}^{\tilde{z}})-N({\mbf{w}}_h,\mbf{\eta}^{\tilde{z}})\right)|\\
\le \; &
 4C^p_{\mathcal{B},{\mbf{u}}}\big(\sum_{j=1}^2(\ltwo{\mbf{\rho}|_{\Omega^j}}{\Gamma_{\inter}}^2+\ltwo{\mbf{\rho}_t|_{\Omega^j}}{\Gamma_{\inter}}^2)\big)^{1/2}
\big(\sum_{j=1}^2\ltwo{\mbf{\eta}^{\tilde{z}}|_{\Omega^j}}{\Gamma_{\inter}}^2\big)^{1/2}\\
\le \; & C({\ndg{\mbf{\eta}}}_{\mathcal{B},\lambda}+{\ndg{\mbf{\eta}_t}}_{\mathcal{B},\lambda}){\ndg{\mbf{\eta}^{\tilde{z}}}}_{\mathcal{B},\lambda}
\end{aligned}
\end{equation}
As for the last term, we proceed as follows.
Spelling out the definition of $\mathcal{P}$,
we have
\begin{align*}
\langle {\mathcal{P}}|_{(t,\smbf{\rho}(t))}, {\tilde{{\mbf{z}}}} \rangle
&=\Bigl\langle \int_0^1{\mathcal N}'({\mbf{w}}^{\theta})(\mbf{\rho})\,d\theta,\tilde{{\mbf{z}}}\Bigr\rangle\\
&=\int_0^1\left\langle \left(\tilde{\mbf{p}}({\mbf{w}}^\theta)+L
 \right)'[\mbf{\rho}],\jumptwo{\tilde{{\mbf{z}}}}\right\rangle d\theta\\
&=\int_0^1\left(\int_{\Gamma_{\inter}}\left(
\tilde{\mbf{p}}'({\mbf{w}}^\theta)[\mbf{\rho}])+L[\mbf{\rho}]\right)
:\jumptwo{\tilde{{\mbf{z}}}}\right) d\theta.
\end{align*}
Thus, using the assumption that the Hessian $\tilde{\mbf{p}}''$ is bounded, the embedding of $H^s(\Omega^j)$ into $L^{\infty}(\Omega^j)$, $j=1,2$,
and then using~\eqref{eq:dual_elliptic_regularity}, Lemma~\ref{broken_trace}, and~\eqref{ell_proj_bound_t}, we get:

\begin{equation}
\begin{aligned}
|\langle \partial_t{\mathcal{P}}|_{(t,\smbf{\rho}(t))}, {\tilde{{\mbf{z}}}} \rangle|
&=
\left|\int_0^1\left( \int_{\Gamma_{\inter}} \bigl(
\tilde{\mbf{p}}''({\mbf{w}}^\theta) [\partial_t {\mbf{w}}^\theta, \mbf{\rho}]\bigr) : \jumptwo{\tilde{{\mbf{z}}}}\right)d\theta \right|\\
&\le \int_0^1 \left( \int_{\Gamma_{\inter}}
|\tilde{\mbf{p}}''({\mbf{w}}^\theta)| |\partial_t {\mbf{w}}^\theta| |\mbf{\rho}|
|\jumptwo{\tilde{{\mbf{z}}}}|\right)d\theta\\
&\le C(\sum_{j=1}^2\ltwo{\tilde{{\mbf{z}}}}{\infty,\Omega^j})\int_{\Gamma_{\inter}}(|\partial_t \mbf{\rho}|+|\partial_t{\mbf{w}}_h|)|\mbf{\rho}|\\
&\le C \ltwo{\mbf{\rho}_t}{}
(\ltwo{\mbf{\rho}_t}{H^{1/2}(\Gamma_{\inter})}+\ltwo{\partial_t{\mbf{w}}_h}{H^{1/2}(\Gamma_{\inter})})\ltwo{\mbf{\rho}}{H^{1/2}(\Gamma_{\inter})^*}\\
&\le C_{{\mbf{u}}} \ltwo{\mbf{\rho}_t}{}\ltwo{\mbf{\rho}}{H^{1/2}(\Gamma_{\inter})^*}.
\end{aligned}
\label{eq:dt_rho_z}
\end{equation}

Hence we are left in need of an estimate of $ \ltwo{\mbf{\rho}}{H^{1/2}(\Gamma_{\inter})^*}$. This can be obtained by the following duality argument.

Let $\hat{\mbf{z}}$ be the solution \eqref{eq:weak_dual_prob} with $\mbf{\alpha} = 0$ and $\mbf{\beta} = \delta(\mbf{\rho})$. Here $\delta$ is the duality map \cite[IIB, p.~860]{Zeidler} between $H^{1/2}(\Gamma_{\inter})$ and $H^{1/2}(\Gamma_{\inter})^*$ such that $\ltwo{\delta(\mbf{\rho})}{H^{1/2}(\Gamma_{\inter})}=\ltwo{\mbf{\rho}}{H^{1/2}(\Gamma_{\inter})^*}$ and $\langle \mbf{\rho},\delta(\mbf{\rho})\rangle_{\Gamma_{\inter}}=\ltwo{\mbf{\rho}}{H^{1/2}(\Gamma_{\inter})^*}^2$. Testing in~\eqref{eq:weak_dual_prob} with ${\mbf{v}}=\mbf{\rho}$ yields
\begin{equation}
\begin{aligned}
\ltwo{\mbf{\rho}}{H^{1/2}(\Gamma_{\inter})^*}^2&= \langle \mbf{\rho},\delta(\mbf{\rho})\rangle_{\Gamma_{\inter}}\\
&={B}(\mbf{\rho},\hat{\mbf{z}})+\lambda \ltwoin{\mbf{\rho}}{\hat{\mbf{z}}}
+\langle{\mathcal{P}}\mbf{\rho},{\hat{\mbf{z}}}\rangle,\\
&={B}(\mbf{\rho},\hat{\mbf{z}})+\lambda \ltwoin{\mbf{\rho}}{\hat{\mbf{z}}}
+N({\mbf{u}},\hat{\mbf{z}})-N({\mbf{w}}_h,\hat{\mbf{z}})\\
&=
B(\mbf{\rho},\mbf{\eta}^{\hat z})
+\lambda \ltwoin{\mbf{\rho}}{\mbf{\eta}^{\hat z}}
+N({\mbf{u}},\mbf{\eta}^{\hat z})-N({\mbf{w}}_h,\mbf{\eta}^{\hat z}),
\label{eq:dual_l2_tt}
\end{aligned}
\end{equation}
with $\mbf{\eta}^{\hat z}=\hat{\mbf{z}}-\Pi\hat{\mbf{z}}$.

Now a bound for $\ltwo{\mbf{\rho}}{H^{1/2}(\Gamma_{\inter})^*}$ can be derived by following the same steps used above to get~\eqref{rho_l2_1}, yielding:
\begin{equation}\label{dual_l2_ttt}
\begin{aligned}
	\ltwo{\mbf{\rho}}{H^{1/2}(\Gamma_{\inter})^*}^2 &\le C(1+h^2_{\max}\lambda)^{1/2} h^{s-1}_{\max}\,
	{\ndg{\mbf{\eta}}}_{\mathcal{B},\lambda}\sum_{j=1}^2\ltwo{\hat{\mbf{z}}}{H^s(\Omega^j)}\\
	&\le C  (1+h^2_{\max}\lambda)^{1/2} h^{s-1}_{\max}\,
	{\ndg{\mbf{\eta}}}_{\mathcal{B},\lambda}
	\ltwo{\delta(\mbf{\rho})}{H^{1/2}(\Gamma_{\inter})},

	\end{aligned}
\end{equation}
having used~\eqref{eq:dual_elliptic_regularity} once again. Finally, by inserting in~\eqref{dual_l2_ttt}  the definition of $\delta(\mbf{\rho})$ we conclude that
\begin{equation}\label{dual_l2_tttt}
	\ltwo{\mbf{\rho}}{H^{1/2}(\Gamma_{\inter})^*}\le
	C  (1+h^2_{\max}\lambda)^{1/2} h^{s-1}_{\max}\,
	{\ndg{\mbf{\eta}}}_{\mathcal{B},\lambda}.
\end{equation}

Using~\eqref{dual_l2_tttt} in~\eqref{eq:dt_rho_z} we bound last term on the right-hand side of~\eqref{eq:dual_l2_t}. Now, recalling~\eqref{eq:Ntbound}, the bound~\eqref{ell_proj_l2_bound_t} easily follows from~\eqref{eq:dual_l2_t}.
\vspace{-.2cm} \\* \mbox{} \hfill $\Box$

\section{DG method for the parabolic system and its error analysis}\label{err_par_prob}
The above discussion motivates the introduction of the following IPDG-in-space method for the system (\ref{pde}), (\ref{modelbc}), (\ref{bcvec}), and (\ref{transmission}).

 For $t=0$, let ${\mbf{u}}_h(0)={\mbf{w}}_h(0)$. For $t\in(0,T]$, find ${\mbf{u}}_h\equiv {\mbf{u}}_h(t)\in \mathbb{V}_h$ such that
\begin{equation}\label{eq:semidiscrete}
\ltwoin{({\mbf{u}}_h)_t}{{\mbf{v}}_h}+B({\mbf{u}}_h,{\mbf{v}}_h)+N({\mbf{u}}_h,{\mbf{v}}_h)+ \ltwoin{\mbf{F}({\mbf{u}}_h)}{{\mbf{v}}_h}=l({\mbf{v}}_h),\texte{for all}
{\mbf{v}}_h\in \mathbb{V}_h,
\end{equation}
where
\begin{equation}\label{eq:linear}
l({\mbf{v}}_h):=\!-\!\int_{\Gamma_{\text{\rm D}}}\hspace{-1.5mm}\Big((\mbf{g}_{\text{\rm D}}\otimes\mbf{n}):(A\nabla{\mbf{v}}_h)+(\mathcal{X}^-G_{\text{\rm D}}\mbf{B}):({\mbf{v}}_h\otimes\mbf{n})-\Sigma
\mbf{g}_{\text{\rm D}}\cdot{\mbf{v}}_h\Big)
+\int_{\Gamma_{\text{\rm N}}}\!\mbf{g}_{\text{\rm N}}\cdot{\mbf{v}}_h,
\end{equation}
noting that $(\mbf{q}\otimes\mbf{n}):(\mbf{v}\otimes\mbf{n})= \mbf{q}\cdot \mbf{v}$, for $\mbf{q},\mbf{v}\in\mathbb{R}^n$, and having denoted
 $G_{\text{\rm D}}:=\diag(g_{\text{\rm D}}^1,\dots,g_{\text{\rm D}}^n)$.

We shall make use of the following result to treat the nonlinear reaction term.

\begin{lemma}\label{pf_lemma}
If $d=2$, let $q\in[1,\infty)$ and if $d=3$, let $q\in[1,6]$. Then, there exists a constant $C_{\rm PF}>0$, depending only on the geometry of the subdomains $\Omega^j$, $j \in \{ 1,2 \}$ and of the Dirichlet boundary, such that,
for all $v\in \mathbb{S}^1$, we have
\begin{equation}\label{pf}
\ltwo{v}{q}^2\le C_{\rm PF}\max\{1,\alpha_{\min}^{-1}\}\ndg{v}^2.
\end{equation}
 Moreover, for $0\le\gamma< 2$, if $d=2$, and for $0\le\gamma\le 4/3$, if $d=3$, we have
\begin{equation}\label{pf_useful}
\ltwo{v}{\gamma+2}^{\gamma+2}\le C_{\rm PF}\max\{1,\alpha_{\min}^{-1} \}\ltwo{v}{}^{\gamma}\ndg{v}^2.
\end{equation}
Finally, for $d=2$ and $\gamma=2$ (which corresponds to $q=\infty$), (\ref{pf}) and (\ref{pf_useful}) hold for $v\in\mathbb{V}_h$, with $C_{\rm PF}:=C |\log(\min_{\Omega}\mbf{\tt h})|$ for some constant $C>0$ depending only on $\Omega^i$.
\end{lemma}

{\bf Proof.}
From the assumptions on the topology of the subdomains, we can apply Theorem 3.7 from \cite{lasis-suli:03}, to deduc
\[
\ltwo{v}{q,\Omega^j}^2\le C\max\{1,\alpha_{\min}^{-1}\}\ndg{v}^2,
\]
with the boundary contribution taken as $\partial\Omega^j\backslash \Gamma_{\inter}$, for $1\le q<\infty$ if $d=2$ and for $1\le q \le 6 $ if $d=3$. For the case $q=\infty$ if $d=2$, we make use of the standard inverse estimate
$\ltwo{v}{\infty,\kappa}^2\le C\log|h_{\kappa}|  \ltwo{v}{H^1(\kappa)}^2$.
Setting $p=2/\gamma$, for $1/p+1/q=1$, which gives $q =2/(2-\gamma)$, H\"older's inequality implies (\ref{pf_useful}).\vspace{-.2cm} \\* \mbox{} \hfill $\Box$

We are now ready to prove a bound for the difference between the elliptic projection ${\mbf{w}}_h$ and the dG approximation ${\mbf{u}}_h$.

\begin{theorem}\label{sd_theorem_one}
Consider the notation of Lemma \ref{ell_proj_lemma} and let $\gamma$ as in Lemma \ref{pf_lemma}. Assume that ${\mbf{u}} \in L^2(0,T;\mathbb{H}^s)\cap L^\infty(0,T\times \Omega)$, $s>3/2$, ${\mbf{u}}_t \in L^2(0,T;L^2(\Omega))$, and that $h_{\max}$ is small enough. Then, we have
\[
         \ltwo{\mbf{\mbf{\theta}}}{L^{\infty}(0,T;L^2(\Omega))}^2+C_{coer}
         \ltwo{\mbf{\mbf{\theta}}}{L^2(0,T;\mathbb{S})}^2\le 4\delta^2{\rm e}^{\tilde{C}T},
\]
where  $\mbf{\theta}:={\mbf{w}}_h-{\mbf{u}}_h$, $\tilde{C}:=C(\ltwo{{\mbf{u}}}{\infty,[0,T]\times\Omega},\lambda)$, and
\begin{equation}\label{delta}
\delta^2(t) = \int_0^T  \tilde{C}\ltwo{\mbf{\rho} }{}^2+ C\Big(\ltwo{\mbf{\rho} }{\gamma+2}^{\gamma+2}+
\lambda^{-1}\ltwo{\mbf{\rho} _t}{}^2\Big) \ud t .
\end{equation}
\end{theorem}

{\bf Proof.} Without loss of generality assume that $\gamma>0$; the case $\gamma=0$ follows by a simple modification of the argument presented below.

Let $\mbf{e}:=\mbf{\rho}+\mbf{\theta}$, with $\mbf{\rho}:={\mbf{u}}-{\mbf{w}}_h$ and $\mbf{\theta}:={\mbf{w}}_h-{\mbf{u}}_h$.
We note that the continuity of ${\mbf{w}}_h$ in the time variable is implied by the well-posedness of the elliptic projection problem, and that of ${\mbf{u}}_h$ by a standard local existence argument near $0$ on the resulting system of ordinary differential equations. Hence, $\mbf{\theta}$ is continuous in $[0,T]$ with $\mbf{\theta}(0)=\mathbf{0}$.

Orthogonality implies
\[
\ltwoin{\partial_t \mbf{e}}{\mbf{\theta}}+B(\mbf{e},\mbf{\theta})
+N({\mbf{u}},\mbf{\theta})-N({\mbf{u}}_h,\mbf{\theta})+\ltwoin{\mbf{F}({\mbf{u}})-\mbf{F}({\mbf{u}}_h)}{\mbf{\theta}}
= 0.
\]
From (\ref{ellipitc_projection}), we deduce
\begin{equation}\label{energy_identity}
\frac{1}{2} \frac{\ud}{\ud t}\ltwo{\mbf{\theta}}{}^2+B(\mbf{\theta},\mbf{\theta})=
 \ltwoin{\mbf{F}({\mbf{u}}_h)-\mbf{F}({\mbf{u}})}{\mbf{\theta}}+N({\mbf{u}}_h,\mbf{\theta})-N({\mbf{w}}_h,\mbf{\theta})
+\ltwoin{\lambda \mbf{\rho}-\mbf{\rho} _t}{\mbf{\theta}}.
\end{equation}

Making use of the inequality
\[
\int_{\Omega} a^{\gamma+1}b\le
\frac{\gamma+1}{\gamma+2}\ltwo{a}{\gamma+2}^{\gamma+2}+\frac{1}{\gamma+2}\ltwo{b}{\gamma+2}^{\gamma+2},
\]  for $a,b>0$, the nonlinear reaction term can be bounded as follows:
\begin{equation}
\begin{aligned}
|\ltwoin{\mbf{F}({\mbf{u}})-\mbf{F}({\mbf{u}}_h)}{\mbf{\theta}}|
&\le C \int_{\Omega} (1 + |{\mbf{u}}|^{\gamma}+|{\mbf{u}}_h|^{\gamma})|{\mbf{u}}-{\mbf{u}}_h||\mbf{\theta}|\\
&\le C \int_{\Omega} (1 + |{\mbf{u}}|^{\gamma})|{\mbf{u}}-{\mbf{u}}_h||\mbf{\theta}|+ C \int_{\Omega} |{\mbf{u}}-{\mbf{u}}_h|^{\gamma+1}|\mbf{\theta}|\\
&\le C({\mbf{u}}) \big(\ltwo{\mbf{\rho} }{}^2+\ltwo{\mbf{\theta}}{}^2\big)
+ C\big(\ltwo{\mbf{\rho} }{\gamma+2}^{\gamma+2}+\ltwo{\mbf{\theta}}{\gamma+2}^{\gamma+2}\big).
\end{aligned}
\end{equation}

Also, using the regularity of $\tilde{\mbf{p}}$ and
(\ref{trace_estimate}), we have
\begin{equation}\label{nterms_bound}
\begin{aligned}
|N({\mbf{u}}_h,\mbf{\theta})-N({\mbf{w}}_h,\mbf{\theta})|
&\le C^p_{\mathcal{B}}\sum_{j=1}^2\ltwo{\mbf{\theta}|_{\Omega^j}}{\Gamma_{\inter}}^2
\le \frac{1}{4}C_{\rm coer}\ndg{\mbf{\theta}}^2+\frac{\lambda}{2}\ltwo{\mbf{\theta}}{}^2,
\end{aligned}
\end{equation}
choosing $\epsilon$ and $\lambda$ as in the proof of Lemma \ref{ell_proj_lemma}.

We also have
\begin{equation}\label{rho_and_rhot_bound}
|\ltwoin{\lambda\mbf{\rho} -\mbf{\rho} _t}{\mbf{\theta}}| \le \frac{\lambda}{2} \ltwo{\mbf{\rho}}{}^2+\frac{1}{2\lambda}\ltwo{\mbf{\rho} _t}{}^2+\lambda\ltwo{\mbf{\theta}}{}^2.
\end{equation}
Lemma \ref{pf_lemma} implies
\begin{equation}\label{imbedding_sd}
\ltwo{\mbf{\theta}}{\gamma+2}^{\gamma+2}
\le C_{\rm PF}\max\{1,\alpha_{\min}^{-1}\} \ltwo{\mbf{\theta}}{}^{\gamma} \ndg{\mbf{\theta}}^2.
\end{equation}
Using these bounds in (\ref{energy_identity}), along with the coercivity and continuity bounds from Lemmas \ref{coercivity} and \ref{continuity}, integrating the resulting inequality with respect to $t$ (and multiplying by $2$) between $0$ and $\tau$, we arrive at
\begin{equation}\label{energy_bound_indermediate_two}
\begin{aligned}
& \ltwo{\mbf{\theta}(\tau)}{}^2+C_{coer}\int_0^{\tau}\ndg{\mbf{\theta}}^2\\ \le&
\delta^2+\tilde{C}\int_0^{\tau}\ltwo{\mbf{\theta}}{}^2
+\hat{C} \esssup_{t\in[0,\tau]}\ltwo{\mbf{\theta}(t)}{}^{\gamma}\int_0^{\tau} \ndg{\mbf{\theta}}^2\\
\le& \delta^2+\tilde{C}\!\int_0^{\tau}\ltwo{\mbf{\theta}}{}^2
+\hat{C} \Big(\esssup_{t\in[0,\tau]}\ltwo{\mbf{\theta}(t)}{}^2+C_{coer}\!\int_0^{\tau} \ndg{\mbf{\theta}}^2\Big)^{1+\gamma/2},
\end{aligned}
\end{equation}
for $\delta^2$ and $\tilde{C}$ as in the statement of the theorem,
and $\hat{C}:=C\max\{1,(C_{coer}\alpha_{\min})^{-1}\}$.
From the assumption that the mesh-size $h_{\max}$ is small enough, we can have
\[
{\delta}\le (4\hat{C}{\rm e}^{\tilde{C}T})^{-\frac{2+\gamma}{2\gamma}},
\]
which implies $(4\hat{C}\delta^2{\rm e}^{\tilde{C}T})^{1+\gamma/2}\le \delta^2$. We now consider the set
\[
I:=\big\{\tau\in[0,T]: \esssup_{t\in[0,\tau]}\ltwo{\mbf{\theta}(t)}{}^2+C_{coer}\int_0^{\tau} \ndg{\mbf{\theta}}^2\le 4\delta^2 {\rm e}^{\tilde{C}T}\big\},
\]
which is non-empty and closed due to the continuity of $\mbf{\theta}$ with respect to the time variable, as $\mbf{\theta}(0)=0$. We set $\tau^*=\max I$ and we suppose that
$\tau^*<T$.
Hence, for $\tau\le \tau^*$
, we have
\[
         \esssup_{t\in[0,\tau]}\ltwo{\mbf{\theta}(t)}{}^2+C_{coer}\int_0^{\tau} \ndg{\mbf{\theta}}^2\le 2\delta^2
   +\tilde{C}\int_0^{\tau}\ltwo{\mbf{\theta}}{}^2.
\]
Gronwall's Lemma then implies
\begin{equation}\label{after_gronwall}
         \ltwo{\mbf{\theta}(\tau^*)}{}^2+C_{coer}\int_0^{\tau^*} \ndg{\mbf{\theta}}^2\le 2\delta^2{\rm e}^{\tilde{C}T},
\end{equation}
setting $\tau=\tau^*$, which contradicts the hypothesis $\tau^*<T$, due to the continuity of the left-hand side of (\ref{after_gronwall}). Hence, $I=[0,T]$
and the result holds.\vspace{-.2cm} \\* \mbox{} \hfill $\Box$

\begin{corollary}\label{sd_cor}
Let $0\le\gamma\le 2$ if $d=2$, and $0\le\gamma\le 4/3$ if $d=3$. With the assumptions of Lemma \ref{ell_proj_lemma} and of Theorem \ref{sd_theorem_one} and ${\mbf{u}}|_{\kappa}\in H^1(0,T;[H^{k_{\kappa}+1}(\kappa)]^n)$, $k_{\kappa}\ge 1$, $ \kappa \in \mathcal{T}$, we have
\begin{equation}\label{apriori_energy}
   \ltwo{\mbf{{\mbf{u}}-{\mbf{u}}_h}}{L^{\infty}(0,T;L^2(\Omega))}^2+C_{coer}
         \ltwo{{\mbf{u}}-{\mbf{u}}_h}{L^2(0,T;\mathbb{S})}^2 \le C\mathcal{E}((0,T],\mbf{h},{\mbf{u}},\mathbb{V}_h),
\end{equation}
for $C$ independent of $\mbf{h}$ and
\begin{equation}\label{curly_E}
\mathcal{E}((0,T],\mbf{h},{\mbf{u}},\mathbb{V}_h):=\sum_{\kappa\in\mathcal{T}}}\newcommand{\ud}{\mathrm{d} \int_0^T h_{\kappa}^{2s_{\kappa}}\big(
|{\mbf{u}}|_{[H^{k_{\kappa}+1}(\kappa)]^n}^2
+|{\mbf{u}}_t|_{[H^{k_{\kappa}+1}(\kappa)]^n}^2\big),
\end{equation}
for $s_{\kappa}=\min\{m_{\kappa},k_{\kappa}\}$. Moreover, with the additional assumptions of Lemma \ref{ell_proj_l2_lemma}, we have
\begin{equation}\label{apriori_l2}
   \ltwo{\mbf{{\mbf{u}}-{\mbf{u}}_h}}{L^{\infty}(0,T;L^2(\Omega))}^2 \le C h_{\max}^{s-1}\mathcal{E}((0,T],\mbf{h},{\mbf{u}},\mathbb{V}_h),
\end{equation}
for $s\in(1,2]$ dictated by the regularity of the exact solution of the dual problem \eqref{eq:weak_dual_prob}.
\end{corollary}
\begin{proof}
We begin by using (\ref{pf}) to deduce the estimates
$
\ltwo{\mbf{\rho}}{\gamma+2}^{\gamma+2}\le C \ndg{\mbf{\rho}}^{\gamma+2}
$, $
\ltwo{\mbf{\rho}}{}^{}\le C \ndg{\mbf{\rho}}^{}
$, and $
\ltwo{\mbf{\rho}_t}{}^{}\le C \ndg{\mbf{\rho}_t}^{}
$.
 An application of Lemma \ref{ell_proj_lemma} on the resulting
norms, along with standard approximation estimates for the error of the (orthogonal) $L^2$-projection (see, e.g., \cite{houston-etal:01}) already implies \eqref{apriori_energy}.

For \eqref{apriori_l2}, we set $ \varepsilon = \min\{2,k\gamma/(s-1)\}$. Then, we have the bound
\[
\ltwo{\mbf{\rho}}{\gamma+2}^{\gamma+2}\le \ltwo{\mbf{\rho}}{}^{2-\varepsilon}\ltwo{\mbf{\rho}}{2(\gamma/\varepsilon+1)}^{\gamma+\varepsilon}\le C \ltwo{\mbf{\rho}}{}^{2-\varepsilon}\ndg{\mbf{\rho}}^{\gamma+\varepsilon},
\]
which implies (for this choice of $\varepsilon$, since then $\gamma/\varepsilon+1\le 2$) that
\[
\ltwo{\mbf{\rho}}{\gamma+2}^{\gamma+2}\le C h^{s-1}_{\max}\mathcal{E}((0,T],\mbf{h},{\mbf{u}},\mathbb{V}_h).
\]
Lemma \ref{ell_proj_l2_lemma} and standard approximation estimates now imply \eqref{apriori_l2}.
\end{proof}

As it can be seen from the above proof, the use of the novel elliptic projection, introduced and analyzed in Section \ref{err_ell_problem}, plays crucial role in the treatment of the interface nonlinearity $\mbf{P}(\cdot)$, while the continuation argument, along with the Gagliardo-Nirenberg-type inequality \eqref{imbedding_sd}, treats the nonlinear reaction $\mbf{F}(\cdot)$.

\begin{remark}
We note that the bounds in Theorem \ref{sd_theorem_one} and, correspondingly, in Corollary \ref{sd_cor} do not require any global  mesh quasi-uniformity assumptions.  This result also holds on domains without internal interfaces as in the setting of \cite{Lasis07}. There, a different continuation argument is used for deriving a priori bounds which requires the mesh to be globally quasiuniform (for non-symmetric spatial operators), albeit for a larger range of $\gamma$ than the one considered in the present work.
\end{remark}

\begin{remark}
It is interesting to note that, due to the careful use of $L^2$-projection operators in conjunction with assumption \eqref{wind_assumption}, the constant in \eqref{apriori_energy} depends  on the Pecl\'et number only through the control of the nonlinear interface conditions. Thus, the present bound produces  error control that is Pecl\'et number independent for the dG method applied to a single domain  problem.
\end{remark}

We remark that it is possible to show optimal error estimates with less restrictive assumptions than \eqref{growth} on the growth of the reaction term. Indeed, assuming only that $F$ is locally Lipschitz, an optimal a priori bound can be proven, subject to certain conditions on the mesh.  This argument, motivated by ideas presented in~\cite{MR1701828,MR2299767}, for different problems will be considered elsewhere.

\section{Numerical examples}\label{sec:numex}

In view of the numerical tests, we introduce a fully discrete discretisation
for the system (\ref{pde}), (\ref{modelbc}), (\ref{transmission}).
The time discretization is based on a second-order linearly implicit method
analysed in~\cite{MR1701828,Akrivis04}.

Let us write the semidiscrete DG formulation~\eqref{eq:semidiscrete} in matrix notation as
\begin{equation}\label{eq:algebric}
{\rm M \rm U_t = LU + F(U)}
\end{equation}
where ${\rm L}$ collects all linear terms and ${\rm F}$ the nonlinear terms in \eqref{eq:semidiscrete}.
Treating the linear terms with the second-order Adams-Moulton method (trapezium rule) and the
nonlinear terms with the second-order Adams-Bashforth method  yields the following
fully discrete method (AB2-AM2):
\[
{\rm M} U^{n+1}= {\rm M}U^n + k (\theta {\rm L} U^{n+1}+(1-\theta) {\rm L} U^{n})+\frac{k}{2}(3 {\rm F}(U^n)-{\rm F}(U^{n-1})).
\]

{\bf Convergence test.} We test the validity of the IPDG method and above error bounds on a system of two equations for which the exact solution is known.

Let the domain $\Omega=[-1,1]^2$ be subdivided into two subdomains interfacing at $x=0$; thus
$\Omega^1=[-1,0]\times [-1,1]$ and $\Omega^2=[0,1]\times [-1,1]$.
We set $\Gamma_{\text{\rm D}}=\{ \pm1\}\times [-1,1]$ and
$\Gamma_{\text{\rm N}}=(-1,1)\times\{ \pm1\}$ and
impose homogeneous Dirichlet and Neumann boundary conditions
on $\Gamma_{\text{\rm D}}$ and $\Gamma_{\text{\rm N}}$, respectively, as shown in Figure~\ref{fig:domain} (left).
\begin{figure}[!t]
\begin{center}
\setlength{\unitlength}{2cm}
\begin{picture}(0,0)
\put(0,0){\line(0,1){2}}
\put(1,0){\line(0,1){2}}
\put(0,0){\line(1,0){1}}
\put(0,1){\line(1,0){1}}
\put(0,2){\line(1,0){1}}
\put(.5,.5){\makebox(0,0){\tiny{$\Omega^1$}}}
\put(.5,1.5){\makebox(0,0){\tiny{$\Omega^2$}}}
\put(.7,1.06){\makebox(0,0){\tiny{{{$\Gamma_{\inter}$}}}}}
\put(.43,1.2){\makebox(0,0){\tiny{$x$}}}
\put(.3,1.1){\makebox(0,0){\tiny{$y$}}}
\put(-.12,1){\makebox(0,0){\rotatebox{90}{\tiny{$\frac{\partial {\mbf{u}}}{\partial\mbf{n}}=\mbf{0}$}}}}
\put(1.1,1){\makebox(0,0){\rotatebox{90}{\tiny{$\frac{\partial {\mbf{u}}}{\partial\mbf{n}}=\mbf{0}$}}}}
\put(.5,1){\vector(-1,0){.2}}
\put(.5,1){\vector(0,1){.2}}
\put(.5,-.1){\makebox(0,0){\tiny{${\mbf{u}}=\mbf{0}$}}}
\put(.5,2.1){\makebox(0,0){\tiny{${\mbf{u}}=\mbf{0}$}}}
\end{picture}
\hspace{25 mm}
\begin{overpic}[scale=0.4]{u1_541.pdf}
 \small{\put(80,3){$x$}
        \put(26,3){$y$}
        \put(4,37){$u_1$}
       }
\end{overpic}
\end{center}
\caption{Convergence test. Solution domain and boundary conditions (left). First component of the exact solution at
the final time $t=1$ (right).}
\label{fig:domain}
\end{figure}
For $t>0$ we consider the system of two advection-diffusion equations
\begin{equation}\label{ex2:system}
\left\{
\begin{array}{l}
\vspace{2 mm}
u_t - \Delta u -u_x   = f^u +
\left\{
\hspace{-2 mm}
\begin{array}{ll}
\vspace{1 mm}
  u^2 - v (1-v) &\textrm{in $\Omega^1$}\\
 - v &\textrm{in } \Omega^2,\\
\end{array}
\right.\\
v_t - \Delta v -v_x = f^v+u \quad\textrm{in $\Omega^1\cup \Omega^2$},
\end{array}
\right.
\end{equation}
and accordingly construct the forcing terms
$f^{u,v}:[0,1]\times\Omega\rightarrow\mathbb{R}$ in order to
yield as exact solution
\begin{equation}\label{ex2:solution}
\left(\hspace{-1 mm}
\begin{array}{c}
\vspace{1 mm}
u\\
v
\end{array}
\hspace{-1 mm}
\right)      =
\left(\hspace{-1 mm}
\begin{array}{c}
\vspace{1 mm}
\cos t\\
\sin t
\end{array}
\hspace{-1 mm}
\right)
e^{(y^2-1)^2}
\times
\left\{
\begin{array}{ll}
\vspace{1 mm}
4 x (1 + x)
 & \textrm{in $\Omega^1$,}\\
(- 4 x^3  + 3 x + 1) &\textrm{in $\Omega^2$}.
\end{array}
\right.
\end{equation}
The first component of the solution at
time $t=1$ is shown in
Figure~\ref{fig:domain} (right).
The functions in~\eqref{ex2:solution} are compatible
with the interface conditions~\eqref{transmission}
with respect to the interface parameters
\[
\mbf{P}=\textrm{diag}(3,3),
\qquad
\Upsilon^{1} =\textrm{diag}(1,1),
\qquad
\Upsilon^{2} =\textrm{diag}(0,0),
\qquad
{\rm R} =\textrm{diag}(1,1),
\]
which are therefore used to close
problem~\eqref{ex2:system} with
conditions~\eqref{transmission}.
We tested the convergence rate of our method
under space discretisation refinement.
A fixed time step of size $.5\times 10^{-3}$ is used throughout while an initial uniform square $4\times 4$ mesh is uniformly refined. The value $C_{\sigma}=10$ was used throughout for the discontinuity-penalization parameter costant.
Numerical results are reported in Table~\ref{table:ex2} in the cases of bilinear and biquadratic spatial discretizations.
The predicted error rates of convergenece are confirmed in  both the $L^2(0,1;\mathbb{S})$, (viz., $\|\cdot\|^2_{L^2(0,1;\mathbb{S})} = \int_0^1\ndg{\cdot}^2$,) and ${L^{\infty}(0,1;L^2(\Omega))}$ norms.
\begin{table}[!t]
    \label{table:ex2}
\begin{center}
\begin{tabular}{|r|r|c|c|c|c|} \hline
\# cells & \# dofs &
\multicolumn{2}{|c|}{${L^2(0,1;\mathbb{S})}$} &
\multicolumn{2}{|c|}{${L^{\infty}(0,1;L^2(\Omega))}$}\\ \hline
\multicolumn{6}{c}{${\tt m}\equiv 1$}\\ \hline
4 & 32 & 1.579e+01 & - & 2.498e+00 & -\\ \hline
16 & 128 & 7.617e+00 & 1.05 & 8.293e-01 & 1.59\\ \hline
64 & 512 & 3.615e+00 & 1.08 & 2.400e-01 & 1.79\\ \hline
256 & 2048 & 1.736e+00 & 1.06 & 6.497e-02 & 1.89\\ \hline
1024 & 8192 & 8.475e-01 & 1.03 & 1.693e-02 & 1.94\\ \hline
4096 & 32768 & 4.182e-01 & 1.02 & 4.324e-03 & 1.97\\ \hline
\multicolumn{6}{c}{${\tt m}\equiv 2$}\\ \hline
4 & 72 & 4.781e+00 & - & 3.812e-01 & -\\ \hline
16 & 288 & 1.013e+00 & 2.24 & 5.422e-02 & 2.81\\ \hline
64 & 1152 & 2.282e-01 & 2.15 & 7.330e-03 & 2.89\\ \hline
256 & 4608 & 5.480e-02 & 2.06 & 1.240e-03 & 2.56\\ \hline
1024 & 18432 & 1.349e-02 & 2.02 & 1.743e-04 & 2.83\\ \hline
4096 & 73728 & 3.427e-03 & 1.98 & 2.259e-05 & 2.95\\ \hline
\end{tabular}
\end{center}
\caption{Convergence test. Errors and convergence rates under uniform mesh refinement.
A fixed time step of size $.5\times 10^{-3}$
was used.
Bilinear (above) and biquadratic (below) dG in space discretizations.
}\vspace{-1mm}
\end{table}

{\bf Advection-dominated test.} With this numerical example we test the robustness of the method in the advection-dominated regime.
We consider once again the domain $\Omega=[-1,1]^2$ subdivided into the  subdomains $\Omega^1$ and $\Omega^2$ already considered in the previous test. In $\Omega$, we solve a scalar equation with no reaction and the
constant diffusion $a_1=10^{-2}$
and the advection field $B_1=(0.5,0.5)$.
On the boundary $\partial\Omega$, we set homogeneous Neumann boundary conditions.
 The parameters
\[
\mbf{P}=2/10,
\qquad
\Upsilon^{1} =5/6,
\qquad
\Upsilon^{2} =1/6,
\qquad
{\rm R} =6/10,
\]
are used in the  interface conditions~\eqref{transmission}. In this case, transfer across the interface is mainly advection-driven. Further, setting the friction coefficient to
less than one models the case in which the interface acts as a filtering wall on the advected quantity; hence a boundary layer in the upwind subdomain in the proximity of the interface is expected.

We solve the problem on a uniform $16\times 16$ mesh using bilinear elements. Such mesh is not fine enough to resolve the layer forming in the proximity of the interface where the solution is also discontinuous, see Figure~\ref{fig:ex2sol:16}. Nevertheless, the numerical solution is stable, and the expected behaviour of the solution is accurately captured, as we can see by comparison with the solution obtained with a layer resolving $64\times 64$ mesh and shown in Figure~\ref{fig:ex2sol:64}.
\begin{figure}
\centering
\includegraphics[scale=0.12]{ref4_sol_0_10_20_30_40_50_60_70.pdf}
\caption{Advection-dominated test. Snapshots of the solution computed on a uniform $16\times 16$ mesh using bilinear elements: the initial condition (top-left) followed by the solution at time intervals of $0.5$.}
\label{fig:ex2sol:16}
\end{figure}
\begin{figure}
\centering
\includegraphics[scale=0.16]{ref6_sol_10_20.pdf}
\caption{Advection-dominated test. Snapshots of the solution computed on a uniform $64\times 64$ mesh using bilinear elements: solution at $t=0.5$ (left) and $t=1$ (right) corresponding to the second and third plots in the first row of Figure~\ref{fig:ex2sol:16}, respectively.}
\label{fig:ex2sol:64}
\end{figure}
We note that the method remains stable when smaller values that $10^{-2}$  are used for the diffusion $a_1$; these results are omitted for brevity.

\section{Concluding remarks}\label{conclusions}
A dG method for the numerical solution of nonlinear interface problems modelling mass transfer through semi-permeable membranes is presented and a priori error bounds are shown under typical regularity  assumptions. The good performance of the method is highlighted  through numerical experiments. A number of extensions of the presented results can be made with modest modifications.  For instance, $hp$-version error bounds can be shown and more general convection coefficients $\mbf{B}$ can be treated. We refrained from doing so in the interest of simplicity of the presentation.
 Interesting directions of further research are the consideration of the variational crimes due to inexact representation of the interface manifold, using, e.g., unfitted finite elements~\cite{BarEll:87,MR1622502} or the related approach in~\cite{Melenk_interface10}, and the treatment of more general interface nonlinearities. These will be considered elsewhere.

\bibliographystyle{siam}

\section{Numerical examples}\label{sec:numex}

In view of the numerical tests, we introduce a fully discrete discretisation
for the system (\ref{pde}), (\ref{modelbc}), (\ref{transmission}).
The time discretization is based on a second-order linearly implicit method
analysed in~\cite{MR1701828,Crouzeix,Akrivis04}.

Let us write the semidiscrete DG formulation~\eqref{eq:semidiscrete} in matrix notation as
\begin{equation}\label{eq:algebric}
{\rm M \rm U_t = LU + F(U)}
\end{equation}
where ${\rm L}$ collects all linear terms and ${\rm F}$ the nonlinear terms in~\eqref{eq:semidiscrete}.
Treating the linear terms with the second-order Adams–Moulton method (trapezium rule) and the
nonlinear terms with the second-order Adams–Bashforth method  yields the following
fully discrete method (AB2-AM2):
\[
M U^{n+1}= MU^n + k (\theta A U^{n+1}+(1-\theta) A U^{n})+\frac{k}{2}(3 F(U^n)-F(U^{n-1})).
\]

{\bf Convergence test.} We test the validity of the IPDG method and above error bounds on a system of two equations for which the exact solution is known.

Let the domain $\Omega=[-1,1]^2$ be subdivided into two subdomains interfacing at $x=0$; thus
$\Omega^1=[-1,0]\times [-1,1]$ and $\Omega^2=[0,1]\times [-1,1]$.
We set $\Gamma_{\text{\rm D}}=\{ \pm1\}\times [-1,1]$ and
$\Gamma_{\text{\rm N}}=(-1,1)\times\{ \pm1\}$ and
impose homogeneous Dirichlet and Neumann boundary conditions
on $\Gamma_{\text{\rm D}}$ and $\Gamma_{\text{\rm N}}$, respectively, as shown in Figure~\ref{fig:domain} (left).
\begin{figure}[!t]
\begin{center}
\setlength{\unitlength}{2cm}
\begin{picture}(0,0)
\put(0,0){\line(0,1){2}}
\put(1,0){\line(0,1){2}}
\put(0,0){\line(1,0){1}}
\put(0,1){\line(1,0){1}}
\put(0,2){\line(1,0){1}}
\put(.5,.5){\makebox(0,0){\tiny{$\Omega^1$}}}
\put(.5,1.5){\makebox(0,0){\tiny{$\Omega^2$}}}
\put(.7,1.06){\makebox(0,0){\tiny{{{$\Gamma_{\inter}$}}}}}
\put(.43,1.2){\makebox(0,0){\tiny{$x$}}}
\put(.3,1.1){\makebox(0,0){\tiny{$y$}}}
\put(-.12,1){\makebox(0,0){\rotatebox{90}{\tiny{$\frac{\partial {\mbf{u}}}{\partial\mbf{n}}=\mbf{0}$}}}}
\put(1.1,1){\makebox(0,0){\rotatebox{90}{\tiny{$\frac{\partial {\mbf{u}}}{\partial\mbf{n}}=\mbf{0}$}}}}
\put(.5,1){\vector(-1,0){.2}}
\put(.5,1){\vector(0,1){.2}}
\put(.5,-.1){\makebox(0,0){\tiny{${\mbf{u}}=\mbf{0}$}}}
\put(.5,2.1){\makebox(0,0){\tiny{${\mbf{u}}=\mbf{0}$}}}
\end{picture}
\hspace{25 mm}
\begin{overpic}[scale=0.4]{u1_541.pdf}
 \small{\put(80,3){$x$}
        \put(26,3){$y$}
        \put(4,37){$u_1$}
       }
\end{overpic}
\end{center}
\caption{Convergence test. Solution domain and boundary conditions (left). First component of the exact solution at
the final time $t=1$ (right).}
\label{fig:domain}
\end{figure}
For $t>0$ we consider the system of two advection-diffusion equations
\begin{equation}\label{ex2:system}
\left\{
\begin{array}{l}
\vspace{2 mm}
u_t - \Delta u -u_x   = f^u +
\left\{
\hspace{-2 mm}
\begin{array}{ll}
\vspace{1 mm}
  u^2 - v (1-v) &\textrm{in $\Omega^1$}\\
 - v &\textrm{in } \Omega^2,\\
\end{array}
\right.\\
v_t - \Delta v -v_x = f^v+u \quad\textrm{in $\Omega^1\cup \Omega^2$},
\end{array}
\right.
\end{equation}
and accordingly construct the forcing terms
$f^{u,v}:[0,1]\times\Omega\rightarrow\mathbb{R}$ in order to
yield as exact solution
\begin{equation}\label{ex2:solution}
\left(\hspace{-1 mm}
\begin{array}{c}
\vspace{1 mm}
u\\
v
\end{array}
\hspace{-1 mm}
\right)      =
\left(\hspace{-1 mm}
\begin{array}{c}
\vspace{1 mm}
\cos t\\
\sin t
\end{array}
\hspace{-1 mm}
\right)
e^{(y^2-1)^2}
\times
\left\{
\begin{array}{ll}
\vspace{1 mm}
4 x (1 + x)
 & \textrm{in $\Omega^1$,}\\
(- 4 x^3  + 3 x + 1) &\textrm{in $\Omega^2$}.
\end{array}
\right.
\end{equation}
The first component of the solution at
time $t=1$ is shown in
Figure~\ref{fig:domain} (right).
The functions in~\eqref{ex2:solution} are compatible
with the interface conditions~\eqref{transmission}
with respect to the interface parameters
\[
\mbf{P}=\textrm{diag}(3,3),
\qquad
\Upsilon^{1} =\textrm{diag}(1,1),
\qquad
\Upsilon^{2} =\textrm{diag}(0,0),
\qquad
{\rm R} =\textrm{diag}(1,1),
\]
which are therefore used to close
problem~\eqref{ex2:system} with
conditions~\eqref{transmission}.

We tested the convergence rate of our method
under space discretisation refinement.
A fixed time step of size $.5\times 10^{-3}$ is used throughout while an initial uniform square $4\times 4$ mesh is uniformly refined. The value $C_{\sigma}=10$ was used throughout for the discontinuity-penalization parameter costant.
Numerical results are reported in Table~\ref{table:ex2} in the cases of linear and quadratic spatial discretizations.
The predicted error rates of convergenece are confirmed in  both the ${L^2(0,1;\mathbb{S})}$ and ${L^{\infty}(0,1;L^2(\Omega))}$ norms.
\begin{table}[!t]
    \label{table:ex2}
\begin{center}
\begin{tabular}{|r|r|c|c|c|c|} \hline
\# cells & \# dofs &
\multicolumn{2}{|c|}{${L^2(0,1;\mathbb{S})}$} &
\multicolumn{2}{|c|}{${L^{\infty}(0,1;L^2(\Omega))}$}\\ \hline
\multicolumn{6}{c}{${\tt m}\equiv 1$}\\ \hline
4 & 32 & 1.579e+01 & - & 2.498e+00 & -\\ \hline
16 & 128 & 7.617e+00 & 1.05 & 8.293e-01 & 1.59\\ \hline
64 & 512 & 3.615e+00 & 1.08 & 2.400e-01 & 1.79\\ \hline
256 & 2048 & 1.736e+00 & 1.06 & 6.497e-02 & 1.89\\ \hline
1024 & 8192 & 8.475e-01 & 1.03 & 1.693e-02 & 1.94\\ \hline
4096 & 32768 & 4.182e-01 & 1.02 & 4.324e-03 & 1.97\\ \hline
\multicolumn{6}{c}{${\tt m}\equiv 2$}\\ \hline
4 & 72 & 4.781e+00 & - & 3.812e-01 & -\\ \hline
16 & 288 & 1.013e+00 & 2.24 & 5.422e-02 & 2.81\\ \hline
64 & 1152 & 2.282e-01 & 2.15 & 7.330e-03 & 2.89\\ \hline
256 & 4608 & 5.480e-02 & 2.06 & 1.240e-03 & 2.56\\ \hline
1024 & 18432 & 1.349e-02 & 2.02 & 1.743e-04 & 2.83\\ \hline
4096 & 73728 & 3.427e-03 & 1.98 & 2.259e-05 & 2.95\\ \hline
\end{tabular}
\end{center}
\caption{Convergence test. Errors and convergence rates under uniform mesh refinement.
A fixed time step of size $.5\times 10^{-3}$
was used.
Linear (above) and quadratic (below) dG in space discretizations.
}
\end{table}

{\bf Advection-dominated test.} With this numerical example we test the robustness of the method in the advection-dominated regime.
We consider once again the domain $\Omega=[-1,1]^2$ subdivided into the  subdomains $\Omega^1$ and $\Omega^2$ already considered in the previous test. In $\Omega$, we solve a scalar equation with no reaction and the
constant diffusion $a_1=10^{-2}$
and the advection field $B_1=(0.5,0.5)$.
On the boundary $\partial\Omega$, we set homogeneous Neumann boundary conditions.
\eqref{transmission}
 The parameters
\[
\mbf{P}=2/10,
\qquad
\Upsilon^{1} =5/6,
\qquad
\Upsilon^{2} =1/6,
\qquad
{\rm R} =6/10,
\]
are used in the  interface conditions~\eqref{transmission}. Thus transfer across the interface is mainly advection-driven. Further, notice that setting a friction coefficient to
less than one models the case in which the interface acts as a filtering wall on the advected quantity; hence a boundary layer in the upwind subdomain in the proximity of the interface is expected.

We solve the problem on a uniform $16\times 16$ mesh using linear elements. Such mesh is not fine enough to resolve the layer forming in the proximity of the interface where the solution is also discontinuous, see Figure~\ref{fig:ex2sol:16}. Nevertheless, the numerical solution is stable, and the expected behaviour of the solution is accurately captured, as we can see by comparison with the solution obtained with a layer resolving $64\times 64$ mesh and shown in Figure~\ref{fig:ex2sol:64}.
\begin{figure}
\centering
\includegraphics[scale=0.12]{figures/ref4_sol_0_10_20_30_40_50_60_70.pdf}
\caption{Advection-dominated test. Snapshots of the solution computed on a uniform $16\times 16$ mesh using linear elements: the initial condition (top-left) followed by the solution at time intervals of $0.5$.}
\label{fig:ex2sol:16}
\end{figure}
\begin{figure}
\centering
\includegraphics[scale=0.2]{figures/ref6_sol_10_20.pdf}
\caption{Advection-dominated test. Snapshots of the solution computed on a uniform $64\times 64$ mesh using linear elements: solution at $t=0.5$ (left) and $t=1$ (right) corresponding to the second and third plots in the first row of Figure~\ref{fig:ex2sol:16}, respectively.}
\label{fig:ex2sol:64}
\end{figure}