% !TEX encoding = IsoLatin
%% Antes de processar este arquivo LaTeX (LaTeX2e) deve
%% verificar que o arquivo TEMA.cls esta na mesma
%% pasta. O arquivo TEMA.cls pode ser obtido do
%% endereco http://tema.sbmac.org.br.

\documentclass{TEMA}

\usepackage[brazil]{babel}      % para texto em Português
%\usepackage[english]{babel}    % para texto em Inglês
 
\usepackage[latin1]{inputenc}   % para acentuação em Português IsoLatin1
%\usepackage[utf8]{inputenc}     % para acentuação em Português UTF-8

% Atenção: O presente arquivo usa a codificação IsoLatin1. Para
% utilizar a codificação UTF-8, converter o arquivo usando seu
% programa favorito, então troque o argumento do inputenc 
% adequadamente.

\usepackage[dvips]{graphics}
\usepackage{subfigure}
\usepackage{graphicx}
\usepackage{epsfig}
\usepackage{hyperref}
\usepackage{framed}
\usepackage{psfrag}
\usepackage{tikz} 
%%%%% PACKAGES INCLUDED %%%%%
\usepackage{amsmath,amssymb}


\newcommand{\B}{{\tt\symbol{92}}}
\newcommand{\til}{{\tt\symbol{126}}}
\newcommand{\chap}{{\tt\symbol{94}}}
\newcommand{\agud}{{\tt\symbol{13}}}
\newcommand{\crav}{{\tt\symbol{18}}}

%%%%%%%%% DEFINITIONS %%%%%%%%%%%%%%%%%%%
\newcommand{\bv}{{\bf v }}
\newcommand{\Ec}{{\cal{E}}}
\newcommand{\Tc}{{\cal{T}}_h}
\newcommand{\Kc}{{\cal{K}}}
\newcommand{\N}{{\cal{N}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

%********************************************************
\title{
     The influence of velocity field approximations in tracer injection processes%
%%     \thanks{Agradecimentos por auxílio. Baseado no documento
%%     proposto por A. Sri Ranga e E. Wendland.} 
}

%\author{
%     Y.~R.~N\'u\~nez%
%     \thanks{yoisellrn@id.uff.br},
%     Instituto de Ciências Exatas,
%     UFF, Universidade Federal Fluminense, 
%     Rua Desembargador Elllis Hermydio Figueira, 783, 
%     27255-125, Volta Redonda, RJ, Brasil.
%     \\ \\
%     C.~O.~Faria%
%     \thanks{cofaria@ime.uerj.br},
%     Instituto de Matemática e Estatística,
%     UERJ, Universidade do Estado do Rio de Janeiro,
%     Rua São Francisco Xavier, 524, 
%     20550-900, Rio de Janeiro, RJ, Brasil. 
%     \\ \\
%     S.~M.~C.~Malta, A.~F.~D.~Loula
%     \thanks{\{smcm, aloc\}@lncc.br},
%     LNCC-MCTIC, Laboratório Nacional de Computação Científica,
%     Av. Getúlio Vargas, 333, 
%     25651-075, Petrópolis, RJ, Brasil. 
%}

\criartitulo

\runningheads{ % N\'u\~nez, Faria,  Malta and Loula
}%
{The influence of velocity field approximations in tracer injection processes}

\begin{abstract}
{\bf Abstract}  Although the concentration is the most important variable in tracer injection processes, an efficient and accurate velocity field approximation is crucial to obtain a good physical behaviour for the problem.  In this paper we analyse a Stabilized Dual Hybrid Mixed (SDHM)  method  to solve  the  Darcy's system in the velocity and pressure variables that involves the conservation of mass and Darcy's law.  This approach  is locally conservative, free of compromise between the finite element approximation spaces and capable of dealing with heterogeneous media with discontinuous properties. The tracer concentration is solved via a combination of the Streamline Upwind Petrov-Galerkin (SUPG) method in space with an implicit finite difference scheme in time. We also employ a semi-analytical approach (Abbaszadeh-Dehghani analytical solution) to integrate the transport equation. A  numerical comparative study using the SDHM formulation, the Galerkin method and a post-processing technique to calculate the velocity field in combination  with those concentration approximation methodologies are presented. In all comparisons, the SDHM formulation appears as the most efficient, accurate and almost free of spurious oscillations. 

{\bf Keywords}. Miscible displacements, Hybridized method, Oil reservoir simulations

\end{abstract}


%********************************************************
\section{Introduction}

Information extracted from tracer breakthrough profiles at production wells plays an important role in reservoir engineering both in the characterization of reservoir heterogeneities as in the project of recovery techniques. Those profiles can be detected either experimentally or via the solution of a mathematical model, which describes the transport of substances  through a porous medium. Although the concentration is the variable of most interest, approximation of the velocity field is crucial, since it is responsible for the flow displacement.

Standard Mixed Finite Element (MFE) methods have been extensively employed
during the last decades to solve Darcy's system \cite{raviart-thomas:1977,bdm:1985,brezzi-fortin:1991,douglas:1983,darlow:1984}. The main idea of these methods is simultaneously approximate pressure and velocity by using different spaces for each variable. This leads to a compatibility condition between the approximation spaces (the LBB condition \cite{brezzi:1974}), and thus restricts the choice of stable finite element spaces. However, these choices are usually unstable with standard dual mixed formulation, as was illustrated in \cite{Correa:2008}. A very common example of stable mixed method is given by the Raviart-Thomas spaces \cite{raviart-thomas:1977}. To overcome the  compatibility conditions, some stabilized mixed finite element methods were proposed in \cite{douglas:1989,Murad:1995,masud:2002,Barrenechea:2007,Correa:2007,Correa:2008,Iury-tese:2015,Barrios:2015}. In general, these stabilized formulations  use continuous Lagrangian finite element spaces and can be
successfully employed in simulating Darcy flows in homogeneous porous media. However, continuous interpolations are not appropriate to
heterogeneous porous media with discontinuous properties. This is due to the fact that on the interfaces of the discontinuities the normal component of Darcy velocity must be continuous (mass conservation) but the
tangential component is discontinuous. Formulations based on continuous Lagrangian interpolation for velocity fail to represent the tangential discontinuity, producing inaccurate approximations and spurious oscillations. Therefore, based on hybridization techniques the Stabilized Dual Hybrid Mixed (SDHM) method that combines advantages of the Discontinuous Galerkin (DG) methods \cite{riviere:2008} with reduced computational was developed  for calculating accurate velocity fields to miscible displacements in homogeneous and heterogeneous porous media \cite{yoisell:2014, tracer-ijsmpi:2012, nunez:2017}. The SDHM formulation 
consists of a set of local problems defined at the element level coupled to a global system for the Lagrange multipliers. Stabilization terms are added to generate a stable and adjoint consistent formulation allowing greater flexibility in the choice of the approximation spaces. The Lagrange multiplier is identified as the pressure trace on the element interface, which is a natural choice.

In \cite{tracer-ijsmpi:2012} the authors have shown that for regular solutions the SDHM method leads to optimum rates of convergence for the velocity and pressure fields, even when same order interpolations are employed. In addition, the mixture transport in homogeneous and heterogeneous media were recently analysed \cite{nunez:2017} by studying the influence of the mobility ratio and the permeability field variation. In these cases, the SDHM method was applied in the approximation of the velocity field coupled to the SUPG method in the calculation of the concentration field. The good performance of this methodology was verified through the presented numerical simulations. 

In this paper, following the works mentioned in the above paragraph, we study  the efficiency and robustness of the SDHM comparing with the usual Galerkin method and a post-processing technique \cite{sandra:2000,sandra:1998}. To do this, we analyse the influence of the velocity field approximations on the behaviour of the prediction of tracer concentration at the producer well.  The concentration approximation is obtained via the SUPG method combined with an implicit finite difference scheme, as in \cite{tracer-ijsmpi:2012, nunez:2017}. In order to validate the numerical results, a semi-analytical methodology, which consists of combining the analytical solution for the concentration given by Abbaszadeh-Dehghani \cite{Abbaszadeh-Dehghani:1982} with an approximation for the velocity field is also provided.

The outline of this work is as follows. Section \ref{S2} considers the mathematical model of the problem. The stabilized dual hybrid mixed (SDHM) method for the Darcy's system is presented in Section \ref{S3}. Tracer concentration approximations based on the SUPG approach and on the semi-analytical methodology are discussed in Section \ref{S4}, as well as its combination with the usual Galerkin method and a post-processing technique. Numerical experiments are reported in Section  \ref{S5}. Finally, some conclusions are given in Section \ref{S6}.

\section{The model problem}\label{S2}
\vspace{-2pt}

The equations governing tracer flows in miscible displacements \cite{peaceman:1977} can be described by the first order system of partial differential equations in the Darcy velocity $\mathbf{u}=\mathbf{u}(\mathbf{x},t)$ and the pressure $p = p(\mathbf{x})$ coming from the mass conservation of the mixture and  the Darcy's law:
%
\begin{eqnarray}
\operatorname{div}\mathbf{u}=\ f&\textrm{in}&~~\Omega~\times(0,T)~, \label{eq:cmDarcy}\\
\mathbf{u}=-\mathbb{K}\nabla p  &\textrm{in}&~~\Omega~\times(0,T)~,\label{eq:leiDarcy}
\end{eqnarray}
%
combined with a convection dominated diffusion-convection equation, which expresses the conservation of mass of the injected fluid (the tracer concentration equation), $c=c(\mathbf{x},t)$, gives by
%
\begin{equation}
\phi\frac{\partial c}{\partial t}+\mathbf{u}\cdot\nabla c-\operatorname{div}
(\mathbb{D}\nabla c)=g\text{ \ in }\Omega\times(0,T),\label{eqc}
\end{equation}
\noindent
%
with the following initial and boundary conditions
%
\begin{eqnarray}
c(\mathbf{x},0)   =\ c_{0}(\mathbf{x})&\text{ \ in }&\Omega,\label{eq:cic}\\
\mathbb{D}\nabla c\cdot \mathbf{n}   =0&\text{ \ on }&\partial \Omega\times(0,T),\label{eq:cfc} \\
\mathbf{u}\cdot \mathbf{n} = 0 &\text{ \ on }&\partial \Omega\times(0,T).\label{eq:ufc}
\end{eqnarray}
%
In (\ref{eq:cmDarcy}-\ref{eq:ufc}) $\mathbb{K}=\mathbb{K}(\mathbf{x})$ is the permeability tensor and $\phi$ is the porosity. Functions $f$ and $g$ are the source terms, $\mathbf{n}$ is the exterior normal to $\partial \Omega$ and
%
\begin{equation}
\mathbb{D}=\mathbb{D}(\mathbf{u})=\big(\alpha_{mol}+ \alpha_t \left\vert \mathbf{u}\right\vert\big)\mathbb{I}+
\frac{\alpha_l-\alpha_t}{\left\vert \mathbf{u}\right\vert }\mathbf{u\otimes u} \label{eq:ddsp}
\end{equation}
%
is the dispersion-diffusion tensor \cite{peaceman:1977} 
where $\alpha_{mol},$ $\alpha_{l}$ and $\alpha_{t}$ are the molecular diffusion, longitudinal and transversal dispersion coefficients,
respectively, with the domain $\Omega \subset \mathbb{R}^{n}$, $\partial \Omega$ its boundary and $T \in \mathbb{R}, T>0$.

\vspace{-6pt}
\section{The Hybrid Method For Darcy Flow Problem}\label{S3}
\vspace{-2pt}


Let $\Tc=\left\{  \Kc\right\}  $ be a regular finite element mesh  on the domain $\Omega$ $\left( \Omega=\bigcup\limits_{\Kc}\Kc\right)$. The set of all edges of all elements $\Kc$ is $\Ec_h = \{ e: e  \ {\rm is \ an \ edge \ of} \ \Kc \ {\rm for \ all} \ \Kc \in \Tc \}$ with $\Ec_h^0$ denoting the set of the interior edges. To introduce the hybrid formulation we first multiply equations (\ref{eq:cmDarcy})-(\ref{eq:leiDarcy}) by their respective weighting functions and integrate by parts on each element $\mathcal{K}$, getting the following local weak form
\begin{align}
( \mathbb{K}^{-1} \mathbf{u} , \mathbf{v})_{\mathcal{K}} -
( p , \operatorname{div} \mathbf{v} )_{\mathcal{K}} +  \int_{\partial \mathcal{K}} p (\mathbf{v} \cdot \mathbf{n}_{\Kc}) d s &= 0,
\ \forall \mathbf{v} \in \mathcal{U}_\mathcal{K}, \label{eq:leiDarcy_int.por.partes}\\
- ( \operatorname{div} \mathbf{u} ,  q )_{\mathcal{K}} + (f, q)_{\mathcal{K}} &= 0,
\ \forall q \in \mathcal{Q}_{\mathcal{K}}. \label{eq:conserv.massa_forma.fraca}
\end{align}
%
$\mathcal{U}_{\Kc}   =\{  \mathbf{v}\in\left[  L^{2}\left(  \Kc\right)  \right]^{2}, \operatorname{div}\mathbf{v}\in L^{2}\left(  \Kc\right), \  \forall \Kc\in\Tc \}$ and $\mathcal{Q}_{\Kc}   =\left\{  q\in L^{2}\left(\Kc\right), \forall \Kc\in\Tc\right\}$ are the local function spaces on each element $\mathcal{K}$. For each element $\Kc\in \Tc$, $(w,v)_\mathcal{K} = \int_\mathcal{K} w \ v d\mathcal{K}$ denotes the usual $L^2(\Kc)$ inner product, $ \partial \mathcal{K}$ is set of all edges of element $\mathcal{K}$ and $\mathbf{n}_{\Kc}$ is the external normal to $\partial \Kc$.
%%
For a given $ p = \bar{p}$ on $ \partial \mathcal{K}$, we can solve the set of local problems:

\noindent For each $ \mathcal{K} \in  \mathcal{T}_h$, find $\{\mathbf{u},p\} \in \mathcal{U}_{\mathcal{K}} 
\times \mathcal{Q}_{\mathcal{K}}$, $\forall \ \{\mathbf{v},q\} \ \in \ \mathcal{U}_{\mathcal{K}} \times \mathcal{Q}_{\mathcal{K}}$ such that
\begin{align*}
(\mathbb{K}^{-1} \mathbf{u},\mathbf{v})_{\mathcal{K}} - (p,\operatorname{div}\mathbf{v})_{\mathcal{K}}
-(\operatorname{div} \mathbf{u},q)_{\mathcal{K}}
 & = - (f,q)_{\mathcal{K}} 
 - \int_{\partial \mathcal{K}} \bar{p} (\mathbf{v} \cdot \mathbf{n}_{\Kc}) d s.
\end{align*}
%
Following the ideas of Arnold et al. \cite{arnold:1985}, an approximation for the pressure trace, $\bar{p}$, can be obtained by solving a global problem associated with the dual hybrid mixed formulation. Defining \ the \ function \ spaces $\mathcal{M} = \{ \mu \in L^2(e), \ \forall e \in \Ec_h \}$, $\ \mathcal{U} =  \prod_{\Kc} \mathcal{U}_{\Kc}$ and $\mathcal{Q}=\prod_{\Kc} \mathcal{Q}_{\Kc}$ the dual hybrid formulation  consists in:

\noindent Find  $\mathbf{u} \in \mathcal{U} $, $p \in \mathcal{Q}$ and  $\lambda \in \mathcal{M}$,  such that
\begin{align}
\displaystyle \sum_{\mathcal{K} \in \mathcal{T}_h} \bigg[ (\mathbb{K}^{-1} \mathbf{u} , \mathbf{v})_{\mathcal{K}}
- ( p ,\operatorname{div} \mathbf{v})_{\mathcal{K}} + \int_{\partial \mathcal{K}} \lambda (\mathbf{v} \cdot \mathbf{n}_{\Kc}) d s \bigg]&= 0,
\forall \mathbf{v} \in \mathcal{U}, \label{eq:gv}\\
\displaystyle \sum_{\mathcal{K} \in \mathcal{T}_h} \big[ - ( \operatorname{div} \mathbf{u},q)_{\mathcal{K}} + (f,q)_{\Kc} \big]&= 0,
\forall q \in \mathcal{Q}, \label{eq:gq}\\
\displaystyle \sum_{\mathcal{K} \in \mathcal{T}_h} \int_{\partial \mathcal{K}} \mu (\mathbf{u} \cdot \mathbf{n}_{\Kc}) d s & = 0,
\forall \mu \in \mathcal{M}. \label{eq:gm}
\end{align}
The Lagrange multiplier  $\lambda$ is identified with trace of the pressure on the all edges of the elements $\Kc$, $\lambda = \bar{p}$,
unlike the classical  primal hybrid formulation of Raviart-Thomas \cite{raviart-thomas:1977} and the stabilized formulation of Ewing et al. \cite{ewing-wang-yang:2003}, where the multiplier is identified with the flux. The third equation, (\ref{eq:gm}), weakly imposes the continuity of the normal component of the velocity field (flux continuity) and the flux boundary condition $\mathbf{u} \cdot \mathbf{n} = 0$ on $\partial \Omega$.

\subsection{The SDHM Formulation}\label{S3b}
To generate a stable and adjoint consistent formulation, allowing greater flexibility in the choice of the finite element approximation spaces for velocity and pressure fields and the Lagrange multiplier, we add to the system (\ref{eq:gv})-(\ref{eq:gm}) inner stabilization terms associated with least square residual forms, coming from the mass balance, the Darcy's law and the curl of Darcy's law, as in \cite{Correa:2007}. Furthermore, we add a stabilization term for the multiplier according to \cite{arnold:2002} obtaining the Stabilized Dual Hybrid Mixed (SDHM) formulation:\\
\noindent Find $\mathbf{u} \in \mathcal{U}$, $p \in \mathcal{Q}$ and $\lambda \in \mathcal{M}$ 
$\forall \mathbf{v} \in \mathcal{U}$, $q \in \mathcal{Q}$ and $\mu \in \mathcal{M}$, such that
%%
\begin{align}
\sum_{\mathcal{K} \in \mathcal{T}_h} & \bigg[ (\mathbb{K}^{-1} \mathbf{u} , \mathbf{v})_{\mathcal{K}}
- ( p ,\operatorname{div} \mathbf{v})_{\mathcal{K}}
+ \int_{\partial \mathcal{K}} \lambda (\mathbf{v} \cdot \mathbf{n}_{\Kc}) ds
+ \frac{1}{2} (\|\mathbb{K}^{-1}\|_\infty (\operatorname{div} \mathbf{u} - f), \operatorname{div} \mathbf{v})_{\mathcal{K}} \nonumber \\
&- \frac{1}{2}( \mathbb{K}^{-1} \mathbf{u} + \nabla p, \mathbf{v})_{\mathcal{K}}
+ \frac{1}{2} (\|\mathbb{K}\|_\infty (\nabla \times \mathbb{K}^{-1} \mathbf{u}), \nabla \times \mathbb{K}^{-1} \mathbf{v})_{\mathcal{K}} \bigg] = 0, \label{eq:sv}\\
\sum_{\mathcal{K} \in \mathcal{T}_h} & \bigg[ \! - \! ( \operatorname{div} \mathbf{u},q)_{\mathcal{K}} + (f,q)_{\Kc}
- \frac{1}{2} (\mathbb{K}^{-1} \mathbf{u} +  \nabla p, \mathbb{K}\nabla q)_{\mathcal{K}}
+ \! \int_{\partial \mathcal{K}} \! \|\mathbb{K}\|_\infty \beta (p - \lambda) q  ds \bigg] = 0, \label{eq:sq}\\
\sum_{\mathcal{K} \in \mathcal{T}_h} & \bigg[ \int_{\partial \mathcal{K}} \mu(\mathbf{u} \cdot 
\mathbf{n}_{\Kc}) d s
+ \int_{\partial \mathcal{K}} \|\mathbb{K}\|_\infty \beta (\lambda - p) \mu ds  \bigg] = 0. \label{eq:sm}
\end{align}
\noindent
%
where $\mathcal{U} = \prod_{\mathcal{K}}( H^1(\mathcal{K})\times H^1(\mathcal{K}))$,
$\mathcal{Q} = \prod_{\mathcal{K}} H^1(\mathcal{K})$,
$\mathcal{M} = \{ \mu \in L^2(e), \ \forall e \in \mathcal{E}^{0}_h \}$,
$\nabla \times$ is the curl operator and
$\beta=\frac{\beta_{0}}{h}$ is the stabilization parameter associated with the Lagrange multiplier $\lambda$,
and  $\beta_{0} \in \mathbb{R}$ is independent of $h$.

\subsection{ Finite Element Approximations}
%

Let $\mathcal{U}^m_h = \{ \bv \in \mathcal{U}_{\delta}: \ \bv|_{\Kc} \in \mathcal{R}^m\times\mathcal{R}^m \
\ \forall \Kc \in \Tc \}$, $\mathcal{Q}^l_h = \{ q \in \mathcal{Q}_{\delta}: \ q|_{\Kc} \in \mathcal{R}^l \ \ 
\forall \Kc \in \Tc \}$
and $\mathcal{M}^s_h = \{ \mu \in \mathcal{M}: \ \mu|_{e} \in \mathcal{P}^s \ \ \forall e\in \Ec_h^0 \}$ be the discontinuous Lagrangian finite element spaces where $\mathcal{R}^r$ is the polynomial set with degree less than or equal to  $r$ if $\Kc$ is a triangle, or less than or equal  to $r$ in each cartesian variable if $\Kc$ is a quadrilateral ($r=l$ or $m$), and $\mathcal{P}^s$ is the polynomial set  of degree less than or equal to $s$ on each edge $e$. We can now present a finite element approximation for the stabilized dual hybrid mixed formulation introduced in the last section. Considering that $\{\mathbf{u}_h,p_h\}$, belonging to 
the broken function spaces, are defined independently on  each element $\mathcal{K} \in \Tc$, we observe that system (\ref{eq:sv})-(\ref{eq:sm}) can be split into a set of local problems defined on each element $\mathcal{K}$ coupled to the global problem defined on $\Ec_h$, as follow:\\
\noindent
{\bf \it Local problems:} Find  $\mathbf{{u}}_h \in \mathcal{U}^m_h, \ p_h \in \mathcal{Q}^l_h$,
for each $\mathcal{K} \in \mathcal{T}_h$, $\forall  \mathbf{v}_h \in \mathcal{U}^m_h$ and $\forall q_h \in \mathcal{Q}^l_h$ such that
\begin{align}
&(\mathbb{K}^{-1} \mathbf{u}_h , \mathbf{v}_h)_{\mathcal{K}}
- ( p_h ,\operatorname{div} \mathbf{v}_h)_{\mathcal{K}}
+ \int_{\partial \mathcal{K}} \lambda_h (\mathbf{v}_h \cdot \mathbf{n}_{\Kc}) d s
- \frac{1}{2} ( \mathbb{K}^{-1} \mathbf{u}_h + \nabla p_h, \mathbf{v}_h)_{\mathcal{K}} \nonumber \\
&+\frac{1}{2} (\|\mathbb{K}^{-1}\|_\infty ( \ \operatorname{div} \mathbf{u}_h - f),\ \operatorname{div} \mathbf{v}_h)_{\mathcal{K}} 
+ \frac{1}{2} (\|\mathbb{K}\|_\infty \ (\nabla \times \mathbb{K}^{-1} \mathbf{u}_h), \nabla \times \mathbb{K}^{-1} \mathbf{v}_h)_{\mathcal{K}}  =   0, \label{eq:svs}\\
&- \!( \!\operatorname{div} \mathbf{u}_h,q_h)_{\mathcal{K}} \! + \! (f,q_h)_{\Kc}
\!-\! \frac{1}{2} (\mathbb{K}^{-1} \!\mathbf{u}_h\! +\! \nabla p_h, \mathbb{K}\nabla q_h)_{\mathcal{K}}
+ \! \int_{\partial \mathcal{K}} \!\|\mathbb{K}\|_\infty \beta (p_h \!-\! \lambda_h) q_h ds \!= \!0,  \label{eq:sqs}
\end{align}
%%
\noindent
{\bf \it Global Problem:} Find $ \lambda_h\in \mathcal{M}^s_h $, $\forall \mu_h \in \mathcal{M}^s_h$ such that
\begin{equation}
 \sum_{\mathcal{K} \in \mathcal{T}_h} \bigg[ \int_{\partial \mathcal{K}} \ \mu_h (\mathbf{{u}}_h \cdot \mathbf{n}_{\Kc}) ds
+ \int_{\partial \mathcal{K}} \|\mathbb{K}\|_\infty  \beta (\lambda_h - p_h) \mu_h ds  \bigg] = 0. \label{eq:sms}
\end{equation}
%

We observe that the  SDHM method is consistent and provides optimal rates of  convergence 
\cite{tracer-ijsmpi:2012}, ensuring flexibility in the choice of the approximation spaces and interpolation 
functions, unlike the classical dual formulation which the stability is restricted to appropriate choices of 
the finite elements spaces, such as Raviart-Thomas \cite{raviart-thomas:1977} and Brezzi-Douglas-Marini (BDM) 
 \cite{bdm:1985} that requiring compromises between the approximation spaces.
Moreover, the  SDHM method is locally conservative for equal order approximations of all fields ($l = m = s $) 
and stable for any value of the edge stabilization parameter $\beta$,
including $ \beta = 0$ \cite{tracer-ijsmpi:2012,yoisell:2014, nunez:2017}. The choose of the multiplier 
as the trace of pressure is crucial to assure that the local problems
(\ref{eq:svs})-(\ref{eq:sqs}) are solvable for the variables
$\{\mathbf{u}_h, p_h\}|_{\mathcal{K}} \in \mathcal{U}^m_h \times \mathcal{Q}^l_h$
as a function of the multiplier $\lambda_h$. 


\subsubsection{Solver Strategy} \label{decp}

 Here we define on each element the  operators
\begin{align*}
\displaystyle a_{\Kc}([\mathbf{{u}}_h,p_h];[\mathbf{v}_h,q_h]) &= 
(\mathbb{K}^{-1} \mathbf{{u}}_h ,  \mathbf{v}_h)_{\mathcal{K}}
- ( p_h ,\operatorname{div} \mathbf{v}_h)_{\mathcal{K}} 
+ \frac{1}{2} \|\mathbb{K}^{-1}\|_\infty (\operatorname{div} \mathbf{{u}}_h, \operatorname{div} \mathbf{v}_h)_{\mathcal{K}}\\
&- ( \operatorname{div} \mathbf{u}_h,q_h)_{\mathcal{K}} -\frac{1}{2} ( \mathbb{K}^{-1}\mathbf{{u}}_h + \nabla p_h, \mathbf{v}_h + \mathbb{K}\nabla q_h)_{\mathcal{K}} \\
&+\frac{1}{2} (\|\mathbb{K}\|_\infty \nabla \times (\mathbb{K}^{-1} \mathbf{{u}}_h), \nabla (\times \mathbb{K}^{-1} \mathbf{v}_h))_{\mathcal{K}},\\
\displaystyle b_{\Kc}(\lambda_h;[\mathbf{v}_h,q_h]) &= 
\int_{\partial \mathcal{K}} \lambda_h (\mathbf{v}_h \cdot \mathbf{n}) ds
+ \int_{\partial \mathcal{K}} \beta \ \|\mathbb{K}\|_\infty \ \lambda_h q_h ds,\\
 \displaystyle c_{\Kc}(\lambda_h,\mu_h) &= \int_{\partial \mathcal{K}} \|\mathbb{K}\|_\infty \ \beta \lambda_h \mu_h ds, 
\end{align*}
and the functional 
\begin{align*}
\displaystyle f_{\Kc}([\mathbf{v}_h,q_h]) = (f,q_h)_{\mathcal{K}}
+  \frac{1}{2} \|\mathbb{K}^{-1}\|_\infty \ (f, \operatorname{div} \mathbf{v}_h)_{\mathcal{K}}. \quad \qquad \qquad \qquad
\end{align*}
Then, the SDHM method is now reformulated as
%

\noindent Find  $\mathbf{{u}}_h \in \mathcal{U}^m_h, \ p_h \in \mathcal{Q}^l_h$,
for each $\Kc \in \mathcal{T}_h$ and $ \lambda_h\in \mathcal{M}^s_h $
such that
%
\begin{align}
 a_{\Kc}([\mathbf{{u}}_h,p_h];[\mathbf{v}_h,q_h]) + b_{\Kc}(\lambda_h;[\mathbf{v}_h,q_h]) &= f_{\Kc}([\mathbf{v}_h,q_h]), \label{eq:localp}\\
 \sum_{\mathcal{K} \in \Tc} b^{T}_{\Kc}([\mathbf{{u}}_h,p_h],\mu_h) + \sum_{\mathcal{K} \in \Tc} c_{\Kc}(\lambda_h,\mu_h) &= 0, \label{eq:globalp}
\end{align}
%
$\forall  [\mathbf{v}_h,q_h] \in \mathcal{U}^m_h \times \mathcal{Q}^l_h$ and $\forall \mu_h \in \mathcal{M}^s_h$.   Considering $\mathbf{A}_{\Kc}$, $\mathbf{B}_{\Kc}$ and $\mathbf{C}_{\Kc}$ 
matrices generated respectively by
the local operators $a_{\Kc}(\cdotp,\cdotp)$, $b_{\Kc}(\cdotp,\cdotp)$ and $c_{\Kc}(\cdotp,\cdotp)$
and $\mathbf{F}_{\Kc}$ the vector given by $f_{\Kc}(\cdotp)$,
we can rewrite (\ref{eq:localp})-(\ref{eq:globalp}) in the following matrix form,
%
\begin{align}
 \mathbf{A}_{\Kc} \mathbf{U} + \mathbf{B}_{\Kc} \Lambda &= \mathbf{F}_{\Kc}, \ \forall \Kc \in \Tc \label{localp-matricial}\\
 \sum_{\mathcal{K} \in \Tc} \mathbf{B}^{T}_{\Kc} \mathbf{U} + \sum_{\mathcal{K} \in \Tc} \mathbf{C}_{\Kc} \Lambda &= \mathbf{0},
 \label{globalp-matricial}
\end{align}
%
Using the static condensation strategy, and given that the  local matrix $\mathbf{A}_\Kc$ is positive definite, 
we have from (\ref{localp-matricial}) that $\mathbf{U}=\{u_h, p_h\}$ can be expressing in terms of $\Lambda = \{ \lambda_h \}$ as
\begin{equation}\label{solucao-localp}
\mathbf{U}=\mathbf{A}_\Kc^{-1}(\mathbf{F}_\Kc-\mathbf{B}_\Kc \Lambda).
\end{equation}
Then, replacing (\ref{solucao-localp}) in (\ref{globalp-matricial}), 
we obtain the global system only for $\Lambda = \{ \lambda_h \}$
\begin{align}\label{mf-tot}
 \sum_{\mathcal{K} \in \Tc}(\mathbf{C}_{\Kc} - \mathbf{B}^{T}_{\Kc}\mathbf{A}^{-1}_{\Kc}\mathbf{B}_{\Kc})\Lambda
 = -\sum_{\mathcal{K} \in \Tc} \mathbf{B}^{T}_{\Kc}\mathbf{A}^{-1}_{\Kc}\mathbf{F}_{\Kc}.
\end{align}
After solving (\ref{mf-tot}), the vector $\mathbf{U}$ is
obtained from (\ref{solucao-localp}).
A great advantage of this methodology is the size reduction of
the overall system, now involving only the degrees of freedom associated with the multipliers $\boldsymbol{\Lambda}$, 
leading to a reduced computational cost, since the time needed to solve all local problems is negligible 
compared to the time to solve the global system.

\vspace{-6pt}
\section{Tracer Injection Approximations}\label{S4}
\vspace{-2pt}

\subsection{Fully discrete SUPG approximation}

A common approach to transient problems is based on fully discrete formulations obtained by combining finite difference approximations in time with finite element methods in space.
According to the Rothe method (or horizontal method of lines) \cite{harari:2004}
of first discretizing in time and then in space on each discrete time level,
we choose the partition $I_{\Delta } = \{0 = t^{0} < t^{1} < \ldots < t^{N} = T\}$
of the interval $I = [0,T]$ with $\Delta t^{n} = t^{n} - t^{n-1}$ and $\Delta t = \max_{n} \Delta t^{n}$. 
Then, using a backward finite difference scheme to  approximate the time derivative in equation
(\ref{eqc}), we have the following sequentially implicit time-stepping algorithm:  
for  $n=0,1,\ldots,N-1$, given $f$,
$g$ and $c_{0}(\mathbf{x})$, find $c^{n+1}$  satisfying
\begin{align}
\phi\frac{c^{n+1} - c^{n}}{\Delta t}+\mathbf{u}^n\cdot\nabla c^{n+1}-\operatorname{div}
(\mathbb{D}(\mathbf{u}^n)\nabla c^{n+1}) &= g^{n+1}  \text{ \ in }\Omega,\label{eq:eqc3}\\
\mathbb{D}(\mathbf{u}^n)\nabla c^{n+1}\cdot \mathbf{n} &=0 \text{ \ on }\partial \Omega,\label{eq:cfc3}
\end{align}
with 
%
\begin{align}
\mathbf{u}^n &= - \mathbb{K}\nabla p^{n} \text{ \ in} ~~\Omega, \label{eq:leiDarcy3}\\
\operatorname{div}\mathbf{u}^n &= f^n \text{ \ on} \ \Omega, \label{eq:cmDarcy3}
\end{align}
%
%
and $\mathbf{u}^n \cdot \mathbf{n} = \ 0  \text { \ on} \ \partial\Omega$ and  $\int_{\Omega} p^{n} dx = 0, \ t \in (0,T)$.
The concentration is achieved at time $n + 1$ and the velocity and pressure are given at time $n$.
A complete numerical analysis, demonstrating existence and uniqueness of solution for the above semi-discrete system can be found in \cite{sandra:2000,sandra:1998}. Note that this sequentially implicit method can be written in predictor-corrector form and the original system becomes partially uncoupled and linearized.


For the tracer injection processes, the velocity and pressure approximations
are calculated just once at the beginning of the process, $\mathbf{u}^n \equiv \mathbf{u}$ $\forall n=0,1, \ldots,N-1$ and, then the concentration is obtained for all instants. Hereafter, to simplify the notation, the superscript $n$ in the velocity and pressure approximations  are dropped out.


We combine the semi-discrete approximation (\ref{eq:eqc3})-(\ref{eq:cfc3})
with a stabilized finite element method in space (the SUPG method), and introduce the following
fully discrete approximation for the concentration equation \cite{tracer-ijsmpi:2012,yoisell:2014,nunez:2017}: \\
\noindent For time levels $n=0,1,2,\ldots$,
given $\mathbf{u}_h$ and $c^0$
find $c_h^{n+1} \in X_h^r$, such that
\begin{align}
 B(\mathbf{u}_h; c_h^{n+1}, \eta_h) = F(c_h^n; \eta_h), \qquad \forall \eta_h \in X_h^r,
\end{align}
with
$(c_h^0, \eta_h) = (c^{0}, \eta_h), \ \forall \eta_h \in X_h^r$,
where $X_h^r$ is a continuous Lagrangian finite element space of degree $r \geq 1$; and
\begin{align}
 B&(\mathbf{u}_h; c_h^{n+1}, \eta_h) = (c_h^{n+1}, \eta_h) + \Delta t (\mathbf{u}_h \cdot \nabla c_h^{n+1}, \eta_h)
+ \Delta t (D(\mathbf{u}_h)\nabla c_h^{n+1}, \nabla \eta_h)\nonumber \\
&+ \displaystyle \sum_{\Kc} (c_h^{n+1} + \Delta t \ \mathbf{u}_h \cdot \nabla c_h^{n+1} , \delta_K \mathbf{u}_h \cdot \nabla \eta_h)_{\Kc}\nonumber\\
&+ \displaystyle \sum_{\Kc} ( - \Delta t \operatorname{div}(D(\mathbf{u}_h)\nabla c_h^{n+1}),
\delta_K \mathbf{u}_h \cdot \nabla \eta_h)_{\Kc},\label{eq:bc}\\
 F&(c_h^n, \eta_h) = (c_h^n + \Delta t g^{n+1}, \eta_h)
+ \displaystyle \sum_{\Kc} (c_h^n + \Delta t g^{n+1}, \delta_K \mathbf{u}_h \cdot \nabla \eta_h)_{\Kc}. \label{eq:fc}
\end{align}
%
The stabilization parameter $\delta_K$ is defined on each $\Kc \in \mathcal{T}_h$ by \cite{brooks:1982, sandra:2000,sandra:1998} as
\begin{equation}
\delta_{\Kc}=\left\{
\begin{array}[c]{cl}
\dfrac{h_{\Kc}}{2 \| \mathbf{u_h}\|_{L^{\infty}(\Kc)}},& \textrm{Pe}_{\Kc} \geq 1,\\
0,& 0<\textrm{Pe}_{\Kc} <1,
\end{array}
\right.
\end{equation}
where
\begin{equation}
\textrm{Pe}_{\Kc} = \frac{m_{\Kc} \|\mathbf{u}_h\|_{L^{\infty}(\Kc)}h_{\Kc}}{D^2_{\Kc}/d_{\Kc}}
\end{equation}
is the mesh-dependent P\'eclet number, 
$\displaystyle{ d_{\Kc} = \alpha_{\textrm{mol}} + \alpha_t \inf_{x \in \Kc} | \mathbf{u}_h(\mathbf{x})|}$, 
 $\displaystyle{ m_{\Kc} = \frac{2}{3} \min (\frac{1}{2},c_{\textrm{inv}})}$,
 $\displaystyle{ D_{\Kc} = \sqrt{2(\alpha_{\textrm{mol}} + \alpha_l \|
 \mathbf{u}_h\|_{L^{\infty}(\Kc)})^2+ 2(3 \alpha_l - 2 \alpha_t )^2
 \|\nabla \mathbf{u}_h\|_{L^{\infty}(\Kc)}^2 h_{\Kc}^2 c_{\textrm{inv}}}}$,
where $c_{\textrm{inv}}$ is the typical inverse constant of finite element spaces.
The terms in the right-hand side of (\ref{eq:bc})--(\ref{eq:fc}) multiplied by $\delta_{\Kc}$ are responsible for the additional stability of this method \cite{brooks:1982}.
In the numerical experiments this proposal will be identified as {\it
the SUPG approach}.

\subsection{ A semi-analytical methodology}

 Following \cite{datta-gupta:1995}, we present a semi-analytical methodology  which
combines  an analytical solution for the concentration given by \cite{Abbaszadeh-Dehghani:1982}
with a finite element approximation for the velocity field.
%
According to \cite{Abbaszadeh-Dehghani:1982}, neglecting transversal flux dispersion ($\alpha_{mol}=\alpha_t=0$), the tracer concentration $c(t)$ at the producer well
is expressed by
 %
\begin{align}
c(t) = \frac{\hat{c}V_{tr}}{2q_{t}(\pi \alpha_l)^{1/2}} \sum^{N_{sl}}_{n=1}
\frac{1}{I^{1/2}_n} \exp\left( -\frac{(t_{btn}-t)^2}{4\alpha_l I_n}\right), \label{eq:solucao-analitica-conc}
\end{align}
%
where $N_{sl}$ is the number of streamlines arriving at the production well, $\hat{c}$ is the injection 
tracer concentration, $V_{tr}$ is the injected tracer volume, $q_{t}$ is the total injection rate, $\alpha_l$ 
is the longitudinal dispersion coefficient, $t_{btn}$ is the breakthrough time for the streamline $n$ given by
$\displaystyle 
t_{btn}=\int_{x_{iw}}^{x_{pw}}\frac{ds}{v_{n}} 
$
%
and
%
$
\displaystyle I_n=\int_{x_{iw}}^{x_{pw}}\frac{ds}{{v_{n}}^2}, \label{In-linhas.de.fluxo}
$
%
where ${v}_{n}$ is the Darcy velocity on the streamline $n$, with $x_{iw}$ and $x_{pw}$ the coordinates of the injection and production wells, respectively.

In \cite{Abbaszadeh-Dehghani:1982} the Darcy's velocity ${v}_{n}$ is also obtained analytically for
particular tracer injection problems in homogeneous porous media.
However, from the definitions of $t_{btn}$ and $I_n$ it is clear that the {\it semi-analytical methodology} 
(\ref{eq:solucao-analitica-conc}) can be naturally reformulated replacing $v_n$ by an approximation of the velocity field. In Section \ref{S5} we will observe, through numerical experiments, the influence of the velocity approximations on some tracer transport simulations in homogeneous and heterogeneous porous media 
solving  the Darcy system (2.1)-(2.2) via the SDHM method, the classical Galerkin method and a post-processing technique.

{\bf Classical Galerkin method.} Substituting (\ref{eq:leiDarcy})  into (\ref{eq:cmDarcy}) the Galerkin approximation of the elliptic sub-system reads:
 Find $p_h\in \N_h^r$, such that
\begin{equation}
 (\mathbb{K}^{-1}\nabla p_h ,\nabla\varphi)=(f,\varphi_h), \quad \forall \varphi_h\in \N_h^r, \label{eq:Galerkinpressao}
\end{equation}
with $\N_h^r$ a continuous Lagrangian finite element space.  From the pressure approximation  (\ref{eq:Galerkinpressao}) and
Darcy's law (\ref{eq:leiDarcy}),  it is {\it natural}  to calculate the velocity field giving by
 \begin{equation}
 \mathbf{{u}_{G}} = - \mathbb{K}\nabla p_{h}. \label{eq:Galerkinvelocidade}
 \end{equation}
This approach generates a discontinuous velocity field at the element interfaces which do not satisfy the boundary condition $\mathbf{{u}_{G}}\cdot \mathbf{n} $ and, in addition, has a sub-optimal convergence rate.


{\bf Post-processing.} To approximate the velocity field with improved accuracy we employed a
post-processing technique \cite{sandra:1998}
based on a variational formulation of Darcy's law combined with the residual of the balance equation:\\
For $p_h$ given by (\ref{eq:Galerkinpressao}), find $\mathbf{{u}_{PP}}\in \SS_h^q=\{\mathbf{v}_h\in \N_h^r\times\N_h^r, \
\mathbf{v}_h\cdot \mathbf{n}=0 \ {\rm on} \ \partial\Omega\}$ such that
\begin{equation}
(\mathbb{K}^{-1}\mathbf{{u}_{PP}}+\nabla p_h,\mathbf{v}_h )+\tau(\operatorname{div} \mathbf{{u}_{PP}} - f , \operatorname{div}
 \mathbf{v}_h)=0, \quad \forall \mathbf{v}_h\in \SS_h^q. \label{eq:postprocessing}
\end{equation}
For homogeneous media with regular solution
this post-processing presents  a gain of $\mathbb{O}(h^{1/2})$ \cite{sandra:1998}
compare to the classical Galerkin approximation (\ref{eq:Galerkinvelocidade}). However, it does not
produce accurate solution when applied to heterogeneous porous media as we will observe
in the numerical experiments.

\vspace{-6pt}
\section{Numerical Simulations}\label{S5}
\vspace{-2pt}

To illustrate the performance of the methodology proposed in Sections \ref{S3} and \ref{S4} we present the results of tracer injection simulations in a quarter of a repeated five-spot pattern in two dimensions consisting of a square domain (unit thickness) with side $L = 1000 ${\it ft}. The injector well is located at the lower-left corner ($x = y = 0$) and the producer well at the upper-right corner ($x = y = L$).
The coefficients have been chosen as $\alpha_{mol} = 0.0$,  $\alpha_l = 1.0$ {\it ft},
$\alpha_t = 0.0$ {\it ft}  and the porosity $\phi=0.1$.
The tracer slug is $0.25 \%$ of the porous volume.
These data were taken from \cite{sandra:2000}.

In all simulations we use uniform meshes of $80 \times 80$ bilinear quadrilateral elements to approximate concentration
as well as velocity and pressure.
Equal order interpolation functions
($m=l=r=s=1$) are employed to all variables.
Discontinuous Lagrangian interpolation to the SDHM method
and continuous Lagrangian interpolation to the Galerkin method and post-processing technique are used.
The post-processing stabilization parameter is fixed as $ \tau = 1.0$ and the stabilization
parameter associated with the Lagrange multiplier in
SDHM is $\beta_0 = 0$. The semi-analytical methodology parameters are
 $\hat{c} = 1.0$,  $V_{tr} = 1000$ {\it ft}$^{2}$, $q_{t}=200$ {\it ft}$^{2}$/day and $N_{sl}=399$.

In order to compare the differents approaches described in Section \ref{S3}, we take into account three scenarios: Case 1, Case 2 and Case 3.
In Case 1 the porous medium is  homogeneous with the permeability constant,  $K = 1000$ {\it mD} and the other scenarios consider heterogeneous porous media with subregions having different permeability values (Figure \ref{definicao.casos}).

% Figura 1
\begin{figure}[!htb]
\subfigure[Case 2]{\includegraphics[scale=0.3]{figures/cenario02.eps}} %
\subfigure[Case 3]{\includegraphics[scale=0.3]{figures/cenario03.eps}} %
\caption{Heterogeneous scenarios: (a) Case 2 and (b) Case 3.}
\label{definicao.casos}
\end{figure}

\subsection{Homogeneous Porous Medium}


Figure \ref{historia-concentracao.tracador-meio.hom02}
exhibits comparisons between the results of the SUPG approach and the semi-analytical methodology for the concentration with velocity approximations obtained by the SDHM formulation, the post-processing technique (PP) and the Galerkin method.
For both approximations we see that the
SDHM method has the  closest profile to the analytical solution obtained by
Abbaszadeh-Dehghani \cite{Abbaszadeh-Dehghani:1982} (peak at 0.01)
as showed in Figure \ref{historia-concentracao.tracador-meio.hom02}-(b).
The conservative property of the SDHM method is illustrated in Figure \ref{volume.tracador-meio.hom-SUPG}.
Therefore, we can conclude that the  SDHM formulation is the more stable and provides accurate results in all approaches considered here
for the transport equation.

%========================================================================================================
%Concentração do Traçador no poço de produção - Cenário Homogêneo --> Metodologia SUPG e Metodologia semianalítica
%
%Figura 2
\begin{figure}[!htb]
\centering
\subfigure[SUPG approach]{\includegraphics[height=4.0cm, width=6.1cm]{figures/TimeHistoryCase1SUPG.eps}}
\hspace{0.25cm}
\subfigure[Semi-analytical methodology]{\includegraphics[height=4.0cm, width=6.1cm]{figures/TimeHistoryCase1Semi.eps}}
\caption{Case 1 - Time history of the tracer concentration to different velocity field approximations.}
\label{historia-concentracao.tracador-meio.hom02}
\end{figure}
%========================================================================================================

% Figura 3
% Volume do Traçador - Cenário Homogêneo - SUPG e Semianalítica
\begin{figure}[!htb]
\centering
\subfigure[SUPG approach]{\includegraphics[height=4.0cm,width=6.1cm]{figures/TraceVolumeCase1SUPG.eps}}
\hspace{0.25cm}
\subfigure[Semi-analytical methodology]{\includegraphics[height=4.0cm,width=6.1cm]{figures/TraceVolumeCase1Semi.eps}}
\caption{Case 1 - Tracer volume of the tracer concentration to different velocity field approximations.}
\label{volume.tracador-meio.hom-SUPG}
\end{figure}


%==================================================================
% Aproximação SUPG + {Galerkin,PP, MDHE} --> Diferentes malhas
%
Next, we study the influence of the P\'eclet number (Pe) on the tracer injection processes.
Figure \ref{historia-concentracao.tracador-diferentes.malhas-meio.hom}
shows the time history of the tracer concentration in the producer well with the SUPG approach for P\'eclet numbers corresponding 
to $100$, $50$, $25$ and $12.5$.
As the P\'eclet number decrease the peak concentration value tends to the
Abbaszadeh-Dehghani exact solution. For high P\'eclet numbers ($Pe \geq 50$)
spurious oscillations associated with the use of the SUPG approach are identified. These 
oscillations become less pronounced
when we combine the SUPG approach with the SDHM method, as can be best viewed in the zoom regions plotted in
Figure \ref{historia-concentracao.tracador-diferentes.malhas-meio.hom} (right side).

%Figura 4
\begin{figure}[!htb]
\centering
\subfigure{\includegraphics[height=4.0cm,width=6.0cm]{figures/TimeHistoryGalerkin.eps}} %
\hspace{0.25cm}
\subfigure{\includegraphics[height=4.0cm,width=6.0cm]{figures/TimeHistoryGalerkinZoom.eps}}\\
\subfigure{\includegraphics[height=4.0cm,width=6.0cm]{figures/TimeHistoryPP.eps}} %
\hspace{0.25cm}
\subfigure{\includegraphics[height=4.0cm,width=6.0cm]{figures/TimeHistoryPPZoom.eps}}\\
\subfigure{\includegraphics[height=4.0cm,width=6.0cm]{figures/TimeHistorySDHM.eps}} %
\hspace{0.25cm}
\subfigure{\includegraphics[height=4.0cm,width=6.0cm]{figures/TimeHistorySDHMZoom.eps}}
\caption{Case 1 - Time history of the tracer concentration (left) and a zoom in the spurious oscillations region (right). The P\'eclet influence study.
Galerkin method, post-processing technique and the SDHM method (from top to bottom).}
\label{historia-concentracao.tracador-diferentes.malhas-meio.hom}
\end{figure}


\subsection{Heterogeneous Porous Media}

In this subsection, we consider two problems defined
in the heterogeneous domains plotted  in Figure \ref{definicao.casos} (Case 2 and Case 3).
Figures \ref{maps-concentration-heterogeneo-cenario01-1000dias} and
\ref{maps-concentration-heterogeneo-cenario03-1000dias}
show the tracer concentration maps at $t = 1000$ days when
we use the SUPG approach combined  with the Galerkin method, the post-processing
technique and the SDHM method to $\Delta t = 5$ days.
Similar transport behaviours are observed for all velocity approximations.
%
%Figura 5
%Case 2
\begin{figure}[!htb]
\centering
\subfigure[]{\includegraphics[scale=0.6]{figures/MapsTracerCase2Galerkin.eps}} %
\subfigure[]{\includegraphics[scale=0.6]{figures/MapsTracerCase2PP.eps}} %
\subfigure[]{\includegraphics[scale=0.6]{figures/MapsTracerCase2SDHM.eps}} %
\caption{Case 2- Tracer concentration maps. (a) Galerkin method, (b) post-processing
technique and (c) the SDHM method.}
\label{maps-concentration-heterogeneo-cenario01-1000dias}
\end{figure}
%
%Figura 6
% Case 3
\begin{figure}[!htb]
\centering
\subfigure[]{\includegraphics[scale=0.6]{figures/MapsTracerCase3Galerkin.eps}} %
\subfigure[]{\includegraphics[scale=0.6]{figures/MapsTracerCase3PP.eps}} %
\subfigure[]{\includegraphics[scale=0.6]{figures/MapsTracerCase3SDHM.eps}} %
\caption{Case 3 - Tracer concentration maps. (a) Galerkin method, 
(b) post-processing technique and (c) the SDHM method.}
\label{maps-concentration-heterogeneo-cenario03-1000dias}
\end{figure}
%
%=====================================================================================
% Concentração de traçador no poço produtor e volume de traçador para os Cenários 1-6
%=====================================================================================
%
% Cenário 1
%

The time history of the tracer concentration in the producer well for Case 2 and Case 3 
using the SUPG approach are showed in
Figures \ref{historia-concentracao.tracador-meio.het-cenario01} and \ref{historia-concentracao.tracador-meio.het-cenario03} for the concentration with velocity approximations obtained by the SDHM formulation, the post-processing technique (PP) and the Galerkin method.
Figure \ref{historia-concentracao.tracador-meio.het-cenario01} exhibits two concentration peaks due to the influence of the 
lower permeability region, which act like a barrier to the flow. 
The SDHM formulation produces the highest peak according to the physical expected behavior.
Note that the tracer concentration reaches the producer well and, the peak concentration is achieved, in an earlier
time compared with that observed in Case 1 (see Figure \ref{historia-concentracao.tracador-meio.hom02}).
This is due to the fact that there are preferential
flow paths where the velocity flow  is increased as a result of the presence of a subregion of lower permeability.
%
%
% Figura 7
\begin{figure}[!htb]
\centering
\includegraphics[height=5.0cm,width=8.0cm]{figures/TimeHistoryCase2SUPG.eps}
\caption{Case 2 - Time history of the tracer concentration to different velocity field approximations.}
\label{historia-concentracao.tracador-meio.het-cenario01}
\end{figure}


%

% Cenário 3
%
Next in Figure \ref{historia-concentracao.tracador-meio.het-cenario03} we observe a similar behavior to the homogeneous scenario (Case 1)  but with higher concentration values. This can be better understand in Case 1 and Case 3 results plot in Figure \ref{historia-concentracao.tracador-comparacao-cenario.hom.vs.cenario3}, where we compare the time history of the tracer concentration obtained by the semi-analytical solution combined with the post-procesing technique and the SDHM formulation. The subregions with differents permeabilities have the same domain width. Therefore, they act as a delay to the flow. 
Thus, the concentration reaches the producer well later than the Case 1.


% Figura 8
\begin{figure}[!htb]
\centering
\includegraphics[height=5.0cm,width=8.0cm]{figures/TimeHistoryCase3SUPG.eps}
\caption{Case 3 - Time history of the tracer concentration to different velocity  field approximations.}
\label{historia-concentracao.tracador-meio.het-cenario03}
\end{figure}

%
% Figura 9
\begin{figure}[!htb]
\centering
\subfigure[]{\includegraphics[height=4.0cm,width=6.1cm]{figures/Semianalitico-Comparacao-PP.eps}} %
\hspace{0.25cm}
\subfigure[]{\includegraphics[height=4.0cm,width=6.1cm]{figures/Semianalitico-Comparacao-SDHM.eps}}
\caption{Time history of the tracer concentration. Comparisons among Abbaszadeh-Dehghani solution with the semi-analytical methodology combined with the (a) Post-Processing technique and the (b) SDHM method.}
\label{historia-concentracao.tracador-comparacao-cenario.hom.vs.cenario3}
\end{figure}

% Figura 10
\begin{figure}[!htb]
\centering
\subfigure[SUPG approach]{\includegraphics[height=4.0cm,width=6.1cm]{figures/TraceVolumeCase3SUPG.eps}}
\hspace{0.25cm}
\subfigure[Semi-analytical methodology]{\includegraphics[height=4.0cm,width=6.1cm]{figures/TraceVolumeCase3Semi.eps}}
\caption{Case 3 - Tracer volume for the tracer concentration. Comparison among the velocity field approximations:  (a) and (b) SDHM, Post-Processing (PP) and Galerkin.}
\label{volume.tracador-meio.het-cenario03}
\end{figure}

\vspace{-6pt}
\section{Conclusion}\label{S6}
\vspace{-2pt}

In this paper the tracer injection process into homogeneous and heterogeneous media was solved using two different methodologies for the concentration approximation combined with three formulations for the velocity field. When the SDHM formulation is employed in the calculation of velocity, combined with the SUPG approach their results lead to more accurate approximations for the concentration. This combination is able to capture the expected discontinuity properties of the solution and, consequently, a proper physical solution. The numerical experiments performed, illustrated the flexibility and robustness of this formulation.



\begin{abstract}
{\bf Resumo}. 
Embora a concentração seja a variável mais importante nos processos de injeção de traçadores, uma eficiente e precisa aproximação do campo de velocidades é crucial para obter um bom comportamento físico para o problema. Neste artigo, analisamos o método misto dual híbrido estabilizado (SDHM) para resolver o sistema de Darcy nas variáveis de velocidade e de pressão a partir da equação de conservação de massa e da lei de Darcy. Esta abordagem é localmente conservativa, livre de comprometimento entre os espaços de aproximação de elementos finitos e capaz de lidar com meios heterogêneos com propriedades descontínuas. A concentração do traçador é resolvida através de uma combinação do método \textit{Streamline Upwind Petrov-Galeklin} (SUPG) no espaço com um método de diferenças finitas implícita no tempo. Também empregamos uma abordagem semi-analítica (solução analítica de Abbaszadeh-Dehghani) para integrar a equação de transporte. Um estudo comparativo numérico utilizando a formulação de SDHM, o método de Galerkin e uma técnica de pós-processamento para calcular o campo de velocidade em combinação com essas metodologias de aproximação da concentração são apresentados. Em todas as comparações, a formulação SDHM aparece como a mais eficiente, precisa e quase sem oscilações espúrias.

{\bf Palavras-Chave}. Deslocamentos miscíveis, Métodos hibridizados, Simulações de reservatórios de óleos

\end{abstract}

\bibliographystyle{ieeetr}
\bibliography{tracer}

\end{document}
