\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small 2014 Madrid Conference on Applied Mathematics in honor of Alfonso Casal,\\ \emph{Electronic Journal of Differential Equations}, Conference 22 (2015), pp. 79--97.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2015 Texas State University.} \vspace{9mm}} \begin{document} \setcounter{page}{79} \title[\hfilneg EJDE-2015/Conf/22 \hfil Two-phase flow models] {Mathematical modeling and numerical simulation of two-phase flow problems at pore scale} \author[P. Luna, A. Hidalgo \hfil EJDE-2015/conf/22 \hfilneg] {Paula Luna, Arturo Hidalgo} \address{Paula Luna \newline Escuela T\'ecnica Superior de Ingenieros de Minas y Energ\'ia, Universidad Polit\'ecnica de Madrid, Madrid, Spain} \email{p.luna@alumnos.upm.es} \address{Arturo Hidalgo \newline Escuela T\'ecnica Superior de Ingenieros de Minas y Energ\'ia, Universidad Polit\'ecnica de Madrid, Madrid, Spain} \email{arturo.hidalgo@upm.es} \thanks{Published November 20, 2015.} \subjclass[2010]{35J20, 35J60, 35J70, 47J10, 46E35} \keywords{Reservoir simulation; numerical simulation; \hfill\break\indent saturation profile; finite volume method} \begin{abstract} Mathematical modeling and numerical simulation of two-phase flow through porous media is a very active field of research, because of its relevancy in a wide range of physical and technological applications. Some outstanding applications concern reservoir simulation and oil and gas recovery, fields in which a great effort is being paid in the development of efficient numerical methods. The mathematical model used in this work is written as a system comprising an elliptic equation for pressure and a hyperbolic one for saturation. Our aim is to obtain the numerical solution of this model by combining finite element and finite volume techniques, with a second-order non-oscillatory reconstruction procedure to build the values of the velocities at the cell interfaces of the FV mesh from pointwise values of the pressure at the FE nodes. The numerical results are compared to those obtained using the commercial code ECLIPSE showing an appropriate behavior from a qualitative point of view. The use of this FE-FV procedure is not the usual numerical method in petroleum reservoir simulation, since the techniques most frequently used are based on finite differences, even in standard commercial tools. \end{abstract} \maketitle \numberwithin{equation}{section} \allowdisplaybreaks \section{Introduction} Multiphase flow in porous media is a very active field of research since this type of problems arise in many practical situations in fluid dynamics, many of them linked to groundwater pollution, oil and gas recovery or $CO_2$ storage in geological formations, to name but a few. Because of the interest of multiphase flow in porous media, many relevant works have been developed to this aim, among which some classical references can be mentioned as \cite{Adler88,bear88,kinzelbach86,Wooding76,vazquez-book07}. The particular context of this work lays on the field of mathematical modeling and numerical simulation of underground flow of oil-water systems. In this context, the fluid in motion is usually formed by a mixture of oil, water and gas moving through the porous media due to an existing connected network of pores and also to the presence of fractures in the geological medium. The most common mathematical models that represent this kind of problems are given by a system of parabolic partial differential equations where the unknowns are the pressure and saturation of the different phases involved in the process, namely the Black-Oil model (see \cite{Tra89}). Many of the simulation tools used to obtain the solution of the mathematical models use finite differences, see \cite{Pea77}, or finite element methods. Other formulations, based on finite volume methods and combinations of finite element and finite volume methods, are also available in the literature \cite{Bergamaschi98}. In this work, we follow a strategy, see \cite{Mon12,Friis12}, where the parabolic system of PDEs is transformed into a system formed by an elliptic equation for the pressure and a hyperbolic equation for the saturation. The latter equation is similar to the well-known Buckley-Leverett equation (first introduced in \cite{Buckley42}). Theoretical analysis of this kind of formulation have been developed in \cite{Coclite12,Scroll2000}. An interesting work in this study corresponds to the treatment of the two phase model based on a conversion into an elliptic-hyperbolic system, developing both analytical study and numerical solution is \cite{Coclite12}. Also in \cite{Coclite12} a numerical resolution of a simplified model is performed. In \cite{Scroll2000} the authors prove local existence and uniqueness of a classical solution for the original hyperbolic-elliptic system arising in the modeling of oil-recovery processes, using the Arzela-Ascoli theorem. Another relevant contribution concerning existence and uniqueness of the solution in models of filtration of immiscible fluids in a porous a media is \cite{Antontsev90}. The elliptic (pressure) equation and the hyperbolic (saturation) equation are coupled by means of the saturation which appears in both of them, more precisely, in the coefficients of mobility. The numerical scheme developed in this work consists of solving the pressure equation (the elliptic part of the model) by means of a finite element method and solving the saturation equation (the hyperbolic part) via a finite volume scheme. The reason behind the use of these two different techniques in the same problem is that finite element methods have been devised specifically to be used in the context of elliptic (and parabolic) problems, but they fail miserably, at least in its classical formulation, when applied to hyperbolic problems, mainly when discontinuous solutions appear. However, finite volume schemes work quite well in these particular situations since they are able to propagate adequately the discontinuities of the solution. In this work we obtain the solution for the pressure at the nodes in which the domain is discretized and, afterwards, these nodes are treated as intercell boundaries for the finite volume scheme. We remark that finite volume scheme is based on integral averages of the solution (saturation in this case). Because we need to know velocities at cell interfaces, and they depend on the gradient of the pressure, according to Darcy's law, a sort of \textit{reconstruction} must be applied. In this work we use a piecewise linear reconstruction, considering the stencil that minimizes the absolute value of the slope. Time integration for the hyperbolic equation is fulfilled via an explicit Euler scheme. We note that these procedure is different to the classical IMPES method, in which both equations are evolution ones and an implicit scheme is used for pressure and an explicit formulation is use to solve the saturation \cite{Pea77,Kle14,Chen04,fanchi06}. We show results of the one-phase case, in which we consider a parabolic equation for the pressure and also results for pressure and saturation are given for the two-phase case. The values of the parameters are realistic and they have been taken from the commercial code ECLIPSE \cite{eclipse} by Schlumberger. The results obtained are compared with those comming from ECLIPSE. In the two-phase case, a qualitative comparison with results obtained from using ECLIPSE is performed. The rest of the paper is organized as follows: Section 2 contains a background on reservoir simulation and on the Black-Oil model, Section 3 deals with the particular mathematical models used in this work, both for the single-phase case and for the two-phase case. Section 4 introduces the numerical scheme used to obtain the solution of the single-phase model which is achieved using a finite element method. Section 5 presents a detailed description of the numerical scheme heading to solve the two-phase model, based on a combination of a finite element scheme for the elliptic part and a finite volume scheme for the hyperbolic part. Section 6 deals with discussion of the results obtained in the previous section. The final part of this work is devoted to conclusions and further research. \section{Physical model} A petroleum reservoir is geological formation in which a mixture of oil, water and gas is trapped. Originally, the hydrocarbon reservoir is in equilibrium containing reservoir fluids such as gas, oil and water separated by gravity at no-flow conditions. When a well is drilled reaching the upper reservoir layer this equilibrium state is immediately disturbed, consequently pressure declines at that specific location for being progressively extended through all radial directions from the well zone to the whole reservoir domain. The study of the movement of the fluid throughout the porous media is very complex. Then, certain simplifying assumptions must be introduced. In this work we establish the following hypothesis: \begin{itemize} \item The flow takes place in one-dimension. \item Fluids are immiscible and their composition is constant over time, therefore there is no mass transfer between phases. \item The fluids are incompressible which means that the density is constant. \item The system is isothermal. \item Capillary pressure is neglected. \item Horizontal flow is assumed, therefore the effect of gravity is neglected. \end{itemize} \subsection*{The Black-oil model} The Black-oil thermodynamic model is one of the most widely used models in petroleum reservoir simulation, see for instance \cite{Tra89}. It assumes the hydrocarbon can be described into two components, namely a heavy component (oil) and a light component (gas), for which its composition remains constant in time along with the assumption that no mass transfer is held between the phases, which means that they are immiscible in each other. At reservoir conditions, both components may appear completely or partially dissolved depending on the reservoir pressure and temperature regime at the reservoir, to form either one or two phases, namely liquid and gas. The Black-oil equations consist of: \begin{itemize} \item Conditions of thermodynamic equilibrium. \item Equation of state which describes a fluid in terms of its fundamental physical properties (constitutive equations). \item Darcy's law for the volumetric flow rates (conservation of momentum). \item Mass conservation equation for each component. \end{itemize} The Black-Oil model describes the reservoir for a constant temperature, what means fluid properties are only function of the pressure response in the reservoir. Therefore it will only be required a table showing this variation in the data respect to a corresponding pressure evolution in time, a PVT table will be defined for each of the phases. The behavior of the dead oil is represented in Figure \ref{fig1}, where it is described oil at pressures above the bubble point for a constant temperature. \begin{figure}[htb] \begin{center} \includegraphics[width=0.7\textwidth]{fig1} % BO_PTdiagram \end{center} \caption{Black-Oil phase diagram} \label{fig1} \end{figure} The bubble point pressure (point 2 in the diagram) is a special case of saturation at which the first bubble of gas is formed. For pressures below the bubble point, there are two hydrocarbon phases present in the reservoir, a liquid phase (gas saturated oil) and gas phase (liberated solution gas). Solution gas starts releasing from oil while crossing the bubble point curve during which its concentration remains constant for the whole process. It is considered dead as dissolved gas does not affect its behavior for a constant temperature remaining always below the critical point. For pressures above the bubble point only one hydrocarbon phase is present in the reservoir, that is liquid oil. This oil is called dead oil or undersaturated oil as it will be subjected to states of pressure always evolving above the bubble point and no gas can be released from the oil. \section{Mathematical models} In this work we deal with a two-phase model in one space dimension useful in the field of reservoir simulation and oil recovery. In this model, oil and water are the only active phases in the reservoir. The aim of this model is to describe the immiscible displacement of oil, which will be moving towards the production well by the action of the water injection into the reservoir. Oil and water are immiscible in each other. The reason for drilling an injector well is to maintain a high pressure in the reservoir by injecting fluids such as water or gas, and to ensure the displacement of hydrocarbons towards the well when the current pressure gradient between the formation and the wellbore is too low to allow the continuity of the production of hydrocarbons. In this section, we establish a mathematical model describing the behavior of fluid flow in a petroleum reservoir, which are mostly a set of nonlinear PDEs. The filtration of fluids in porous materials is represented by a set of mass, momentum, energy conservation and constitutive equations which arranged all together, comprise what is known as flow equations. Isothermal conditions are assumed for simplicity, thus the formulation will not require an additional equation for the energy conservation. The conservation of momentum is determined by the Navier-Stokes equations for fluid mechanics, though simplified for describing a low velocity flow behavior are actually representing the semi-empirical Darcy's equation, which in fact is governing the velocity of flow at the pore scale. Denoting with the subscript $n$ the non-wetting phase, usually oil, and with $w$ the wetting phase, usually water, the addition of both saturations must be the unity, that is \begin{equation} s_w+s_n=1. \end{equation} Another relevant physical magnitude in oil reservoir is the capillary pressure, $p_c$, which depends on water saturation and can be expressed as \begin{equation} p_c(s_w)=p_n-p_w, \end{equation} where $p_n$ is the pressure of the non-wetting phase and $p_w$ is the pressure of the wetting phase. Particularly useful flow in porous media is Darcy's law, which in the case of multiphase flow reads \begin{equation}\label{eq:velocity} v_i=-\Big(\frac{K k_{ri}}{\mu_i}\Big) \frac{\partial P_i}{\partial x}, \quad \text{$i$ is the fluid phase.} \end{equation} In this case gas is present as an active phase in the reservoir, it would always be considered a non-wetting phase. The variable $s$ is related to saturation, $p$ is the phase pressure and $p_c$ stands for the capillary pressure. The relative permeability of phase $i$ is given by $k_{ri}$. $K$ represents the total permeability. In Darcy's law also the total pressure of the fluid, $P_i$ is present. The corresponding PDEs of mass conservation are obtained by the same procedure as single-phase flow equations but without considering a term referring to the rate of accumulation. Otherwise, for obtaining this term, the current volume ($\phi\rho$), where $\phi$ is the mass fraction and $\rho$ is the density, must be multiplied by the phase saturation because of the presence of more than one phase in the system for this case. The general equations for two-phase flow (water and oil respectively) read \begin{equation}\label{eq:2phase} \begin{gathered} \frac{\partial}{\partial x} \Big(\frac{K k_{rw}}{\mu_w B_w}\frac{\partial P_w}{\partial x}\Big)+q_w = \frac{\partial}{\partial t}\Big(\frac{\phi s_w}{B_w}\Big), \\ \frac{\partial}{\partial x} \Big(\frac{K k_{ro}}{\mu_o B_o}\frac{\partial P_o}{\partial x}\Big)+q_o = \frac{\partial}{\partial t}\Big(\frac{\phi s_o}{B_o}\Big), \end{gathered} \end{equation} where $\mu_{i}$ is the viscosity of phase $i$ and $B_i$ is the volumetric factor for phase $i$. The volumetric factor represents the relationship between the volume occupied by the phase $i$ at reservoir conditions and the volume at surface conditions, due to the difference in pressures in both cases. In addition, $q_w$ and $q_o$ are the injection/production rate of water and oil respectively. \subsection*{Single-phase flow model}\label{singlephase} We first consider the mathematical model representing a single phase flowing, whose solution provides the pressure distribution in the reservoir. These one-phase models are extremely useful to identify flow directions, connections between production and injection wells, etc. Also, we use the numerical solution of this model to perform a comparison with results obtained with the code ECLIPSE. In this situation, oil is the only active phase in the reservoir and since we are not considering the presence of gas, oil must be undersaturated (dead oil), what means that the reservoir pressure remains higher than the bubble point pressure. In this case, most of the energy retained in the reservoir is due to the fluid and rock compressibility. The pressure declines as fast as fluids are extracted from a subsaturated reservoir until bubble point is reached. Thereafter, solution gas will become the energy source allowing its displacement. This reservoirs are good candidates for water injection, the purpose is to maintain a high pressure and therefore enhance the oil recovery. When dealing with single-phase flow, nor saturation or relative permeabilities appear in the equations and, therefore, system \eqref{eq:2phase} becomes \begin{equation}\label{SingleP} \frac{\partial}{\partial x} \Big(\frac{K}{\mu B}\frac{\partial P}{\partial x}\Big)+q =\frac{\partial}{\partial t}\big(\frac{\phi}{B}\big). \end{equation} For the special case of an incompressible fluid, porosity of the rock is assumed constant in time, also the incompressibility of the fluid implies its density must be constant in time so the time derivative vanishes obtaining an elliptic equation for the fluid pressure distribution. The general elliptic equation for an incompressible single phase flow in a horizontal domain is \begin{equation} -\frac{\partial}{\partial x}\Big(\frac{K}{\mu B}\frac{\partial P}{\partial x}\Big)=q. \end{equation} \subsubsection*{Radial flow model} In oil reservoir simulation it is very common to work in radial coordinates. The equation describing the flow of an incompressible fluid in radial coordinates because of the presence of a nearby well is \begin{equation}\label{eq:radial} \frac{1}{r} \frac{\partial}{\partial r}\Big(r\frac{\partial P}{\partial r}\Big) =\frac{\phi\mu c}{K}\frac{\partial P}{\partial t} \end{equation} where, $r$ stands for the drainage radius of the well and $c$ represents the compressibility of the total system, both rock and fluid. \subsection*{Two-phase flow model} In this work we transform the parabolic system \eqref{eq:2phase} into an elliptic-hyperbolic system (see \cite{Mon12,Friis12,Coclite12} for more details on this way to proceed). In this approach, the elliptic equation models the distribution of pressure whereas the hyperbolic one is used for the evolution of the saturation. The latter one is similar to the well-known Buckley-Leverett model. An important factor when dealing with multi-phase problems is the mobility of the phase $i$, usually denoted as $\lambda_i$ whose expression is: \begin{equation} \lambda_i=\frac{k_{ri}}{\mu_i}, \end{equation} which represents the ability of a certain phase to \textit{move} with respect to the other phases. In the next subsections the general equations for pressure and saturations are presented. \subsubsection*{General pressure equation} The general pressure equation in this model reads as follows \begin{equation}\label{Pressure} \begin{split} &-\nabla [K\lambda_w\nabla p_w+K\lambda_n\nabla p_n] + c_r\phi\frac{\partial p}{\partial t}\\ &- c_w\big[\nabla p_w K\lambda_w(\nabla p_w)-\phi s_w\frac{\partial p_w}{\partial t} \big] - c_n\big[\nabla p_n K\lambda_w(\nabla p_n)-\phi s_n\frac{\partial p_n}{\partial t} \big]= \tilde{Q}, \end{split} \end{equation} where $\lambda_i$ is the mobility of phase $i$, $c_r$ is the compressibility of the rock, $c_w$ is the compressibility of the wetting phase (usually water), $c_n$ is the compressibility of the non-wetting phase (usually oil) and $\tilde{Q}$ is a specific volumetric injection/production rate term given by $\tilde{Q}=\frac{q_n}{\rho_n}+\frac{q_w}{\rho_w}$. In the particular case of an \emph{immiscible and incompressible flow} ($c_r=c_n=c_w=0$) and one-space dimension, equation \eqref{Pressure} becomes \begin{equation} -\frac{\partial}{\partial x}\Big(K \lambda_w\frac{\partial p_w}{\partial x} + K \lambda_n\frac{\partial p_n}{\partial x}\Big)=\tilde{Q}. \end{equation} We now express the total velocity as function of the global pressure $P$ \begin{equation} v=-K(\lambda_w+\lambda_n)\frac{\partial P}{\partial x} \end{equation} and introduce the term of total mobility $\lambda=\lambda_w+\lambda_n$ to obtain \begin{equation} -\frac{\partial}{\partial x}\Big(K \lambda\big(\frac{\partial P}{\partial x}\big)\Big) =\tilde{Q}, \quad \tilde{Q}=\frac{q_n}{\rho_n}+\frac{q_w}{\rho_w} \end{equation} which can be solved for the global pressure $P$. \subsubsection*{General saturation equation} The analysis of the two-phase flow model is assimilated to the Buckley-Leverett displacement problem for which saturation profiles are calculated when capillary and gravity effects are not taken into account. In the recent work \cite{Adrianov13} there is an interesting study, numerical and analytical, of the Buckley-Leverett equation. Considering that there is an expression relating both saturations; $s_w+s_n=1$, it is only necessary to calculate one of them, being commonly defined for the water saturation $s_w$. \begin{equation} \frac{\partial \phi \rho_w s_w}{\partial t}=\nabla(\rho_w h_w \nabla s_w) -\nabla(\rho_w(f_w v))+q_w \end{equation} where $h_w=-f_w K \lambda_n \frac{\partial P_c}{\partial s_w}$ and the fractional flow function $f_w=\lambda_w/(\lambda_w+\lambda_n)$. For the particular case of an incompressible flow $\phi, \rho_n, \rho_w$ are constant, therefore we can write \begin{equation} \phi\frac{\partial s_w}{\partial t}+\frac{\partial}{\partial t}(f_w(s_w) v) -\frac{\partial}{\partial x}\Big(h_w\frac{\partial s_w}{\partial x}\Big) =\frac{q_w}{\rho_w} \end{equation} In general, the saturation equation has a parabolic nature. Yet, on the pore scale, term $f_w(s)v$ usually dominates the term $-\nabla (h_w\nabla s_w)$, consequently acquiring an hyperbolic character so it needs to be discretized in an alternative way than the pressure equation. \subsubsection*{Equations of the two-phase flow model} An immiscible and incompressible flow is considered for this study, as well as capillary forces are neglected for simplicity. \begin{itemize} \begin{subequations} \label{eq:2phasemodel} \item Pressure equation (diffusive elliptic equation) \begin{equation} -\frac{\partial}{\partial x}\Big(K\lambda_w\frac{\partial P}{\partial x}\Big)=Q_t \label{eq:press-eq} \end{equation} \item Saturation equation (diffusive hyperbolic equation) \begin{equation} \phi \frac{\partial s_w}{\partial t}+\frac{\partial}{\partial x}(f_w(s_w) v) =\frac{q_w}{\rho_w} \label{eq:sat-eq} \end{equation} \end{subequations} \end{itemize} Both equations are non-linearly coupled through the saturation dependent mobility, as well as the pressure dependent velocity. The unknown variables are the global pressure and the saturation of the non-wetting phase (water). The system \eqref{eq:2phasemodel} is equivalent to \eqref{eq:2phase}. Equation \eqref{eq:sat-eq} is quasi-linear since the fractional flow function $f_w$ depends on the saturation $s_w$, via the mobilities $\lambda_w$, $\lambda_n$. However, this formulation is different to the quasi-linear Buckley-Leverett model, where the fractional flow function $f_w^{BL}$ is given by $f_w^{BL}=s_w^2/(s_w^2+\alpha(1-s_w)^2)$ (see for instance \cite{toro-book,leveque02}) where $\alpha\in(0,1)$ to account for the fact that the viscosity of the oil is higher than that of the water. The reason for choosing the model \eqref{eq:sat-eq} is that it is more realistic in the context of reservoir simulation since it takes into consideration the water/oil mobility ratio, which is a key parameter in determining the efficiency of the water/oil displacement process. \section{Numerical approximation of the single-phase flow model} \subsection*{Formulation of the single-phase flow problem} Previously, in Section~\ref{singlephase}, it was established the PDE representing the single-phase flow problem was given in \eqref{SingleP}. After some manipulations, see \cite{Kle14}, it can be written, with the addition of suitable initial and boundary conditions, as \begin{gather*} \frac{\partial}{\partial x}\Big(\frac{K}{\mu B}\frac{\partial P}{\partial x} \Big)=\frac{\phi c_r}{B}\frac{\partial P}{\partial t} \\ \text{I.C.} \quad P(x,0)=P_0(x) \\ \text{B.C. Neumann conditions}\\ (x=0) \quad -\frac{KA}{\mu B}\Big(\frac{\partial P}{\partial x}\Big)_{x=0}=0, \quad\text{no flow boundary}\\ (x=L) \quad -\frac{KA}{\mu B}\Big(\frac{\partial P}{\partial x}\Big)_{x=L}=Q_o(t) , \quad\text{oil production rate}, \end{gather*} where $P_0(x)$ is a prescribed initial pressure. In this approach we consider that the viscosity $\mu_o$ and the volumetric factor $B_o$ are pressure-dependent which confers more realism to the problem. In table~\ref{tab:Table1} values of viscosity and volumetric factor are shown for different values of pressure, which have been taken from ECLIPSE and, in this work, they are used to carry out interpolation for the simulated values of the pressure. Because of this nonlinearity, it is necessary to apply a numerical technique to solve the problem. In this case a finite element technique has been implemented, since this type of numerical methods are very convenient when dealing with parabolic problems. To obtain the weak formulation we multiply the PDE by test functions $\varphi_i(x)$ ($i=1,\dots, n$) and integrate on the domain $[0,L]$ to yield \begin{equation} \int_0^L{\varphi_i(x)\frac{\partial}{\partial x} \Big(\tilde{D}\frac{\partial P}{\partial x}\Big)dx} =\int_0^L{\varphi_j(x)\tilde{a}\frac{\partial P}{\partial t}dx} \end{equation} where we have introduced the coefficients \begin{equation} \tilde{D}:=\frac{K}{\mu_oB_o}, \quad \tilde{a}:=\frac{\phi c_r}{B_o}. \end{equation} Both $\tilde{D}$ and $\tilde{a}$ are functions of pressure, according to table~\ref{tab:Table1}. After integration by parts and introducing the Neumann boundary conditions we obtain \begin{equation} -\int_0^L{\frac{\partial \varphi_i(x)}{\partial x} \Big(\tilde{D}\frac{\partial P}{\partial x}\Big)dx} -Q_o(t)/A=\int_0^L{\varphi_j(x)\tilde{a}\frac{\partial P}{\partial t}dx} \end{equation} We can obtain an approximate solution of the weak formulation by introducing the linear combination $P(x)=\sum_{j=1}^n{\varphi_j(x)P_j(t)}$, where $P_j(t)$ are the nodal values of the pressure which vary with time. Therefore, rearranging terms we have \[ \sum_{j=1}^n{\Big(\int_0^L{\tilde{a}\varphi_j(x)\varphi_i(x)dx}\Big) \frac{dP_j(t)}{dt}} +\sum_{j=1}^n\Big(\int_0^L\tilde{D}\frac{\partial\varphi_j(x)}{\partial x} \frac{\partial\varphi_i(x)}{\partial x}\Big)P_j(t)dx = Q_o(t)/A \] which can be written in a more compact way as \begin{equation} C_{ij}\frac{dP_j(t)}{dt}+R_{ij}P_j(t)=V_i \quad (i,j=1,\dots,N) \end{equation} where, \[ C_{ij}=\int_0^L\tilde{a}\varphi_i(x)\varphi_j(x)dx,\quad R_{ij}=\int_0^L\tilde{D}\frac{\partial\varphi_j(x)}{\partial x} \frac{\partial\varphi_i(x)}{\partial x}dx \quad V_i=Q_o(t)/A\,. \] \begin{table}[htb] \caption{Data of viscosity and volumetric factor for different values of pressure. Data taken from ECLIPSE code.} \label{tab:Table1} \begin{center} \begin{tabular}{|c|c|c|} \hline $P$ (bar) & $B_o $ & $\mu_o$ (cp) \\ \hline 180.0 & 1.2601 & 1.041 \\ \hline 227.0 & 1.2600 & 1.042 \\ \hline 253.4 & 1.2555 & 1.072 \\ \hline 281.6 & 1.2507 & 1.096 \\ \hline 311.1 & 1.2463 & 1.118 \\ \hline 343.8 & 1.24173 & 1.151 \\ \hline 373.5 & 1.2377 & 1.174 \\ \hline 395.5 & 1.2356 & 1.2 \\ \hline \end{tabular} \end{center} \end{table} This system of equations must be solved by means of an appropriate ODE solver. In this work we have used a \emph{Backward Euler scheme} which is unconditionally stable, therefore we have the scheme \begin{equation} (C_{ij}+\Delta tR_{ij})P_j^{n+1}=C_{ij}P_j^n+\Delta t V_i(t^{n+1}) \end{equation} Compacting this expression, we obtain the approximate solution of the implicit pressure, computed for each time step by solving the system of equations \begin{equation} A_{ij}P_j^{n+1}=B_i^n, \end{equation} where \[ A_{ij} =C_{ij}+\Delta tR_{ij}, \quad B_{i}^n =C_{ij}P_j^n+\Delta tV_i(t^{n+1}). \] In this work we have used Gaussian elimination adapted to a tridiagonal system. Because of the dependency of viscosity and volumetric factor with pressure, we perform linear interpolation, based on table~\ref{tab:Table1}, taking into account the values of the pressure of the previous time step. \subsection*{Formulation of the radial flow model} We consider an infinite reservoir (constant pressure at the boundaries) at an initial pressure state. \begin{gather*} \frac{1}{r} \frac{\partial}{\partial r} \Big(r\frac{\partial P}{\partial r}\Big) =\frac{\phi\mu c}{K}\frac{\partial P}{\partial t} \\ \text{I. C.\quad $P=172.35$ bar}\\ \text{B. C. \quad Dirichlet conditions}\\ (r=0) \quad P=124.1\text{ bar}\\ (r=L) \quad P=172.35\text{ bar}. \end{gather*} This problem provides a solution to the radial flow example appearing in \cite{Ahm10} page 347. Same boundary conditions and initial pressure values have been applied to this model in order to compare the solution obtained with the one corresponding to the example. It is necessary to point out that this model has been solved according to different values for certain parameters, that is viscosity of oil $\mu_o$, which in this case does not correspond to a constant value but as function of pressure, resulting a more realistic approach (we are solving for real values). The field is producing at a variable rate (in contrast with the constant production rate of the example in \cite{Ahm10}) at a stabilized bottom-hole flowing pressure of 124.1 bar. This problem is solved by the Finite Element method for obtaining an approximate solution of the parabolic equation \eqref{eq:radial}. For the weak formulation, solve $P(r,t)\in H^1(0,L) / P(0,t)=P_0(t)$ and $P(r,0)=P_0(r)$, solution of \begin{equation} \int_0^L \dfrac{\partial w}{\partial r} \Big(r\frac{\partial P}{\partial r}\Big)dr =\int_0^L w \: r\frac{\phi\mu c}{K}\dfrac{\partial P}{\partial r}dr \quad \forall w \in V_h. \end{equation} The approximate solution is expressed as follows: $P_h(r,t) \in V_h / P_h(0,t)=P_0(t)$ and $P_h(r,0)=P_0(r)$, solution of \begin{equation} \int_0^L \dfrac{\partial w}{\partial r} \Big(r\frac{\partial P_h}{\partial r}\Big)dr =\int_0^L w \: r\frac{\phi\mu c}{K}\dfrac{\partial P_h}{\partial r}dr, \end{equation} where \begin{equation} P_h(r,t)=\sum_{j=1}^N P_j(t)\varphi_j(r) \label{eq:P_h(x,t)} \end{equation} where $P_j$ stands for the approximate pressure solution at the nodes and $\varphi_j$ representing the basis functions (piecewise polynomial functions). Test functions can be chosen as \begin{equation} w(r)=\varphi_i(r), \quad (i=1,\dots,N) \label{eq:w(x)} \end{equation} Using \eqref{eq:P_h(x,t)} and \eqref{eq:w(x)} in the expression of the approximate solution for the pressure, we obtain \begin{equation} \sum_{j=1}^N \Big(\int_0^L r\varphi'_i \varphi'_j dr\Big) P_j(t) =\sum_{j=1}^N \Big( \int_0^L r \frac{\phi\mu c}{K} \varphi_i \varphi_j dr \Big) \frac{dP_j(t)}{dt} \end{equation} which in fact can be expressed as \begin{equation} R_{ij}P_j(t)=C_{ij}\frac{dP_j(t)}{dt}\quad (i,j=1,\dots,N), \end{equation} where \begin{equation} R_{ij}=\int_0^L r\varphi'_i \varphi'_j dr, \quad C_{ij}=\int_0^L r \frac{\phi\mu c}{K} \varphi_i \varphi_j dr\,. \end{equation} This system of equations is solved by the \emph{Backward Euler scheme} for time discretization. \begin{equation} (C_{ij}+\Delta t R_{ij})P_j^{n+1}=C_{ij}P_j^n \end{equation} Compacting this equation, we obtain the approximate solution of the implicit pressure, computed for each time step by solving the following system of equations, \begin{equation} A_{ij}P_j^{n+1}=B_i, \end{equation} where \[ A_{ij} = C_{ij}+\Delta t R_{ij} ,\quad B_i = C_{ij}P_j^n \] \section{Numerical approximation of the two-phase flow problem} We consider the two-phase problem \begin{gather*} -\frac{\partial}{\partial x}(K \lambda(s) \frac{\partial P}{\partial x})=\hat{Q} \\ \phi\frac{\partial s}{\partial t}+\frac{\partial}{\partial x}(f_w(s)v) =\frac{q_w}{\rho_w} \\ \text{I. C.}\quad s_w=\begin{cases} 0.044624 & (x=0.04\text{ km}) \\ 0 & (0.04 \text{ km} 0 \,,\\ \overline{s}_{i+1} & \text{if } v(x_i,t^n) < 0 \,, \end{cases} \end{equation} where positive velocity means that flow goes right. Since the scheme is explicit a \eqref{eq:esqnumsat} a restrictive $cfl$ condition to keep stability of the numerical solution must be applied. Another relevant issue is that saturation depends on the velocity, which is obtained from the pressure equation, as described in the next subsection. \subsection*{Coupling pressure-saturation} In each time step the pressure equation and saturation equation are solved sequentially, according to the following strategy: \smallskip \noindent\textbf{First stage.} Starting with an initial profile of the saturation: $s(x,0)$ the pressure equation is solved using the FEM considering the intervals $[x_i,x_{i+1}]$ as finite elements. After applying the numerical scheme the nodal values of the pressure, $P_i^{n+1}$, are computed. \smallskip \noindent\textbf{Second stage.} Values of the velocity at each node $x_i$ are calculated, according to Darcy's law \eqref{eq:velocity}. The gradients of the pressure are computed using a \textit{reconstruction} procedure based, in this work, on obtaining a first degree polynomial, $R(x,t^n)=P_i^n+m_i(x-x_i)$, for each nodal point whose derivative is used to approximate the gradient of the pressure, \begin{equation}\label{eq:gradP} \frac{\partial P}{\partial x}(x_i,t^n)\approx m_i \end{equation} This is a second order reconstruction which in principle is non-monotone thus, to avoid oscillations, we resort to a TVD scheme based on choosing the slope according to: $|m_i|=min({|m_{i,-1}|,|m_{i,0}|,|m_{i,1}|})$ where $m_{i,-1}=(P_i^{n+1}-P_{i-1}^{n+1})/\Delta x$, $m_{i,0}=(P_{i+1}^{n+1}-P_{i-1}^{n+1})/(2\Delta x)$ and $m_{i,1}=(P_{i+1}^{n+1}-P_i^{n+1})/\Delta x$ where we have considered an equally spaced mesh of size $\Delta x$ for simplicity. \smallskip \noindent\textbf{Third stage.} Once the nodal values of the velocity are computed by using \eqref{eq:gradP} in \eqref{eq:velocity} we compute the interface fluxes by introducing the interface velocities in $f_w(s_{i+1}^n)v_{i+1}$. \smallskip \noindent\textbf{Fourth stage.} The interface fluxes calculated in the previous stage are plugged into \eqref{eq:esqnumsat} to get the final numerical scheme for the saturation. \smallskip After applying the corresponding discretization techniques, a time resolution method is required for solving the problem. In this study we used a sequential IMPES-type method (implicit pressure - explicit saturation). \section{Results and discussion} We have solved the mathematical models describing the single-phase flow problem, both in cartesian and radial coordinates which are solved using a finite element method as described in previous sections. The computed codes are written in Fortran90 and they will be referred to in the following as SPFLOW and SPFLOW-R respectively. We have also solved the two-phase flow mathematical model with the combined finite element-finite volume formulation as described in previous sections, generating a Fortran code which will be referred to in the following as TPFLOW. \subsection*{Single-phase model} The results obtained are depicted in figure \ref{fig3} where one can notice the evident similarities revealed in the declining tendency of the curves representing the results with SPFLOW and ECLIPSE, also pointing out the decay in pressure always remains above the bubble point, which is set at 180 bar. \begin{figure}[htb] \begin{center} \includegraphics[width=0.6\textwidth]{fig3} %FPR2 \end{center} \caption{Average Reservoir Pressure distribution obtained from ECLIPSE and the SPFLOW code} \label{fig3} \end{figure} Results show the reservoir pressure evolves from an initial value of 393.6 bar (I. C. of the problem) to a final state reaching 294.49 bar; for the SPFLOW code. ECLIPSE is initialized with a pressure of 393.6 bars, ending up with a pressure of 296.73 bar after 540 days. It is also essential to point out that the data handled for ECLIPSE and so for the SPFLOW, SPFLOW-R and TPFLOW codes is real data. In this way, as an example, ECLIPSE requires a table alluding to the properties of the fluids contained, called PVT table, which shows the dependency of the formation volume factor, compressibility and viscosity of the fluids with pressure, for a known temperature. For the SPFLOW code, we have used those tables by applying several interpolations of those parameters with its corresponding pressures, so as to match the higher the possible the behavior of the problem analyzed by both tools. \subsubsection*{Radial flow model} \begin{figure} \begin{center} \includegraphics[width=0.5\textwidth]{fig4} % radial1 \end{center} \caption{Pressure profile in radial coordinates for the single-phase flow model (SPFLOW-R). Output time 540 days.} \label{fig:radial} \end{figure} \begin{figure} \begin{center} \includegraphics[width=0.5\textwidth]{fig5} % wellbore_area1 \end{center} \caption{Pressure profile in radial coordinates for the single-phase flow model (SPFLOW-R) in the wellbore area. Output time 540 days.} \label{fig:wellbore} \end{figure} In figure \ref{fig:radial} it is represented the pressure distribution (pressure profile) over the drainage radius. Permeability value used is 120 mD, matching the examples. As indicated in the previous section, viscosity data does not match the corresponding value in \cite{Ahm10}. The applied range of variation for viscosity $\mu_o$ goes from 1.041 cp to 1.2 cp, instead of the 2.5 cp constant value used in the example. The drainage radius is 227 m, matching the data in the example. By performing a qualitative analysis it is noticeable that, in terms of behavior, the resulting pressure profile is almost matching the pressure distribution of the example, setting aside the differences arising due to the use of different input data. To see the resulting pressure values compared with those appearing in the example, we made a quantitative analysis for which we have displayed a table showing the pressures obtained by the SPFLOW-R in contrast with the pressures resulting in example from \cite{Ahm10} for a corresponding radius. \begin{table}[htb] \caption{Pressure 1 represents the pressure from the example analyzed in \cite{Ahm10} for a given radius. Pressure 2 represents the corresponding pressure values obtained from the SPFLOW-R code for a given radius.} \label{table:Pcomparison} \begin{center} \begin{tabular}{|l|l|c|c|} \hline Radius (ft) & Radius (m) & Pressure 1 (bar) & Pressure 2 (bar) \\ \hline 0.25 & 0.0762 & 124.11 & 124.10 \\ 1.25 & 0.381 & 133.90 & 132.46 \\ 4 & 1.2192 & 141.00 & 139.57 \\ 5 & 1.524 & 142.31 & 140.98 \\ 19 & 5.7912 & 150.44 & 149.36 \\ 20 & 6.096 & 150.72 & 149.68 \\ 99 & 30.1752 & 160.51 & 159.70 \\ 100 & 30.48 & 160.58 & 159.76 \\ 744 & 226.7712 & 172.79 & 172.34 \\ 745 & 227.076 & 172.80 & 172.35 \\ \hline \end{tabular} \end{center} \end{table} \begin{figure}[htb] \begin{center} \includegraphics[width=0.7\textwidth]{fig6} % Pcomparison \end{center} \caption{Logarithmic trend line to fit the data from table~\ref{table:Pcomparison}. In blue circles results from \cite{Ahm10}, in red squares results obtained with SPFLOW-R. Color version on-line} \label{fig:Pcomparison} \end{figure} \begin{table}[htb] \caption{For a given radius interval, there are represented both, the pressure drop 1 corresponding to the values of the example from\cite{Ahm10} and pressure drop 2 representing the corresponding values obtained from the SPFLOW-R code.} \label{table:Pdropcomparison} \begin{center} \begin{tabular}{|c|c|c|c|} \hline Radius (ft) & Radius (m) & Press. drop 1 (bar) & Press. drop 2 (bar) \\ \hline 0.25 - 1.25 & 0.0762 - 0.381 & 9.790555352 & 8.357534417 \\ 4 - 5 & 1.2192 - 1.524 & 1.310003885 & 1.407243027 \\ 19 - 20 & 5.7912 - 6.096 & 0.275790292 & 0.321533055 \\ 99 - 100 & 30.1752 - 30.48 & 0.068947573 & 0.062983543 \\ 744 - 745 & 226.7712 - 227.076 & 0.006894757 & 0.008417394 \\ \hline \end{tabular} \end{center} \end{table} \subsection*{Two-Phase model} \begin{figure}[htb] \begin{center} \includegraphics[width=0.6\textwidth]{fig7} % swfortrancode \end{center} \caption{Water saturation profile for 3600 days of performance; TPFLOW code} \label{fig:swfortran} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.6\textwidth]{fig8} % sw \end{center} \caption{Water saturation profile after 3600 days of performance; ECLIPSE code} \label{fig:sweclipse} \end{figure} In figure \ref{fig:swfortran} it is represented the final water saturation profile for the domain after 3600 days of production, obtained by the TPFLOW code. It is actually displaying the saturation profile for each time step until a final saturation profile is reached, corresponding to the outermost part of the colored area. The performance of the saturation profile over the domain after 3600 days by the TPFLOW code shows less diffusion. Our goal was not to analyze the behavior of the problem quantitatively but qualitatively, what actually reproduces the expected behavior of diffusive problems. ECLIPSE: the saturation front moves up to 0.8 km of the domain extension with a water saturation value of 0.005598. Capillary pressure has not been considered for this study. On the other hand both cases have been run, that is, with and without considering the capillary pressure effect in ECLIPSE, in order to compare results so to be aware of the influence the capillary pressure has in the distribution of the saturation profiles. \section{Conclusions and further work} In the single phase flow problem, pressure evolution shows a parallel tendency in both results from ECLIPSE and the solution obtained from the SPFLOW routine. In the study of the two-phase flow, the treatment of the problem has proven effective and appropriate, as in terms of the saturation solution, results show evident similarities with respect to the solution obtained by ECLIPSE. The advance of the saturation front is regulated by the fractional flow rates derivatives to water saturation. By a qualitative analysis it is shown a constant progress of the saturation front. According to future considerations, it will be interesting to make a study in detail about the effect the capillary pressure has in the global pressure, which in fact has not been considered for this study. In this work porosity $\phi$ and permeability $K$ have been assumed constant, though in real cases this parameters are discontinuous functions. Production flow rates are not constant over time because of the fact that reservoirs contain limited volume of reserves, so a decline in the flow rates will take place during the final stages of production (because of the depletion of reserves), that is why it has been considered a transient source term in the formulation, which is actually an output obtained after running the model in ECLIPSE. The fractional flow formulation has proven to be a good choice in terms of performance and approach to the solution obtained by ECLIPSE. For this case (two-phase flow); the Buckley-Leverett model, has been solved for a variable flow rate instead of a constant rate which in fact is what this formulation is established to work with. Which means it is not correct to mention we are solving for the Buckley-Leverett model strictly. \begin{thebibliography}{00} \bibitem{Adler88} P. M. Adler, H.~Brenner; \emph{Multiphase flow in porous media}, Annual Review of Fluid Mechanics \textbf{20} (1988), no.~1, 35--59. \bibitem{Adrianov13} B.~Adrianov, C.~Canc\`es; \emph{Vanishing capillarity solutions of buckley-leverett equation with gravity in two-rocks medium}, Computational Geosciences \textbf{17(3)} (2013), no.~4, 551--572 (eng). \bibitem{Ahm10} T.~Ahmed; \emph{Reservoir engineering handbook, fourth edition}, Gulf Professional Publishing, 2010. \bibitem{Antontsev90} S. N.~Antontsev A. V.~Kazhikhov, V. N. Monakhov; \emph{Boundary value problems in mechanics of nonhomogeneous fluids}, Studies in Mathematics and Its Applications, vol.~22, North-Holland, 1990. \bibitem{bear88} J.~Bear; \emph{Dynamics of fluids in porous media}, Dover Civil and Mechanical Engineering Series, Dover, 1988. \bibitem{Buckley42} S. E. Buckley, M. C. Leverett; \emph{Mechanism of fluid displacements in sands.}, Transactions of the AIME \textbf{146} (1942), 107--116. \bibitem{fanchi06} J. R. Fanchi; \emph{Principles of applied reservoir simulation}, Gulf Professional Pub., 2006. \bibitem{Friis12} H.~Friis and S.~Evje; \emph{Numerical treatment of two-phase flow in capillary heterogeneous porous media by finite-volume approximations}, Int. J. Numer. Anal. Mod. \textbf{9} (2012), no.~3, 505--528. \bibitem{Coclite12} S.~Mishra G.~Coclite, K.~Karlsen, N.~Risebro; \emph{A hyperbolic-elliptic model of two-phase flow in porous media-existence of entropy solutions.}, Int. J. Numer. Anal. Mod. \textbf{9} (2012), no.~3, 562--583. \bibitem{kinzelbach86} W.~Kinzelbach; \emph{Groundwater modelling: An introduction with sample programs in basic}, Developments in Water Science, Elsevier Science, 1986. \bibitem{Kle14} J.~Kleppe; \emph{Reservoir simulation lecture notes}, 2014. \bibitem{Bergamaschi98} S.~Mantica L.~Bergamaschi and G.~Manzini; \emph{A mixed finite element--finite volume formulation of the black-oil model}, SIAM Journal on Scientific Computing \textbf{20} (1998), no.~3, 970--997. \bibitem{leveque02} R. J. LeVeque; \emph{Finite volume methods for hyperbolic problems}, Cambridge Texts in Applied Mathematics, Cambridge University Press, 2002. \bibitem{Mon12} F.~Monkeberg; \emph{Finite volume methods for fluid flow in porous media}, (2012). \bibitem{Pea77} D. W. Peaceman; \emph{Fundamentals of numerical reservoir simulation}, Elsevier, 1977. \bibitem{eclipse} Schlumberger; \emph{Eclipse reference manual 2009.1}. \bibitem{Scroll2000} A.~Schroll, A.~Tveito; \emph{Local existence and stability for a hyperbolic-elliptic system modeling two-phase reservoir flow}, Electron. J. Differential Equations \textbf{2000} (2000), no.~4, 1--28 (eng). \bibitem{toro-book} E. F. Toro; \emph{{Riemann} solvers and numerical methods for fluid dynamics}, third ed., Springer, 2009. \bibitem{toro09} E. F. Toro; \emph{Riemann solvers and numerical methods for fluid dynamics: A practical introduction}, Springer, 2009. \bibitem{Toro-Hidalgo09} E. F. Toro, A.~Hidalgo; \emph{Ader finite volume schemes for nonlinear reaction-diffusion equations}, Applied Numerical Mathematics \textbf{59} (2009), no.~1, 73 -- 100. \bibitem{Tra89} J.~Trangenstein, J.~Bell; \emph{Mathematical structure of the black-oil model for petroleum reservoir simulation}, SIAM Journal on Applied Mathematics \textbf{49} (1989), no.~3, 749--783. \bibitem{vazquez-book07} J. L. V{\'a}zquez; \emph{The porous medium equations. mathematical theory}, Oxford Mathematical Monographs, Oxford Science Publications, 2007. \bibitem{Wooding76} R. A. Wooding, H. J. Morel-Seytoux; \emph{Multiphase fluid flow through porous media}, Annual Review of Fluid Mechanics \textbf{8} (1976), no.~1, 233--274. \bibitem{Chen04} G.~Huan Z.~Chen, B.~Li; \emph{An improved impes method for two-phase flow in porous media}, Transport in Porous Media \textbf{54} (2004), no.~3, 361--376. \end{thebibliography} \end{document}