## 1 Introduction
\[\mathbf{v}_{a}=-\lambda_{a}\mathbf{K}(\
abla p_{a}-\rho_{a} \mathbf{g}) \tag{1}\] \[\frac{\partial}{\partial t}(\phi A_{a})+\
abla\mathbf{u}_{a}+q_{ a}=0 \tag{2}\]  
where the subscript \(a\) denotes the pseudocomponent oil (o), gas (g), or water (w). For each one, the variables are defined as follows: \(\mathbf{v}_{a}\) is the phase flux, and \(\lambda_{a}\) is the mobility of the phase. This is defined as \(\lambda_{a}=\frac{k_{r,a}}{\mu_{a}}\), where \(\mu_{a}\) and \(k_{r,a}\) are the viscosity and the relative permeability of each phase, respectively. In short, relative permeability is a measure of the effective permeability reduction of a phase in the presence of other phases and is a function of saturation. \(p_{a}\) is the pressure, \(\rho_{a}\) is the density, and \(q_{a}\) is the well outflux/influx of each phase, i.e., the volumetric flow rate at which each pseudocomponent is produced/injected from the system. \(A_{a}=f(b_{w},b_{w},b_{g},s_{w},s_{w},s_{g})\) are accumulation terms dependent on phase saturation (\(s_{a}\)) and shrinkage phase factors (\(b_{a}\)). The flux terms \(\mathbf{u}_{a}=f(\mathbf{v}_{w},\mathbf{v}_{o},\mathbf{v}_{g},b_{w},b_{o},b_{g})\) are known functions of shrinkage factor and phase fluxes. Finally, the terms \(\phi\) and \(\mathbf{K}\) are the porosity and permeability, respectively. The other popular simulation model is the compositional one [18; 19]. The difference between the black-oil and compositional models lies in that the former traces the flow of each separate phase; hence, it only considers single-phase (gas, oil, and water) density and compressibility effects, as well as dissolution of gas in the oil phase. The fact that the volumetric and dissolution factors can be easily described by simple functions of the pressure is based on the assumption that the overall composition does not change significantly during reservoir depletion or waterflooding. On the other hand, compositional modeling, facilitated by means of an EoS model, is used when some \"alien\" composition is incorporated in the reservoir fluids. This is the typical case when handling EOR cases of CO\({}_{2}\), natural or acid gas injection, and injection of dry gas in gas condensate reservoirs. Equations (1) and (2) are combined in a second-order partial differential equation (PDE) which can be analytically solved under simplifying assumptions (incompressible, single-phase system with constant porosity and permeability). However, when trying to model realistic subsurface formations, these assumptions do not hold. In this case, the discretization of the equation through linearization is a key step in numerical reservoir simulation. This process involves dividing the reservoir into a grid to represent the spatial distribution of rock properties and fluid flow and using a linearization scheme [20]. Common methods for discretization include finite differences, finite volumes, and finite elements. Finite volumes [21], in particular, are widely used in reservoir simulation due to their simplicity and efficiency. Gridding plays a vital role in the discretization process, involving defining the size and shape of the cells within the reservoir grid. Various types of grids, such as Cartesian, corner-point, and unstructured grids, are used based on the geological complexity of the reservoir. The choice of grid impacts the accuracy and computational efficiency of the simulation. Petrophysical and fluid properties, such as porosity, permeability, fluid saturation, compressibility, and relative permeability, are crucial inputs for reservoir simulation [22; 23]. Accurate characterization of these properties is essential for realistic simulation results. Furthermore, incorporating phase behavior models is essential to accurately representing the complex interactions among different phases of fluids and dissolved particles in multiphase flow simulations.