Semi-Implicit Finite Volume Procedure for Compositional Subsurface Flow Simulation in Highly Anisotropic Porous Media

Abstract

Subsurface multiphase flow in porous media simulation is extensively used in many disciplines. Large meshes with non-orthogonalities (e.g., corner point geometries) and full tensor highly anisotropy ratios are usually required for subsurface flow applications. Nonetheless, simulations using two-point flux approximations (TPFA) fail to accurately calculate fluxes in these types of meshes. Several simulators account for non-orthogonal meshes, but their discretization method is usually non-conservative. In this work, we propose a semi-implicit procedure for general compositional flow simulation in highly anisotropic porous media as an extension of TPFA. This procedure accounts for non-orthogonalities by adding corrections to residual in the Newton-Raphson method. Our semi-implicit formulation poses the guideline for FlowTraM (Flow and Transport Modeller) implementation for research and industry subsurface purposes. We validated FlowTraM with a non-orthogonal variation of the Third SPE Comparative Solution Project case. Our model is used to successfully simulating a real Colombian oil field.

1. Introduction

Subsurface flow in porous media remains an important area in many biological, industrial and scientific applications. Numerical simulations of coupled multiphase-multicomponent flow in porous media, using continuum mechanics approach [1], have been extensively applied in various fields including, hydrology [2,3], soil remediation [4,5], ${\mathrm{CO}}_2$ storage [6], oil and gas, composite materials manufacturing [7,8]. One important application is enhanced oil recovery (EOR) in which various multicomponent phenomena arise, e.g., foam transport and generation [9,10], chemical flooding [11,12], in-situ upgrading [13], and asphaltene precipitation [14,15], flocculation, and remediation [16,17]. Due to this high number of applications, and the continuous development of new techniques for EOR based on chemical and thermal processes, robust and flexible numerical simulation frameworks are required to evaluate and optimize the field deployment of such technologies.

The solution of the partial differential equations (PDE) originated from flow and transport models in porous media, is approximated using conservative numerical methods such as finite volume method (FVM) [18,19]. Subsurface spacial domains have complex geometries that are usually modeled by unstructured or non-orthogonal grids (e.g., PEBI grids [20] or corner-point geometries [21,22]). The advective term in FVM is often discretized following the two-point flux approximation (TPFA) [23]. The TPFA is the industry-standard methodology, in which the local K-orthogonality effect must be considered to achieve numerical convergence [24]. Numerical meshes in many cases, have highly anisotropic ratios inducing an additional source of error [25,26]. Consequently, several authors proposed methods for computing numerical gradients at element surfaces to enable accurate solutions with TPFA in non-orthogonal grids [27,28,29,30,31,32], but most of these methods were developed primarily for isotropic tensorial field applications. In most cases, those are limited to 2-dimensional geometries or simple single-phase fluid problems for anisotropic porous media applications [33]. Thus, the numerical method formulation to be developed, must be general, flexible, and easily adaptable to admit the various physical phenomena involved in many environmental, industrial and biological applications in porous materials.

There are currently several simulators in development for multiphase transport in porous media, mainly to subsurface applications. On the industry side, the leading commercial simulators are ECLIPSE [34], GEM [35], and STARS [36], which are used in oil and gas applications. The first uses a finite difference method (FDM); the second and third use FVM with finite differences to approximate the numerical fluxes at the faces. On the other hand, there are research purpose simulators like Finite Element Heat and Mass Transfer Simulator (FEHM) [37] and OpenGeoSys [38]. These are based primarily on the finite element method (FEM), which is also locally non-conservative for any fluid phases, because of the imposed continuity of the solution through element faces [39]. Another more elaborated academic simulators such as Stanford General Purpose Research Simulator (GPRS) [40,41] uses FVM with a multi-point flux approximation (MPFA) technique instead of TPFA to deal with non-orthogonality and full tensor effects in arbitrary polyhedral grids. There are also parallel oriented simulation platforms like Integrated Parallel Accurate Reservoir Simulators (IPARS) [42], which is a framework for developing parallel models of subsurface flow and transport through porous media; and DuMu${}^{\mathrm{x}}$ [43], a compositional simulator built on top of DUNE [44,45], a modular toolbox for solving partial differential equations. Besides, it is worthy of mention Open Porous Media simulator (OPM) [46] and the MATLAB Reservoir Simulation Toolkit MRST [47], designed for rapid prototyping and evaluation for reservoir simulation problems [48]. The last one uses a Mimetic finite difference method (MFDM) to achieve very accurate solutions in unstructured geometries at the cost of higher computational complexity. Moreover, the methodology developed in this paper is the theoretical foundation of FlowTraM (Flow and Transport Modeller), a compositional subsurface flow simulator developed at the National University of Colombia for research purposes in developing and testing new EOR technologies.

In this work, we present an unstructured finite volume method (UFVM), with an improved TPFA, which allows for managing mesh non-orthogonality with full tensor permeability effect. Mesh skewness [49] is taken into account by a direct consequence of our methodology formulation, which is a generalization to the 3-dimensional case of the previous work of Loudyi [32], applied to non-linear multiphase flow. Besides, the loss of convergence rate, due to the K-orthogonality effect, is compensated by introducing additional corrective fluxes in Newton’s method formulation. These corrections, non-orthogonality, and skewness, require gradient reconstruction techniques. This gradient reconstruction is done efficiently by constructing preprocessed stencils in the weighted least-squares sense (WLS). We use a generalized Crank Nicholson scheme [50] (also known as implicit-explicit or IMEX method [51,52]) for time discretization. Furthermore, our semi-implicit procedure serves to approximate the corrective fluxes, taking advantage of the IMEX method formulation, the current time solution, and the internal iterations of Newton’s method for the resulting non-linear algebraic equations system. The discretization presented here for the advective transport term employs the second-order anisotropic elliptic operator [53,54], which enables a trivial extension of the proposed methodology to any physical transport phenomena originated by a gradient law [55] in anisotropic tensorial fields. Also, our methodology conception allows for naturally extending the functionality of already implemented flow simulators based on FVM for cartesian meshes, increasing its functionality to non-orthogonal structured meshes, as is the case for corner point geometries in a natural way.

This paper is organized as follows: First, we present the mathematical model, emphasizing the physical principles of transport phenomena, using a continuum mechanics approach. In this context, the conservation equations for an arbitrary number of components are presented. Subsequently, we show our UFVM formulation, including some remarks about our implementation (FlowTraM). Furthermore, we present three numerical examples, validating our method and verifying the effectiveness of corrective fluxes. Finally, we demonstrate the robustness and flexibility of this numerical approach simulating an actual Colombian reservoir under EOR operation.

2. Mathematical Model

2.1. Transport Models

The total flux of transport of species k in phase p is expressed in terms of pressure, molecular diffusion, mechanical dispersion, and other factors as [55]:

$${\mathbf{J}}_{k,p}=\underset{\mathrm{Advection}}{\underbrace{{\mathbf{J}}_{k,p}^A}}+\underset{\mathrm{Hydrodynamic}\mathrm{dispersion}}{\underbrace{{\mathbf{J}}_{k,p}^M+{\displaystyle {\mathbf{J}}_{k,p}^D}}}+\underset{\text{Other factors}}{\underbrace{{\displaystyle \sum _O}{\mathbf{J}}_{k,p}^O}},$$

(1)

The transport of components in a fluid medium can occur due to the effect of various internal forces [56], which, in general, are a consequence of the presence of gradients, either from pressure, concentration, or temperature [57]. Advection and hydrodynamic dispersion are the main mechanisms dominating transport in porous media [55]. However, other types of effects may be relevant depending on the case study, such as the Soret effect induced by temperature gradients [58,59] or the chemotaxis phenomenon, in which certain types of microorganisms move through the porous medium due to the presence of nutrient concentration gradients [60,61,62]. However, although this work is limited to the modeling of advection transport, the discretization presented here can be extended to any phenomenon expressed by a gradient law. Thus, the transport of species by advection is given by:

$${\mathbf{J}}_{k,p}^A=w_{k,p}{\rho }_p{\mathbf{u}}_p,$$

(2)

where $w_{k,p}$ is the molar fraction of the component k in the phase p, and ${\mathbf{u}}_p$ the velocity of each phase given by the Darcy equation for multiphase flow, which establishes a relationship of proportionality with the phase potential ${\Phi }_p(\mathbf{x},t)$ as:

$${\mathbf{u}}_p=-\frac{\mathbb{K}{k}_{r,p}}{{\mu }_p}\nabla {\Phi }_p,$$

(3)

and $\mathbb{K}(\mathbf{x},t)$ is the symmetric positive definite tensor of absolute permeability, $k_{r,p}\left(S_p\right)$ the relative permeability of the phase p and ${\mu }_p(\mathbf{x},t)$ the phase viscosity. The ratio ${\lambda }_p:=k_{r,p}/{\mu }_p$ is known as the mobility of the phase p. Additionally, the potential is defined in terms of pressure as:

$${\Phi }_p(\mathbf{x},t)={\int }_{P_p^0}^{P_p}\frac{d{P}_p}{{\rho }_p}+g(z-z_0),$$

(4)

2.2. Compositional System of Equations

Let $\Omega \subset {\mathbb{R}}^3$ an open bounded set and $\Theta :=(t_0,t_f)\subset \mathbb{R}$. The system of continuity equations, for $n_k$ number of components, distributed in $n_p$ phases is

$$\frac{\partial }{\partial t}{N}_k(\mathbf{x},t)+\mathrm{div}\left({\mathbf{F}}_k(\mathbf{x},t)\right)+{\dot{q}}_k(\mathbf{x},t)=0,\forall k\in \mathcal{K},(\mathbf{x},t)\in \Omega \times \Theta ,$$

(5)

where $\mathcal{K}:=\left\{1,2,\dots ,n_k\right\}$ is the family of indices for each component, ${\mathbf{F}}_k={\sum }_p{\mathbf{J}}_{k,p}^A$ is the molar flux by advection, and $N_k$ is the number of moles per unit of volume. Additionally, let ${\Gamma }_{\mathrm{N}},{\Gamma }_{\mathrm{D}}\subset \partial \Omega$ be open sets such that ${\Gamma }_{\mathrm{N}}\cap {\Gamma }_{\mathrm{D}}=\varnothing$ and $\overline{{\Gamma }_{\mathrm{N}}}\cup \overline{{\Gamma }_{\mathrm{D}}}=\partial \Omega$, the Neumann, and Dirichlet type boundaries respectively. Equation (5) satisfies the following boundary conditions

$$\left\{\begin{array}{c}\begin{array}{ccccc}\hfill -{\displaystyle \frac{k_{r,p}}{{\mu }_p}}〈\mathbb{K}\nabla P_p,\mathbf{n}〉& =& g_{p,\mathrm{N}}^0,& \mathrm{on}& {\Gamma }_{\mathrm{N}},\\ \hfill P_p& =& g_{p,\mathrm{D}}^0,& \mathrm{on}& {\Gamma }_{\mathrm{D}},\\ \hfill S_p& =& g_{s,\mathrm{D}}^0,& \mathrm{on}& {\Gamma }_{\mathrm{D}}.\end{array}\hfill \end{array}\right.$$

(6)

Usually, to indicate that the domain is physically closed, in a portion of the boundary, $g_{p,\mathrm{N}}^0=0$ is set, indicating that $\nabla {\Phi }_p=\mathbf{0}$ or that the permeability at the border is zero i.e., $\mathbb{K}={\mathbf{0}}_{3\times 3}$. Since Equation (5), subject to boundary conditions (6), models the evolution of the system in time interval $\Theta$, an initial configuration of the variables is required, which is

$$\left\{\begin{array}{c}\begin{array}{cccc}{P}_p(\mathbf{x},t_0)& =& P_p^0\left(\mathbf{x}\right),& \forall p\in \mathcal{P},\\ S_p(\mathbf{x},t_0)& =& S_p^0\left(\mathbf{x}\right),& \forall p\in \mathcal{P},\\ w_{k,p}(\mathbf{x},t_0)& =& w_{k,p}^0\left(\mathbf{x}\right),& \forall (p,k)\in \mathcal{P}\times \mathcal{K},\end{array}\hfill \end{array}\right.$$

(7)

where $\mathcal{P}$ is the family of phase indexes. The functions $P_p^0\left(\mathbf{x}\right)$, $S_p^0\left(\mathbf{x}\right)$, and $w_{k,p}^0$ must satisfy the compatibility equations shown in the next section.

The main unknowns of the problem are the pressure functions for each of the phases $P_p$, the saturations $S_p$, and the molar fractions $w_{k,p}$ of the components distributed in the phases. Including this, the total number of unknowns is $n_{\mathrm{IN}}:=2n_p+n_p{n}_k$, of which $2n_p$ term corresponds to the phase, pressure and saturation unknowns, and $n_p{n}_k$ is the number of unknowns corresponding to the molar fractions.

2.3. Compatibility Equations

Properly speaking, the gradient laws for transport phenomena are compatibility equations, because they satisfy the material objectivity axioms [63]. Nevertheless, they are assumed here as part of the governing equations, and the constitutive relations or compatibility equations are stated as follows:

2.3.1. Thermodynamic Equilibrium

Transport equation terms and particularly their non-linearities are highly related to the thermodynamic state and properties of the existing phases in a specific space-time location. Nowadays, most of the models working with local thermodynamic equilibrium assumption are capable of reproducing a wide range of compositional phenomena in oil and gas reservoirs [64]. This assumption holds as long as exists an appropriate equation of state calibration for the PVT fluid behavior. Adding flexibility to our model, and making it more general for Black-oil type and fully-compositional fluids, the problem of thermodynamic properties determination was decoupled from the conservation equations, being solved sequentially after determination of primary variables. As a general description, thermodynamic equilibrium can be defined mathematically as the equality of fugacities, ${\widehat{f}}_{k,p}$, of each component k in a mixture for all $n_p$ phases [65].

$$\begin{array}{cc}\hfill {\widehat{f}}_{k,1}={\widehat{f}}_{k,2}=\dots & ={\widehat{f}}_{k,p}=\cdots ={\widehat{f}}_{k,n_p}\hfill \\ & \forall (p,k)\in \mathcal{P}\times \mathcal{K}.\hfill \end{array}$$

(8)

Fugacity is solved by one of the several equations of state available for hydrocarbon modeling, attached to the following inner constraints:

$$\sum _{k=1}^{n_k}{w}_{k,p}=1;\sum _{p=1}^{n_p}{F}_p=1;\sum _{i=1}^{n_k}{z}_k=1$$

(9)

where $F_p$ is the molar fraction of phase p in a point, and $z_k$ is the global molar fraction of component k in that particular spatial location.

This set of equations is internally solved at each instant and space location, applying the Peng & Robinson [66] equation of state. The solution of such a system is then sequentially coupled to the transport equations, assuring a precise calculation of the involved thermodynamic properties. Considering three pseudo-components gas, oil, and water for an extended black-oil model, a modified version of the Wang [67] approach was used to translate the tabulated PVT information to the compositional formulation.

2.3.2. Capillary Pressure

Capillary pressure equations state a pressure discontinuity between two competing phases, a wetting phase $P_{\mathrm{w}}$, and a non-wetting phase $P_{\mathrm{nw}}$ [68,69,70]. Capillary pressure relationship is expressed as a general functional as follows for two arbitrary wetting and non-wetting phases:

$$P_{\mathrm{c}}={\mathcal{L}}_{\alpha }\left(P_{\mathrm{nw}},P_{\mathrm{w}}\right)$$

(10)

where, $n_p-1$ capillary pressure equations, one per binary system $\alpha$, are required in a system of multiple $n_p$ phases. The pressure difference between phases is considered for the transport equations but is assumed to be insignificant for thermodynamic equilibrium calculations.

2.3.3. Volume Constraint

The volume constraint equation states that the sum of the volume occupied for each phase $V_p$ is equal to the porous space. Volume constraint equation is as follows:

$$\sum _{p\in \mathcal{P}}\frac{V_p}{V_T}-\varphi =0$$

(11)

2.4. Second-Order Anisotropic Elliptic Operator

Most of the transport models can be written in terms of the second-order anisotropic elliptic operator defined as $\mathfrak{D}\left(\Lambda ,\Psi \right):=\mathrm{div}\left(\Lambda \nabla \Psi \right)$, where $\Lambda (\mathbf{x},t)$ is a tensorial property and $\Psi (\mathbf{x},t)$ the scalar function of interest. For an advective model, $\Lambda =-{\lambda }_p\mathbb{K}$ and $\Psi ={\Phi }_p$.

3. Numerical Model

In this section, we present the discretization for the continuity Equation (5) for advection, using an unstructured finite volume method (UFVM) maintaining the second order of convergence, considering the correction for anisotropy non-orthogonality and skewness.

3.1. Time Discretization and Mesh Definition

For computing an approximation of the solution in time, a discretization of the simulation interval $\Theta$ is required first. Let ${\mathcal{T}}_t:=\left\{t_0,t_1,\dots ,t_f\right\}$ a partition of $\Theta$. Each temporal element or simulation interval is defined as $\Delta t_n:=(t_n,t_{n+1})$, for $n\in \left\{0,1,\dots ,k-1\right\}$.

For spatial discretization, let $\Omega \in {\mathbb{R}}^3$ the physical domain of interest and ${\Omega }_j\subset \Omega$, with $j\in \mathcal{F}:=\left\{1,2,\dots ,N_{\Omega }\right\}$, open, bounded and simply connected subsets in ${\mathbb{R}}^3$, such that ${\Omega }_i\cap {\Omega }_j=\varnothing$, for all $i\ne j$ and also

$$\bigcup _{j=1}^{N_{\Omega }}{\overline{\Omega }}_j=\overline{{\Omega }^h}\underset{N_{\Omega }\to \infty }{⟶}\overline{\Omega },$$

(12)

where $\overline{{\Omega }^h}$ is the discrete domain, which approximate $\Omega$[71]. The subsets of the partition ${\Omega }_j$ are denoted as elements or cells [72]. Note that in the context of this definition, when considering the notion of an open set and its closure, is implicitly distinguishing between the interior of the elements ${\Omega }_j$ and their respective boundaries $\partial {\Omega }_j$, which are divided in a finite number of simply connected faces whose boundaries are also known as edges [44,45,73]. Two cells ${\Omega }_j$ and ${\Omega }_l$ are neighboring or adjacent if both share a face f, that is, $\overline{{\Omega }_j}\cap \overline{{\Omega }_j}=f\ne \varnothing$.

3.2. Finite Volume Integral Formulation

The integral form of the continuity equations, in terms of the second-order anisotropic elliptic operator, for advection transport, is obtained by integrating Equation (5) over an arbitrary element ${\Omega }_j$ and time interval $\Delta t_n$ as:

$$\begin{array}{l}\underset{\mathrm{Accumulation}\mathrm{term}.}{\underbrace{{\int }_{\Delta t_n}{\int }_{{\Omega }_j}\frac{\partial N_k}{\partial t}dVdt}}+\underset{\mathrm{Transport}\mathrm{term}}{\underbrace{{\displaystyle \sum _{p=1}^{n_p}}{\int }_{\Delta t_n}{\int }_{{\Omega }_j}\mathfrak{D}\left(w_{k,p}{\lambda }_p\mathbb{K},{\Phi }_p\right)dVdt}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\underset{\mathrm{Source}/\mathrm{sink}\mathrm{term}.}{\underbrace{{\int }_{\Delta t_n}{\int }_{{\Omega }_j}{\dot{q}}_{k}dVdt}}=0,\forall k\in \left\{1,2,\dots ,N_k\right\}\end{array}$$

(13)

The discretization of each of the terms of Equation (13) is done separately, resulting in the system of nonlinear algebraic equations, which is solved by the Newton method.

Accumulation Term

For the discretization of time derivative is assumed that each element geometry remains constant along with the simulation, and the fluid properties are homogeneous over each element [74]. By applying the Fubini theorem and the fundamental theorem of calculus, the discretization form of the accumulation term is obtained as:

$${\int }_{\Delta t_n}{\int }_{{\Omega }_j}\frac{\partial N_{k,j}}{\partial t}dVdt={\int }_{{\Omega }_j}\left({\int }_{\Delta t_n}d{N}_{k,j}\right)dV={\left.N_{k,j}\right|}_{\Delta t_n}\left|{\Omega }_j\right|,$$

(14)

where $\left|{\Omega }_j\right|$ is the measure of ${\Omega }_j$.

3.3. Advective Term

This work focuses mainly on the discretization of this term, considering full permeability tensors and irregular polyhedral elements. The first step is to apply the Green-Gauss theorem to the elliptic operator as follows:

$$\begin{array}{l}{\int }_{{\Omega }_j}\mathfrak{D}\left(w_{k,p}{\lambda }_p\mathbb{K},{\Phi }_p\right)dV={\int }_{\partial {\Omega }_j}{\lambda }_p{w}_{k,p}〈\nabla {\Phi }_p,{\mathbb{K}}^{\top }\mathbf{n}〉dS\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\int }_{\partial {\Omega }_j}{\lambda }_p{w}_{k,p}\frac{\partial {\Phi }_p}{\partial \mathbf{w}}dS,\end{array}$$

(15)

where $〈·,·〉$ denotes the inner product in ${\mathbb{R}}^n$, $\mathbf{n}=\mathbf{n}\left(\mathbf{x}\right)$ is the unit vector normal to surface $\partial {\Omega }_j$, and ${\mathbf{w}}_f$ is the co-normal vector defined as ${\mathbf{w}}_f={\mathbb{K}}_f^{\top }{\mathbf{n}}_f$ in which the anisotropy effects reside. To deal with potential field discontinuity due to heterogeneity, the permeability tensor in f is computed in harmonic sense as [32]:

$${\mathbb{K}}_f={\displaystyle \frac{\left|{\mathbf{x}}_f^i-{\mathbf{x}}_j^i\right|+\left|{\mathbf{x}}_f^i-{\mathbf{x}}_l^i\right|}{\frac{\left|{\mathbf{x}}_f^i-{\mathbf{x}}_j^i\right|}{{\mathbb{K}}_j^i}+\frac{\left|{\mathbf{x}}_f^i-{\mathbf{x}}_l^i\right|}{{\mathbb{K}}_l^i}}},$$

(16)

where $i=1,2,3,$ is the coordinate system direction.

The surface integral in Equation (15) can be interpreted as the net flux of moles of component k, in phase p, across the boundary of the element. Then, as the control volume is a polyhedron defined by a finite number of faces, the surface integral also can be decomposed in a sum over his faces, and the flux en each one is approximated with a second-order midpoint rule [75] as:

$${\int }_{\partial {\Omega }_j}{\lambda }_p{w}_{k,p}\frac{\partial {\Phi }_p}{\partial \mathbf{w}}dS=\sum _f{\lambda }_{p,f}{w}_{k,p,f}\frac{\partial {\Phi }_{p,f}}{\partial {\mathbf{w}}_f}{A}_f,$$

(17)

where $A_f$ is the face area, and in general, the subscript f indicates that the property, whether scalar or vector, is evaluated at the face centroid ${\mathbf{x}}_f$. For obtaining the system of algebraic equations, it is necessary to find an approximation of the directional derivative $\frac{\partial {\Phi }_{p,f}}{\partial {\mathbf{w}}_f}$, in terms of the potentials evaluated in the cell centroids.

Let ${\Omega }_j$ and ${\Omega }_l$ adjacent elements through the face f with centroids ${\mathbf{x}}_j$ and ${\mathbf{x}}_l$ respectively and the vector that connects these centroids, defined as ${\mathbf{d}}_f:={\mathbf{x}}_l-{\mathbf{x}}_j$. Then, vectors ${\mathbf{d}}_f$ and ${\mathbf{w}}_f$ are expressed as a linear combination of the elements of the face orthonormal basis ${\mathcal{B}}_f:=\left\{{\mathbf{n}}_f,{\mathbf{T}}_f,{\mathbf{T}}_f^{\prime }\right\}$ as

$$\begin{array}{cc}\hfill {\mathbf{d}}_f=& \langle {\mathbf{d}}_f,{\mathbf{n}}_f\rangle {\mathbf{n}}_f+\langle {\mathbf{d}}_f,{\mathbf{T}}_f\rangle {\mathbf{T}}_f+〈{\mathbf{d}}_f,{\mathbf{T}}_f^{\prime }〉{\mathbf{T}}_f^{\prime }\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\mathbf{w}}_f=& \langle {\mathbf{w}}_f,{\mathbf{n}}_f\rangle {\mathbf{n}}_f+\langle {\mathbf{w}}_f,{\mathbf{T}}_f\rangle {\mathbf{T}}_f+〈{\mathbf{w}}_f,{\mathbf{T}}_f^{\prime }〉{\mathbf{T}}_f^{\prime }\hfill \end{array}$$

(19)

where ${\mathbf{T}}_f$ and ${\mathbf{T}}_f^{\prime }$ are unit vectors tangential to the face such that ${\mathbf{T}}_f\perp {\mathbf{T}}_f^{\prime }$. The geometrical elements defined here are shown schematically in Figure 1. Getting an expression for the normal vector from (18) and replacing it on (19), a vector equation for ${\mathbf{w}}_f$ is obtained in terms of ${\mathbf{T}}_f,{\mathbf{T}}_f^{\prime }$ and ${\mathbf{d}}_f$ as

$${\mathbf{w}}_f={\alpha }_f{\mathbf{d}}_f+{\beta }_f{\mathbf{T}}_f+{\beta }_f^{\prime }{\mathbf{T}}_f^{\prime },$$

(20)

where the coefficients are

$${\alpha }_f={\displaystyle \frac{〈{\mathbf{w}}_f,{\mathbf{n}}_f〉}{〈{\mathbf{d}}_f,{\mathbf{n}}_f〉}},$$

(21)

and

$$\begin{array}{cc}\hfill {\beta }_f& =-{\alpha }_f〈{\mathbf{d}}_f,{\mathbf{T}}_f〉+〈{\mathbf{w}}_f,{\mathbf{T}}_f〉,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\beta }_f^{\prime }& =-{\alpha }_f〈{\mathbf{d}}_f,{\mathbf{T}}_f^{\prime }〉+〈{\mathbf{w}}_f,{\mathbf{T}}_f^{\prime }〉.\hfill \end{array}$$

(23)

Besides, the product ${\alpha }_f{A}_f$ is known as shape factor in classical reservoir simulation literature, and it coincides with the formulation of Islam [74], Abou-Kassem [68] and Aziz [69], for rectilinear orthogonal grids. Thus, applying the vector decomposition (20) to the directional derivative at (17) gives the following expression for the advective term:

$$\begin{array}{l}{\int }_{{\Omega }_j}\mathfrak{D}\left(w_{k,p}{\lambda }_p\mathbb{K},{\Phi }_p\right)dV=\sum _f\underset{\mathrm{Orthogonal\; flux}}{\underbrace{{\lambda }_{p,f}{w}_{k,p,f}{\alpha }_f\frac{\partial {\Phi }_{p,f}}{\partial {\mathbf{d}}_f}{A}_f}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\sum _f\underset{\mathrm{Tan}\mathrm{gential}\mathrm{flux}}{\underbrace{{\lambda }_{p,f}{w}_{k,p,f}\left({\beta }_f\frac{\partial {\Phi }_{p,f}}{\partial {\mathbf{T}}_f}+{\beta }_f^{\prime }\frac{\partial {\Phi }_{p,f}}{\partial {\mathbf{T}}_f^{\prime }}\right)A_f}}.\end{array}$$

(24)

The tangential component takes account of skewness distortion error, as shown below:

First, using a Taylor series expansion for the potential centered at the face centroid, evaluated at ${\mathbf{x}}_j$ and ${\mathbf{x}}_l$, yield to expressions for ${\Phi }_{p,j}$ and ${\Phi }_{p,l}$ respectively, in terms of potential partial derivatives as

$${\Phi }_{p,j}=\sum _{k=0}^m\frac{1}{k!}{〈{\mathbf{d}}_{fj},\nabla 〉}^k{\Phi }_p\left({\mathbf{x}}_f\right)+R_{fj},$$

(25)

$${\Phi }_{p,l}=\sum _{k=0}^m\frac{1}{k!}{〈{\mathbf{d}}_{fl},\nabla 〉}^k{\Phi }_p\left({\mathbf{x}}_f\right)+R_{fl},$$

(26)

where $R_{fj}$ and $R_{fl}$ are the residual truncation errors for the Taylor series, and the face-to-cell director vectors are defined as ${\mathbf{d}}_{fj}:={\mathbf{x}}_f-{\mathbf{x}}_j$ and ${\mathbf{d}}_{fl}:={\mathbf{x}}_f-{\mathbf{x}}_l$. Now, subtracting the expressions (25) and (26) result in:

$${\Phi }_{p,j}-{\Phi }_{p,l}=〈{\mathbf{d}}_f,\nabla {\Phi }_{p,f}〉+\sum _{k=2}^m\frac{1}{k!}\left\{{〈{\mathbf{d}}_{fj},\nabla 〉}^k-{〈{\mathbf{d}}_{fj},\nabla 〉}^k\right\}{\Phi }_{p,f}+R_{fl},$$

(27)

where $R_{fl}=R_{fj}-R_{fl}$, and therefore

$$\frac{\partial {\Phi }_{p,f}}{\partial {\mathbf{d}}_f}=\left({\Phi }_{p,l}-{\Phi }_{p,j}\right)-{\delta }_{p,f}^{\left(m\right)}+\mathcal{O}\left(h^{m+1}\right),$$

(28)

in which

$${\delta }_{p,f}^{\left(m\right)}=\sum _{k=2}^m\frac{1}{k!}\left\{{〈{\mathbf{d}}_{fj},\nabla 〉}^k-{〈{\mathbf{d}}_{fl},\nabla 〉}^k\right\}{\Phi }_{p,f}.$$

(29)

This work is limited to the case $m=2$ in which the operator ${〈\mathbf{d},\nabla 〉}^2$ is expressed as

$${〈\mathbf{d},\nabla 〉}^2={\left(\sum _{i=1}^d{d}_i\frac{\partial }{\partial x_i}\right)}^2=\sum _{k=1}^d{d}_k^2\frac{{\partial }^2}{\partial x_k^2}+2\sum _{j=1}^d\sum _{i=1}^{j-1}{d}_i{d}_j\frac{{\partial }^2}{\partial x_i\partial x_j}.$$

(30)

Finally, the advective term discretization is

$$\begin{array}{l}{\int }_{{\Omega }_j}\mathfrak{D}\left(w_{k,p}{\lambda }_p\mathbb{K},{\Phi }_p\right)dV={\displaystyle \sum _f}{\lambda }_{p,f}{w}_{k,p,f}{\alpha }_f\left({\Phi }_{p,l}-{\Phi }_{p,j}\right)A_f\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _f}{Q}_{p,f}^{\mathrm{NO}}-{\displaystyle \sum _f}{Q}_{p,f}^{\mathrm{skew}}.\end{array}$$

(31)

where

$$Q_{p,f}^{\mathrm{NO}}={\lambda }_{p,f}{w}_{k,p,f}\left({\beta }_f\frac{\partial {\Phi }_{p,f}}{\partial {\mathbf{T}}_f}+{\beta }_f^{\prime }\frac{\partial {\Phi }_{p,f}}{\partial {\mathbf{T}}_f^{\prime }}\right)A_f,$$

(32)

$$Q_{p,f}^{\mathrm{skew}}={\lambda }_{p,f}{w}_{k,p,f}{\alpha }_f{\delta }_{p,f}^{\left(2\right)}{A}_f.$$

(33)

3.4. Source/Sink Term

This term is directly associated with the production/injection rate of the component k in the domain. The production rate ${\dot{q}}_k$ can be understood as a function, of space and time, whose second-order discretization consists only of the application of the mean value theorem [74] as follows:

$${\int }_{\Delta t_n}{\int }_{{\Omega }_j}{\dot{q}}_{k}d\Omega dt\approx \left|\Delta t_n\right|\left|{\Omega }_j\right|{\dot{q}}_k({\mathbf{x}}_j,t)$$

(34)

However, in engineering applications, an approach of type (34) is not practical in most cases, mainly because this term is used for including injection/production well modeling, implying the function ${\dot{q}}_k$ is highly discontinuous on the domain. Additionally, component k can be produced in different phases and is associated with the volumetric rates of each phase $Q_{p,j}$, through the molar fraction $w_{k,p}$. In this case, additional equations are required [76]. Peaceman extended well model [77] is used in this work. In this way, the source/sink term is thoroughly characterized as:

$${\int }_{{\Omega }_j}{\dot{q}}_{k}dV={\mathrm{WI}}_j\sum _{p=1}^{n_p}{w}_{k,p,j}{\rho }_{p,j}{\lambda }_{p,j}\left(P_{p,j}-P_{\mathrm{wf},j}\right).$$

(35)

where $P_{\mathrm{wf},j}:=P_{\mathrm{wf},w}+{\gamma }_w\left(z_j-z_w\right)$ is the bottom hole well pressure and ${\mathrm{WI}}_j$ is the well index upscaling factor of element ${\Omega }_j$ defined as

$${\mathrm{WI}}_j:=\frac{2\pi h_j\sqrt{{\mathbb{K}}_{j,1}{\mathbb{K}}_{j,2}}}{ln\left(r_e/r_0\right)+S_j}.$$

(36)

The coefficient $S_j$ is the skin factor [78], which allows for considering permeability reductions in well-bore damage. The geometrical elements used here are shown schematically in Figure 2.

Also, it is worth mentioning that the variable $Q_{p,j}$ is usually the determining variable for EOR economic models. For this reason, this variable is essential in the validation of the method (see numerical example Section 6.1).

Time Integration

For the time integrals, a generalized Crank Nicholson scheme [50], also known as implicit-explicit or IMEX method [51,52], is used. This scheme can be seen as a weighted trapezoidal rule [79] applied on a continuous function F as

$${\int }_{\Delta t_n}Fdt=\left[\vartheta F\left(t_n\right)+(1-\vartheta )F\left(t_{n+1}\right)\right]\left|\Delta t_n\right|,$$

(37)

where $\left|\Delta t_n\right|$ is the length of the time interval. The classical implicit and explicit schemes can be obtained directly by setting the parameter $\vartheta$ to 0 or 1, respectively. To simplify the notation, write $F\left(t_n\right):=F^n$.

4. Weighted Least Squares

4.1. General Formulation

Let $\mathcal{S}:={\left\{{\mathbf{x}}_i\right\}}_{i\in {\mathcal{F}}_s}$ a family of nodes in space, where ${\mathcal{F}}_s$ is the index collection. Neglecting the truncation error, consider Taylor series centered at a point of interest ${\mathbf{x}}_{\lambda }$ which produces the following overdetermined system of linear equations obtained by evaluating for each ${\mathbf{x}}_i$:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \left\{\begin{array}{c}\varphi \left({\mathbf{x}}_1\right)=\varphi \left({\mathbf{x}}_{\lambda }\right)+\sum _{k=1}^m\frac{1}{k!}{〈{\mathbf{x}}_1-{\mathbf{x}}_{\lambda },\nabla 〉}^k\varphi \left({\mathbf{x}}_{\lambda }\right),\hfill \\ \varphi \left({\mathbf{x}}_2\right)=\varphi \left({\mathbf{x}}_{\lambda }\right)+\sum _{k=1}^m\frac{1}{k!}{〈{\mathbf{x}}_2-{\mathbf{x}}_{\lambda },\nabla 〉}^k\varphi \left({\mathbf{x}}_{\lambda }\right),\hfill \\ ⋮⋮⋮\hfill \\ \varphi \left({\mathbf{x}}_i\right)=\varphi \left({\mathbf{x}}_{\lambda }\right)+\sum _{k=1}^m\frac{1}{k!}{〈{\mathbf{x}}_i-{\mathbf{x}}_{\lambda },\nabla 〉}^k\varphi \left({\mathbf{x}}_{\lambda }\right),\hfill \\ ⋮⋮⋮\hfill \\ \varphi \left({\mathbf{x}}_{n_s}\right)=\varphi \left({\mathbf{x}}_{\lambda }\right)+\sum _{k=1}^m\frac{1}{k!}{〈{\mathbf{x}}_{n_s}-{\mathbf{x}}_{\lambda },\nabla 〉}^k\varphi \left({\mathbf{x}}_{\lambda }\right).\hfill \end{array}\right.\end{array}\end{array}$$

(38)

where the number of nodes $n_s$, are greater than the number of unknown partial derivatives until order m of the scalar function $\varphi$ at ${\mathbf{x}}_{\lambda }$. Moreover, the equation system (38) can be written in matrix form as

$${\mathbb{M}}_{n_s\times n_I}{\dot{\mathbf{\Phi }}}_{n_I\times 1}={\mathbf{b}}_{n_s\times 1},$$

(39)

where $n_I$ is the number of unknowns. The weighted least squares [80,81] problem (WLS) for this formulation is: Find ${\dot{\mathbf{\Phi }}}^h$, approximated solution of (39), such that

$${\dot{\mathbf{\Phi }}}^h=\underset{\dot{\mathbf{\Phi }}\in {\mathbb{R}}^{n_I}}{min}\left\{{∥\mathbf{r}(w,\mathbf{b},\dot{\mathbf{\Phi }})∥}_2\right\}=\underset{\dot{\mathbf{\Phi }}\in {\mathbb{R}}^{n_I}}{min}\left\{\sum _{i=1}^{n_s}\sqrt{w_i{\mathbf{r}}_i^2}\right\}=\underset{\dot{\mathbf{\Phi }}\in {\mathbb{R}}^{n_I}}{min}\left\{\sum _{i=1}^{n_s}{w}_i{\mathbf{r}}_i^2\right\}.$$

(40)

where

$$w_i=\frac{1}{{∥{\mathbf{x}}_i-{\mathbf{x}}_{\lambda }∥}^p}.$$

(41)

Some authors take the p power equal to 1 or 2 [30,32,82].

WLS Solution

There are two equivalent alternatives to solve the WLS problems [83]. The first is to solve the linear system of normal equations $\left({\mathbb{M}}^T\mathbb{W}\mathbb{M}\right)\dot{\mathbf{\Phi }}={\mathbb{M}}^T\mathbb{W}\mathbf{b}$, where the solution is given by:

$${\dot{\mathbf{\Phi }}}^h={\mathbb{M}}_w^{†}\mathbf{b},$$

(42)

where the matrix ${\mathbb{M}}_w^{†}\in {\mathbb{R}}^{n\times m}$, is the Moore-Penrose weighted pseudo-inverse of $\mathbb{M}$, defined as:

$${\mathbb{M}}_w^{†}={\left({\mathbb{M}}^T\mathbb{W}\mathbb{M}\right)}^{-1}{\mathbb{M}}^T\mathbb{W}.$$

(43)

The second alternative is to find an orthogonal projector $\mathbf{P}\in {\mathbb{R}}^{n_s\times n_s}$, on the range of $\mathbb{M}$ such that $\mathbf{P}\mathbf{b}=\mathbb{M}\mathbf{\Phi }$.

4.2. Stencil Definition

Depending on the selection of the point of interest ${\mathbf{x}}_{\lambda }$, there are two possible definitions of the stencil for WLS application in a grid: cell stencils, where ${\mathbf{x}}_{\lambda }={\mathbf{x}}_j$, and face stencils with ${\mathbf{x}}_{\lambda }={\mathbf{x}}_f$.

4.2.1. Cell Stencils

Let $S_j$ the set of nodes for stencil definition of an arbitrary cell ${\Omega }_j$. The most natural way to select these nodes is to take advantage of the connectivity of the cell through their faces, choosing the nodes of the neighboring elements, defining what is called a compact stencil [84]. Also, the stencil can be extended considering, not only the elements neighboring through the faces but those that have common nodes. For example, see Figure 3, in which the compact stencil is ${\mathcal{S}}_j=\left\{{\mathbf{x}}_{s_1},{\mathbf{x}}_{s_2},\dots ,{\mathbf{x}}_{s_6}\right\}$, and the extended stencil is given by ${\mathcal{S}}_j\cup \left\{{\mathbf{x}}_{s_7},{\mathbf{x}}_{s_8}\right\}$.

4.2.2. Face Stencils

Let ${\mathcal{S}}_j$ and ${\mathcal{S}}_l$ nodes for either compact or extended stencils of neighboring cells ${\Omega }_j$ and ${\Omega }_l$. The nodes for the face stencil are defined as

$${\mathcal{S}}_f={\mathcal{S}}_j\cup {\mathcal{S}}_l\cup \left\{{\mathbf{x}}_j,{\mathbf{x}}_l\right\}.$$

(44)

The typical definition of the face stencil is illustrated in Figure 4, in which ${\mathcal{S}}_f=$$\left\{{\mathbf{x}}_j,{\mathbf{x}}_l,{\mathbf{x}}_{s_1},{\mathbf{x}}_{s_2},\dots ,{\mathbf{x}}_{s_{10}}\right\}$.

This stencil construction, in terms of the connectivity of both cells, has its raison d’être; it guarantees a single stencil per face. This fact corrects a possible source of error in the formulation of Loudyi [32] in which the face stencil is relative to each cell. This source of error causes numerical instabilities in Newton’s method when solving the resultant algebraic system from finite volume discretization, mainly because of the difference between fluxes at each cell residual, which can lead to non-convergence when the algebraic equations system is ill-conditioned; a condition that is not unusual in the modeling of multiphase flow in porous media [85].

4.2.3. Stencils at Boundaries

If the point of interest ${\mathbf{x}}_{\lambda }$ lies in a boundary cell with Neumann or Dirichlet faces, then some ${\dot{\mathbf{\Phi }}}^h$ information is known depending on its boundary type. This information allows for simplifying the overdetermined system (38) by adding the boundary face nodes to the stencil, as is shown schematically in Figure 5.

4.3. Gradient Reconstruction

The application of (38) to approximate the gradient derivatives for the potential at time n, in the point of interest ${\mathbf{x}}_{\lambda }$ is

$$\left[\begin{array}{cccc}1& \Delta x_1& \Delta y_1& \Delta z_1\\ 1& \Delta x_2& \Delta y_2& \Delta z_2\\ ⋮& ⋮& ⋮& ⋮\\ 1& \Delta x_i& \Delta y_i& \Delta z_i\\ ⋮& ⋮& ⋮& ⋮\\ 1& \Delta x_{n_s}& \Delta y_{n_s}& \Delta z_{n_s}\end{array}\right]\left[\begin{array}{c}\Phi \left({\mathbf{x}}_{\lambda }\right)\\ \frac{\partial \Phi }{\partial x}\left({\mathbf{x}}_{\lambda }\right)\\ \frac{\partial \Phi }{\partial y}\left({\mathbf{x}}_{\lambda }\right)\\ \frac{\partial \Phi }{\partial z}\left({\mathbf{x}}_{\lambda }\right)\end{array}\right]=\left[\begin{array}{c}\Phi \left({\mathbf{x}}_1\right)\\ \Phi \left({\mathbf{x}}_2\right)\\ ⋮\\ \Phi \left({\mathbf{x}}_i\right)\\ ⋮\\ \Phi \left({\mathbf{x}}_{n_s}\right)\end{array}\right]$$

(45)

where the coefficients $\Delta x_i$, $\Delta y_i$, and $\Delta z_i$ are the components of the vector ${\mathbf{x}}_i-{\mathbf{x}}_{\lambda }$. It is essential to take account of the sign of each coefficient because it considers the orientation of the Taylor series radius. If a cell stencil is used, the procedure requires an additional interpolation [86,87] in order to obtain the gradient at the face required by (24). The most common way to do this is the linear interpolation between the gradient values at neighboring cell centroids, given by

$$\nabla {\Phi }_p^n\left({\mathbf{x}}_f\right)=f_x\nabla {\Phi }_p^n\left({\mathbf{x}}_l\right)+(1-f_x)\nabla {\Phi }_p^n\left({\mathbf{x}}_j\right),$$

(46)

where

$$f_x=\frac{∥{\mathbf{x}}_j-{\mathbf{x}}_f∥}{∥{\mathbf{x}}_j-{\mathbf{x}}_l∥}.$$

(47)

In this case, to guarantee the second-order convergence when (46) is used, an additional correction for skewness is required [88,89].

4.4. Hessian Reconstruction

The application of (38) to approximate the Hessian for the flow potential at time n, in the point of interest ${\mathbf{x}}_{\lambda }$, required for computing the skewness flow correction (33), is:

$$\left[\begin{array}{cccccc}1& \Delta x_1& \dots & \frac{1}{2}\Delta x_1^2& \dots & \Delta y_1\Delta z_1\\ 1& \Delta x_2& \dots & \frac{1}{2}\Delta x_2^2& \dots & \Delta y_2\Delta z_2\\ ⋮& ⋮& & ⋮& & ⋮\\ 1& \Delta x_i& \dots & \frac{1}{2}\Delta x_i^2& \dots & \Delta y_i\Delta z_i\\ ⋮& ⋮& & ⋮& & ⋮\\ 1& \Delta x_{n_s}& \dots & \frac{1}{2}\Delta x_{n_s}^2& \dots & \Delta y_{n_s}\Delta z_{n_s}\end{array}\right]\left[\begin{array}{c}{\Phi }^n\left({\mathbf{x}}_{\lambda }\right)\\ \frac{\partial \Phi }{\partial x}\left({\mathbf{x}}_{\lambda }\right)\\ \frac{\partial \Phi }{\partial y}\left({\mathbf{x}}_{\lambda }\right)\\ \frac{\partial \Phi }{\partial z}\left({\mathbf{x}}_{\lambda }\right)\\ \frac{{\partial }^2\Phi }{\partial x^2}\left({\mathbf{x}}_{\lambda }\right)\\ \frac{{\partial }^2\Phi }{\partial y^2}\left({\mathbf{x}}_{\lambda }\right)\\ \frac{{\partial }^2\Phi }{\partial z^2}\left({\mathbf{x}}_{\lambda }\right)\\ \frac{{\partial }^2\Phi }{\partial x\partial x}\left({\mathbf{x}}_{\lambda }\right)\\ \frac{{\partial }^2\Phi }{\partial x\partial z}\left({\mathbf{x}}_{\lambda }\right)\\ \frac{{\partial }^2\Phi }{\partial y\partial z}\left({\mathbf{x}}_{\lambda }\right)\end{array}\right]=\left[\begin{array}{c}\Phi \left({\mathbf{x}}_1\right)\\ \Phi \left({\mathbf{x}}_2\right)\\ ⋮\\ \Phi \left({\mathbf{x}}_i\right)\\ ⋮\\ \Phi \left({\mathbf{x}}_{n_s}\right)\end{array}\right]$$

(48)

5. Semi-Implicit Implementation

The formulation above is the guideline for FlowTraM semi-implicit implementation, which uses the Newton-Raphson method for algebraic equations linearization. This implementation is presented in Algorithm 1. The core steps in the algorithm, which allow for managing non-orthogonal geometries, are cell stencils and potential stencils calculation, required for multipoint flux approximation (MPFA); TPFA and MPFA calculation, and residual correction.

In cell stencil calculations, the matrix operator is created, which depends on its geometrical information, like $\Delta x$, $\Delta y$, $\Delta z$. This process is done once per simulation unless changes exist in geometry through time. If absolute permeability remains constant, alpha and beta shape factors can be calculated to this stage.

Subsequently, in potential stencil calculations, gradients per cell are reconstructed by solving the overdetermined linear system given by cell-matrix operators and time k potentials $\left({\Phi }_{\Omega }^k\right)$. Later, face gradients are obtained via interpolation using the centroid distance between two neighboring cells. MPFA calculation refers to non-orthogonal flux calculation.

Finally, a correction for residual is obtained by accumulating non-orthogonal fluxes for each face of a cell. The resulting correction is added to residual in Newton iteration, and a corrected solution is obtained by applying a preconditioned iterative solver like GMRES with ILU to Newton’s linear system.

The proposal above is called semi-implicit because, for MPFA calculation in gradient reconstruction, time k potentials are used $\left({\Phi }_{\Omega }^k\right)$. Nonetheless, these fluxes are omitted in Jacobian matrix calculation but added as a correction to TPFA residual. It is worth noting that, with this approach, the sparse nature of the Jacobian matrix is maintained.

Algorithm 1 Semi-Implicit approach Pseudocode

6. Numerical Experiments and Discussion

The proposed methodology was proven successful in Colciencias 273-2017, 272-2017, 064-2018 projects. These were aimed at the development of technologies for EOR processes, to maximize usages of oil reserves in Colombia, thus, demonstrating our model versatility. Nevertheless, in this work, three numerical experiments are presented. The first one is for validating the methodology; the second one presents a pseudo-compositional implementation to reproduce a Black-oil study case, and the third one is for demonstrating practical usage and applying it to a field case.

The first experiment is an adaptation of the Third Comparative Solution Project (SPE3) [90], simulated on a mesh with high non-orthogonalities. The obtained results allow for validating the formulation and help to illustrate the improvement when corrective fluxes are introduced clearly. The second numerical experiment is a three-phase conning study in a radial mesh, which is used to compare multiple commercial simulators and validate newer ones. The third experiment is a history matching of parameters for a Colombian oil field at the llanos basin, located at Piedemonte llanero, a realistic simulation.

6.1. SPE Comparative Solution Project 3

This simulation problem consists of a gas re-circulation, for pressure maintenance, in a heterogeneous reservoir of condensate gas. Dry gas is injected at 4.7 MMSCF/day in the upper area of the reservoir and sustained for ten years, while the lower area of the reservoir has a producing well conditioned to 6.2 MMSCF/day. Fluid properties are listed in

Table 1. Properties of condensate gas evaluated in the third SPE comparative project  [90].

k	Acentric Factor	${\mathit{P}}_{\mathbf{c}}$ (psi)	${\mathit{T}}_{\mathbf{c}}$ (${}^{\circ }$R)	${\mathit{w}}_{\mathbf{mol}}$	${\mathit{V}}_{\mathbf{cr}}$ (ft${}^3$/lbm)
C${}_1$	0.0130	588.0	350.00	16.04	0.0988
C${}_2$	0.0986	708.5	549.77	30.07	0.0788
C${}_3$	0.1524	617.5	665.82	44.10	0.0727
C${}_4$	0.2010	550.8	765.36	58.12	0.0703
C${}_5$	0.2539	489.7	845.28	72.15	0.0690
C${}_6$	0.3007	439.8	914.22	86.18	0.0684
C${}_{7-9}$	0.3613	385.9	1032.21	114.43	0.0652
C${}_{10-11}$	0.4501	340.8	1148.02	144.83	0.0630
C${}_{12-14}$	0.5339	293.8	1234.35	177.78	0.0622
C${}_{15+}$	0.7244	184.5	1347.00	253.63	0.0610

.

The orthogonal and the adapted non-orthogonal grids used to simulate this case are shown in Figure 6. Moreover, Figure 6 shows the cell average of the K-orthogonality angle, defined for each face by

$$cos{\theta }_f=\frac{〈{\mathbf{d}}_f,{\mathbf{w}}_f〉}{∥{\mathbf{d}}_f∥∥{\mathbf{w}}_f∥}.$$

(49)

The pressure distribution for both geometries is presented in Figure 7 for different time steps. Due to the pressure drop caused during gas production, condensate oil banks continuously appear in the producing blocks. These banks become mobile when they reach critical saturation, which in this case, is 16%. Thus, once formed, they can only disappear by contact with very light gases or by an additional depressurization that causes the gradual volatilization of its lighter components.

The oil produced in this reservoir is governed by condensation of intermediate components from the gas phase once it reaches the surface and oil mobilization once the critical saturation is reached. Notice that the most significant variable in reservoir simulation is the oil rate computed at the surface. Figure 8 shows surface oil production for the different test geometries. The obtained results validate the implemented formulation in FlowTraM and show the improvement in convergence obtained by introducing the proposed first and second-order corrective flows.

Convergence Analysis

In order to validate the effect of the correction schemes, we performed a convergence analysis using a significantly refined orthogonal mesh geometry as a reference solution, in the absence of an analytical solution for this problem. The spatial discretization of the reference mesh is represented by $NX=100$, $NY=100$, $NZ=40$, whose result is shown in Figure 9.

The analysis makes use of the normalized experimental error in norm $L^1$, defined as

$$e_h^{L^1}=\frac{\parallel {\Phi }_{\mathrm{ref}}-{\Phi }_h{\parallel }_{L^1(\Omega )}}{\parallel {\Phi }_{\mathrm{ref}}{\parallel }_{L^1(\Omega )}},$$

(50)

where ${\Phi }_{\mathrm{ref}}$ is the resolved potential in the reference mesh, and ${\Phi }_h$ is the solution in a skewed non-orthogonal mesh of size h, which can be defined as

$$h=\underset{j\in \mathcal{F}}{max}\left|{\Omega }_j\right|.$$

(51)

These meshes were constructed analogously to the one shown in Figure 6. The results are shown in

Table 2. Convergence analysis results (UC:Uncorrected; FOC:First Order Corrected; SOC: Second Order Corrected).

h	${\mathit{e}}_{\mathit{h}}^{{\mathit{L}}^1}$ UC	${\mathit{\alpha }}_{\mathit{i}}$ UC	${\mathit{e}}_{\mathit{h}}^{{\mathit{L}}^1}$ FOC	${\mathit{\alpha }}_{\mathit{i}}$ FOC	${\mathit{e}}_{\mathit{h}}^{{\mathit{L}}^1}$ SOC	${\mathit{\alpha }}_{\mathit{i}}$ SOC
70.9514	26.6976	-	7.2053	-	5.3999	-
47.8249	15.6275	1.3577	3.7519	1.6543	2.5007	1.9517
36.1943	10.7943	1.3279	2.3784	1.6358	1.4547	1.9444
29.2541	8.1559	1.3166	1.6669	1.6700	0.9448	2.0273
23.7852	6.5166	1.0843	1.2480	1.3982	0.6631	1.7104
20.9633	5.4088	1.4755	0.9782	1.9287	0.4903	2.3901
18.2621	4.6133	1.1531	0.7934	1.5179	0.3781	1.8835
16.1160	4.0165	1.1081	0.6596	1.4781	0.2992	1.8739
14.6520	3.5535	1.2862	0.5592	1.7327	0.2431	2.1808
12.2806	2.8821	1.1860	0.4213	1.6040	0.1698	2.0322
10.7534	2.4213	1.3119	0.3315	1.8053	0.1252	2.2939

, for which the order of convergence $\alpha$ has been calculated as

$${\alpha }_i=\frac{log(e_{h_i}/e_{h_{i+1}})}{log(h_i/h_{i+1})},$$

(52)

The error behavior is depicted in Figure 10, where the logarithmic scales were chosen strategically to represent the rate of convergence as the slope of the curves. These results prove the effectiveness of the corrective fluxes as expected, reducing the magnitude of the error and improving the rate of convergence with both first and second-order corrections. In each case, corrective fluxes recover the order of convergence previously lost due to the mesh non-orthogonality and skewness effects.

6.2. SPE Comparative Solution Project 2

The SPE second comparative solution project (SCSP) is a coning study of three phases in a single-well radial cross-section with logarithmic distribution in the radial axis. This simulation case consists of one production well having step variations in the production rate, representing a numerical challenge, which serves for testing convergence and stability capabilities [91].

The thermodynamic behavior of the fluid studied in the SCSP is represented by Black Oil fluid model (see

Table 3. PVT Table for SCSP [91].

Pressure (Psia)	Bo (RB/STB)	ρo (lbm/${\mathrm{ft}}^3$)	${\mathit{\mu }}_{\mathbf{o}}$ (cp)	Rs (scf/STB)	Bg (Mcf/STB)	ρg (lbm/${\mathit{ft}}^3$)	${\mathit{\mu }}_{\mathbf{g}}$ (cp)	Bw (Mcf/STB)	ρw (lbm/${\mathit{ft}}^3$)	${\mathit{\mu }}_{\mathbf{w}}$ (cp)
400	1.0120	46.497	1.17	165	5.90	2.119	0.0130	1.01303	62.212	0.96
800	1.0255	48.100	1.14	335	2.95	4.238	0.0135	1.01182	62.286	0.96
1200	1.0380	49.372	1.11	500	1.96	6.397	0.0140	1.01061	62.360	0.96
1600	1.0510	50.726	1.08	665	1.47	8.506	0.0145	1.00940	62.436	0.96
2000	1.0630	52.072	1.06	828	1.18	10.596	0.0150	1.00820	62.510	0.96
2400	1.0750	53.318	1.03	985	0.98	12.758	0.0155	1.00700	62.585	0.96
2800	1.0870	54.399	1.00	1130	0.84	14.885	0.0160	1.00580	62.659	0.96
3200	1.0985	55.424	0.98	1270	0.74	16.896	0.0165	1.00460	62.734	0.96
3600	1.1100	56.203	0.95	1390	0.65	19.236	0.0170	1.00341	62.808	0.96
4000	1.1200	56.930	0.94	1500	0.59	21.192	0.0175	1.00222	62.883	0.96
4400	1.1300	57.534	0.92	1600	0.54	23.154	0.0180	1.00103	62.958	0.96
4800	1.1400	57.864	0.91	1676	0.49	25.517	0.0185	0.99985	63.032	0.96
5200	1.1480	58.267	0.90	1750	0.45	27.785	0.0190	0.99866	63.107	0.96
5600	1.1550	58.564	0.89	1810	0.42	29.769	0.0195	0.99749	63.181	0.96

). Table 3 shows the densities $\left(\rho \right)$, volumetric factors $\left(B\right)$, viscosities ($\mu$), and gas-oil ratio $\left(R_s\right)$, as a function of the fluid pressure. Black Oil PVT data is converted into equilibrium relations, which serve for calculating thermodynamic state, i.e., components and phases composition, and gas and oil molar fractions in oleic and gaseous phases. FlowTraM uses this converted data for solving thermodynamic state, molar balance, and volume conservation equations in a pseudo-compositional model.

SCSP reservoir is initially at gravitational and capillary equilibrium conditions. Thus, fluids locate in specific zones, and pressure increases with depth. Figure 11a–c show gas, oil, and water saturation, respectively, with their transition zones. Additionally, pressure variation with depth is presented in Figure 11d.

We compared FlowTraM simulation results with the predictions of CMG IMEX, which we obtained by running a pre-designed input file in CMG library; and with Shell simulator results, which we took from SCSP paper results [91]. Borehole pressure, oil rate production, and water cut at the production well are shown in Figure 12, Figure 13 and Figure 14, respectively. It is worth noting that both Shell, IMEX, and FlowTraM results are in accordance. Nonetheless, Gas-Oil ratio (GOR) results on FlowTraM, presented in Figure 15, are significantly higher than IMEX results. These results are due to physical model differences. However, FlowTraM GOR results fit better with Shell results. Furthermore, FlowTraM results approximate to the average between other simulators GOR reported in SCSP paper results [91].

6.3. Colombian Volatile Oil Field Case

The field study case is a sector model of a volatile oil field in the Colombian foothills, Llanos basin. A heterogeneous simulation grid with 37,960 blocks was used to simulate the recovery process at a reservoir scale, with permeability changes between 1 and 100 mD and porosity between 1–6.7%. The mesh porosity and permeability distributions of an intermediate layer are shown in Figure 16 and Figure 17. The dynamic model was calibrated using a history matching process in FlowTraM. The vertical dashed lines represent production and injection wells.

Historical well information is the most sensitive information when making a history matching. The dynamic model must be able to reproduce the historical pressure or production data, depending on the operating condition imposed on it. That is, the simulator must reproduce the historical production if it is fed with the flowing well pressure, or reproduce the historical background pressure if it is fed with the production data. For the system evaluated in this work, the simulator is fed with background pressure data for producing wells (seeking to reproduce production) and with injection flow data for injection wells (seeking to reproduce background pressure).

For initializing the simulation model, an initial pressure of 4500 psi is set at 12,000 ft. Initial water saturation is 0.22, and the average oil saturation is 0.78. The distance between wells is about 11,000 ft. Producing well is programmed to produce with an increasing bottom hole pressure from 3000 psi at the beginning, to 3500 psi at day 1110, then the well continues producing at a constant pressure of 3500 psi. The injector well is programmed to start injection with 45 MMSCF/day from day 1110.

For obtaining the history matching, the information supplied by an operator company was directly introduced into the flow simulator. The production was evaluated, taking into account the historical productivity and sensitivity of the flow correlation in pipes, aiming to obtain well constraints that fitted its productivity.

Oil and gas production are presented respectively in Figure 18 and Figure 19, where it is shown that the model can reproduce the effects of recycle gas injection and the production history of the reservoir. At a field scale, the simulation results showed a good agreement with field measurements. Also, the gas saturation distribution is presented in Figure 20.

Once the historical matching has been verified, different gas injection is evaluated from day 1423: CO${}_2$, Flue gas, and recycle gas at a fixed rate of 12 MMSCF/day. Figure 21 and Figure 22 show the simulation results for oil and gas, respectively, with the different injection gases, where it is observed that CO${}_2$ is the gas with the best productivity response.

7. Conclusions and Future Work

This work serves as a basis for the construction of more robust numerical methods. The goal is to account for significant problems regarding flow in highly anisotropic and heterogeneous porous media, as in many applications in the oil industry, whose solutions benefit from conservative numerical methods. In this sense, a new version of the unstructured finite volume method is proposed, with an improved scheme of the traditional two-point flux approximation, allowing for dealing with the phenomenon of non-orthogonality of the meshes and the effect of the full permeability tensor. Additionally, a direct consequence of our formulation is the correction of errors due to the mesh skewness. Moreover, the mathematical formulation based on the second-order anisotropic elliptic operator allows the method to be applied in diverse engineering applications, like enhanced oil recovery, at a low cost. Furthermore, we proposed a modification to Newton’s method to efficiently apply the corrections, maintaining the sparse nature of the Jacobian matrix.

The comparisons done in several numerical examples, between FlowTraM, CMG-IMEX, and Shell simulator, demonstrate that our results follow other commercial simulators’. In the specific case of SCSP, our approximation to Gas-Oil ratio is closer to the average between GOR values reported by other simulators in SCSP paper results [91]. Moreover, on a Colombian field history-matching, we present the capability of FlowTraM for reproducing an actual field operation.

The next stage of this work is to analyze and to compare the performance of the UFVM here presented with respect to the variational methods as DG and FEM. The current development and implementation of UFVM are planned to be extended to chemical transport and thermal recovery applications. Also, more elaborate models can be included for unconventional wells [92] and elements of any geometries [93] in the source/sink term discrete equations.