Methods in Oil Recovery Processes and Reservoir Simulation

The exploitation of an oil field is a complex and multidisciplinary task, which demands a lot of prior knowledge, time, and money. A good reservoir characterization is deemed essential in the accomplishment of Enhanced Oil Recovery (EOR) processes in order to estimate accurately the properties of the porous medium affecting the flow properties. Several techniques at a field scale are currently being used to determine these properties, which are time and money consuming. But these alone do not guarantee the success of the project. Reservoir simulation and numerical techniques were then included in the pre-development and follow-up studies as an effective tool to determine the productivity and future behavior of the oil field. As the computational power increased, more advanced and detailed models were developed, including different chemical and physical phenomena. But alongside this process, there was an active research in the area of reservoir simulation, improving the accuracy and efficiency of the numerical schemes used for the flow, transport, and energy equations. The aim of this review is to address the topics described. Firstly, the origin of an oil recovery process, the economic factors and field tests involved are introduced. Secondly, the oil and porous medium origin and characterization as well as an introduction to the fundamental concepts and equations are associated to reservoir simulation. Finally, a brief description and analysis of the techniques are used in reservoir simulation employing finite difference methods, their downsides and possible ways to overcome these problems.


Introduction
The era of discovery and subsequent exploitation of the denominated "easy oil" (part of the so called conventional reserves) is over [1]- [5]. After primary and secondary recovery, oil companies have begun with more complex Enhanced Oil Recovery (EOR) processes involving fluids or chemical reactions that require a greater detail of the phenomena taking place in the porous medium. In addition to these complications there exist non-technical related factors that might increase the risk of an investment [6]. An example of this may be the development of off-shore platforms, plants located in remote areas of the world or the exploitation of unconventional deposits (heavy oil, shale oil, tar sands). All these projects require a prior investment (in the order of several tens of millions of dollars high risk projects), so companies make use of several techniques in the stage of feasibility analysis to limit this risk and increase the chances of success of the operation (Figure 1) [7] [8].
These techniques, used to predict and optimize the exploitation, include laboratory and field tests (e.g. seismic 2D/3D/4C before starting operations and 4D to follow up the changes during exploitation, geostatistics) to give an idea of the conditions of the porous medium and its production performance [9]- [13]. However, these tools have proven to be insufficient [12] [14] [15] and therefore oil companies have begun using computational tools to predict and optimize their projects and production facilities, which is known as Reservoir Simulation. The latter consists in solving numerically the differential equations describing the fluid flow in porous media, which have no analytical solution, taking into account all geological, physical, chemical and/or mechanical Figure 1. A typical cash-flow of an exploration and production project (adapted from [8]).
phenomena occurring during operation so as to analyze and predict behavior as a function of time. The reservoir simulation can be used as well in inverse engineering problems for optimizing existing numerical models and couple the dynamic/historic data (production) in the simulation [15]- [18].
Generally speaking, reservoir simulation consists of three main parts: the physical characterization of a geological model describing the rock formation; a model characterizing the fluid flow and finally "well models" which describe the conditions under which fluids are injected or extracted from the reservoir [19]. The latter, along with the wellbore and the primary surface facilities, constitutes what is known as upstream.
During the last 30 years numerous theoretical and practical advances have been developed due to the appearance of new numerical techniques and increased computational power, respectively [15] [20]. This led to a new generation of more complex and detailed models. The accurate representation of the reservoir and the fluid contained in it is an issue that still needs to be more carefully addressed in order to reduce risks in exploration and production (E&P) projects. One of the most important points is related to the scales of the models: current grids are still large when compared to the geological characterization, description of chemical processes or fluid flow. Most of current numerical models used in reservoir simulation may contain from 10 5 to 10 8 grid blocks, depending on the model type, complexity, computational power available and fluid behavior.
The development of increasingly complex and detailed models requires the use of numerical techniques to solve these at reasonable times. Moreover, the representation of the properties of the porous medium and the characteristics of crude oil and natural gas at high pressures and temperatures may differ from laboratory tests, causing differences between the results and simulation. Another important topic is how to assess and properly estimate the properties of the rock formation. Geologists use statistical techniques in order to recreate the model properties of a porous medium, which are determined by several tests (e.g. seismic studies, drill core samples and even production data). However, other numerical tools are required (e.g. Monte-Carlo or Stochastic Processes) to take into account the effects of uncertainties in the model [21] [22].
The aim of this review is to briefly describe the characterization of porous media and its main parameters. Then, models governing fluid flow in porous media are introduced. This consists in Darcy's equation and the more complex compositional model for multiphase flow (where the component transfer between phases introduces nonlinearities in the mass transport equation), which is useful to describe chemical EOR operations. The second part comprises a brief explanation of the numerical methods used as well as the problems associated with these schemes (numerical dispersion and dissipation), and possible solutions using advanced numerical schemes.

Origin of the Oil
A reservoir is an underground trap where different fluids (water, oil and gas) have accumulated due to a migration from the source rock where they were originated. The porous medium is generally considered of sedimentary origin and consists of a series of microchannels (about 1 -100 microns diameter) interconnected where these fluids can flow ( Figure 2). These formations may have some tens of meters thickness, but may extend several kilometers in the lateral directions [23].
The origin of crude oil prior to the migration and deposit of hydrocarbons in the porous media is also a long and complex process [25]- [27]. The source of hydrocarbons consists of a series of phenomena, both organic and inorganic, taking place during long periods of time (in the order of million years) [28]. Most of these hydrocarbons originate in organic decomposition processes. The first stage (called diagenesis) involves the sedimentation of remains of dead plants and animals. Under these conditions, at low depths, the action of bacteria produce methane, water and CO 2 , leaving as reaction products of kerogen substances (cyclic and large hydrocarbon molecules containing oxygen, nitrogen and sulfur). The second process, called catagenesis, occurs at greater depths and temperatures. In this region, known as the oil window zone, the kerogen molecules break in smaller, heavy hydrocarbons forming the oil phase. At higher temperatures the cracking of hydrocarbons continues creating lighter compounds, first wet gas and subsequently dry gas. Finally, at even greater depths, the last process takes place (called metagenesis), where the remaining kerogen is out of hydrocarbons and subsequent cracking processes terminate when there is no hydrogen in the compound. The result of these reactions is the formation of graphite ( Figure 3).

Porosity and Permeability
A porous rock formation is composed of a solid part, called solid matrix, and the remaining void space or microchannels whereto oil migrates [31]. The volumetric fraction of these channels is denominated porosity of the porous media. The latter depends on the fluid pressure, if the rock is compressible, or in some other phenomena which may take place (e.g. adsorption of chemical components during EOR processes). The following list provides typical values of porosity according to the origin of the rock formation [32] [33]: Consolidated sandstones 0.05 to 0.3; limestones 0.1 to 0.4; uniform spheres with minimal porosity packing 0.25; uniform spheres with normal packing  The permeability of the formation is a property that characterizes the ease with which fluids can flow when a pressure gradient is applied between two points. Nonetheless, reservoir rocks usually have no uniformity in their properties because of the mechanisms involved on its formation, thus the permeability will have a large dispersity in its values [34] [35]. The fluids used commonly in EOR operations, more mobile than oil, occupy the high-permeability zones (e.g. faults or fractures), with the result that large areas of oil will be bypassed, reducing the efficiency of the process. Mathematically, the permeability can be expressed as a diagonal tensor (K). When the medium is isotropic then the permeability can be represented as a scalar function. Due to transitions between different rock types, the permeability may vary swiftly throughout the reservoir, going from extreme low permeability areas of 1 mD to areas with permeabilities exceeding 10 D.
In order to develop a mathematical model of the porous medium, Corey [32] estab-lished several restrictions: the whole void space of the porous medium is interconnected; the mean free path length of the fluid molecules or molecules contained in the fluid must be negligible when compared to the dimensions of the pore channels, and the dimensions of the void space must be small enough so that the fluid flow is controlled by adhesive forces at fluid-solid interfaces and cohesive forces at fluid-fluid interfaces (in multiphase systems) [31]. These assumptions allow excluding any disconnected channels in which there can be no fluid flow, eliminating the difference between the concepts of total porosity and effective porosity. Furthermore, since the dimensions of the molecules or particles in the fluid are negligible with respect to the microchannels, a suitable replacement for a hypothetical continuous medium can be performed.

Representative Elementary Volume
The flow of reservoir fluids in porous media can be described at several different scales, from a microscopic to a macroscopic/formation scale. In order to perform large-scale reservoir simulations, a microscopic description of the flow channels would be too demanding for the computational power available and besides to characterize a reservoir rock so accurately to determine the geometry of the pore network is beyond the scope of modern techniques and equipment. A continuum scale description is then utilized, and its behavior is governed by forces acting between the different fluids and the rock formation. The goal of a reservoir continuum model is then to average both the fluids and reservoir rock [28] [36]- [39]. In order to develop the mathematical model based on a continuum, the concept of Representative Elementary Volume (REV) ( Figure 4) is introduced. This is based on the hypothesis that certain properties of both the fluid and the rock may be considered constant along a certain range of scale and thus it establishes limits for the physical scales in the numerical models. If a REV cannot be identified for a specific porous medium then this concept cannot be applied and the macroscopic approach should be discarded [40].
The procedure for estimating REV dimensions and establish boundaries between microscopic and macroscopic scales is explained below using the porosity in Figure 4 as an example [31]. A porous medium is then considered occupying the domain Ω, with a volume V(Ω). A subdomain Ω 0 (d) ⊂ Ω with a characteristic dimension d is also defined. Furthermore, the porosity piece-wise function is defined in Equation (1) as follows, (1) Figure 4. Scheme showing the boundaries to determine the REV (adapted from [19]).
Then, the porosity of an element with characteristic dimension d is defined by Equa-

Fluid Flow Models
In underground processes in porous media, fluid flow involves mainly convection (advection and diffusion) of the different phases through a heterogeneous medium.
The equations used to describe the flow at a microscopic level or poralscale are variations of Navier-Stokes (creeping flow) and the mass conservation law. At a macroscopic level, Darcy's law [43] was derived and it is used to describe flow behavior. Also, the effects of the displacing process on the rock formation will be considered negligible even though these mechanisms are sometimes necessary to represent first-order effects (e.g. adsorption) [44] [45]. When dealing with multiphase/ multicomponent flows, some of the processes therein involved are characterized by the chemical and physical interaction among the components present in the fluids.
Therefore, diffusive and/or dispersive mixing of these components is often critical and must be correctly understood and modeled in order to obtain accurate results.
Molecular diffusion is typically quite small in porous media processes. Nevertheless, hydrodynamic dispersion may be important and therefore it should be incorporated in the flow equations. This can be done by means of the standard diffusion/dispersion tensor [46] [47], provided that the degree of the heterogeneity is not too large (Dykstra-Parsons coefficient < 0.50) resulting in a Fickian behavior.
Generally, three fluid phases may exist inside a reservoir ( Figure 2): aqueous/ brine (the phase containing predominantly water), oleous (the phase containing liquid hydrocarbons) and gas phase, which contains lighter gaseous hydrocarbons ( Figure 3). In the case of single phase systems, the void space of the porous medium is filled by a single fluid or by several fluids completely miscible with each other. In multiphase systems the void space is filled by two or more fluids that are immiscible with each other, thus maintaining a distinct boundary between them. Formally speaking, the solid matrix of the porous medium can also be considered as a phase called the solid phase. Each phase may also be composed by many chemical components. For example, oil is a very complex mixture of hundreds of hydrocarbons with different chemical properties. In order to derive a set of equations for the fluid flow some terms should be defined beforehand. Firstly, the term phase stands for matter that has a homogeneous chemical composition and physical state of a system under consideration that is separated from other such portions by a definite physical boundary. Secondly, it is defined as component present in a phase to the matter that is composed of an identifiable homogeneous unique chemical species or of an assembly of species [41]. The number of components needed to describe a phase is given by the conceptual model, i.e. it depends on the physical processes to be modeled and the desired accuracy. In many oil reservoirs (above bubble point pressure) crude oil contains some amount of dissolved gas. It is acceptable to assume that the oil and gas compositions are fixed [15] [20] [26], and that the solubility of the gas in the oil depends on pressure only. Provided these conditions are met, then it is possible to consider the pseudo-components oil and gas.
Both microscopic and macroscopic effects control the movement of fluids in the reservoir. At the pore scale, interfacial tension (IFT) and capillary effects control the fluid behavior. Macroscopically, fluid flow is controlled by reservoir heterogeneity and mobility differences between the fluids. Viscosities, capillary pressure, IFT and mobility differences vary throughout the reservoir and depend mainly on phase saturations, their interactions and molecular compositions. Chemical components may transfer between contacting phases, altering the fluid properties of both. Interactions between the fluids or their components and the reservoir rock may also impact performance (e.g. adsorption of chemical components onto the surface of the rock altering the wettability). Thermal effects are generally very small due to the large heat capacity of the rock. However, in EOR thermal processes (steam flooding or in-situ combustion), the conservation of energy in the REV should be considered.

Single Phase Flow
The governing equations for single phase flow in porous media are the conservation of mass, the Darcy equation and an equation of state (EOS). Considering the flow of a single fluid with density ρ through a REV of a porous medium the differential form of the continuity Equation (3) where φ is the porosity of the rock formation, q represents the fluid source/sink term and u is the velocity vector. The fundamental law of fluid flow in a porous medium is Darcy's law [43]. It gives the effective flow velocity across a REV of porous media and thus does not analyze the flow at a microscopic scale. In its differential form, the relationship between velocity and pressure drop is given by Equ- where K is the absolute permeability tensor of the porous medium, µ the fluid viscosity, g the acceleration field, and z represents physical dimensions. In most of the cases, it is possible to assume that K is a diagonal tensor as presented in When 11 22 33 k k k = = , the porous medium is called isotropic; otherwise, it is anisotropic. Generally in porous media can be considered that both lateral permeabilities are in the same order of magnitude while the vertical permeability is considerably lower (at least one order of magnitude) than the other two components. Originally, Darcy's law was derived experimentally and was thus considered an empirical law [28]. Several authors reported the derivation of Darcy's law based on volume averaging of the Navier-Stokes momentum equations [38] [41] [46] [56]- [61]. The assumptions needed for the derivation of Darcy's law include low flow speeds, Newtonian behavior and that the pore/fluid friction is the dominating force acting on the fluid. Also, the porous medium is assumed to be rigid and not compacted due to fluid flow. Introducing rock and fluid compressibilities in Equation (6), c r and c f respectively, both equations can be coupled in the parabolic Equation (7) for the fluid pressure, In the special case of incompressible rock and fluid (generally acceptable for liquid systems) the partial differential equation (PDE) simplifies to a Poisson elliptic equation.

Two Phase Flow
The space in a reservoir is generally filled by both an oleous phase and brine. In addition, during secondary recovery processes, water is frequently injected in order to improve oil recovery. If the fluids are immiscible, they are referred to as phases. A two-phase system is commonly divided into a wetting and a non-wetting phase, given by the contact angle between the solid surface and the fluid-fluid interface on the microscale. For each pair of phases, one phase will wet the rock more than the other phase, and that phase will be referred to as the wetting phase (j = w). The other phase is then the non-wetting phase (j = nw). Normally, water is the wetting phase in a water-oil system, and oil is the wetting phase in an oil-gas system. In the absence of phase transitions, the saturations change when one phase displaces the other. During the displacement, the ability of one phase to move is affected by the interaction with the other phase at the pore scale. In the mathematical model, at a macroscopic scale, this effect is represented by the relative permeabilities k rj (j = w, nw), which are a dimensionless scaling factor that depends on the saturation and modifies the absolute permeability to account for the rock's reduced capability to make one phase to flow in the presence of the other. Then, the mass conservation Equation (8) for each phase yields [15], And the multiphase extension of Darcy's law is presented in Equation (9), Together, they form the basic system of equations. Because of the interfacial tension (IFT), the pressure in the two phases will differ. This difference is called capil- , which is usually assumed on the macroscale to be a function of saturation. From the formulation exposed previously, the following Equations (10) and (11) can be derived, where T u is the total Darcy velocity, which is useful in schemes employed to solve the system of equations, and ( ) w w f S is the fractional flow of the wetting phase.
The system of equations derived from the formulation of two phases can be solved using various numerical techniques. The most used schemes are: formulation using the pressure of both phases, known as simultaneous solutions (SS); formulation using pressure and saturation phase, known as IMPES or IMPSAT, depending on the numerical treatment of the saturation equation (explicit or implicit, respectively), and the global pressure formulation [62].
The volume fraction occupied by each phase is defined as the saturation of that phase. Thus, for a two-phase system, and considering no phase transitions, the sum of the saturation of both the wetting and non-wetting phases is equal to unity, as presented in Equation (12). Similar to the void space indicator function, the phase indicator piece-wise function is defined by Equation (13), Then, the saturation of the phase j in an REV Ω 0 element with characteristic dimension x 0 will be defined by Equation (14), The relative permeability of each phase depends on the phase saturations but does not depend directly on fluid flow properties [63]. If only a single phase is present the relative permeability has no physical meaning, but in a multiphase system, the flow of one phase interferes with the others, hence this influence is taken into account in the Darcy equation (k rj ≤ 1). It is usual in multiphase systems to use correlations of the relative permeabilities expressed as functions of the wetting phase saturation (Equation (15) and Figure 5), In addition to relative permeability correlations, also analytical capillary pressure functions are needed. In two phase simulations it is standard to use the relations provided by either Brooks-Corey or Van Genuchten [36]. As for the relative permeability, these depend on empirical constants (e.g., if the system is oil-wet or water-wet), so several models have been developed through the years [64]- [70].

Compositional Models
In two phase models it was assumed a no mass transfer condition between the phases. This assumption is valid for two phase flows of water/brine and oil, which is often the case in primary and secondary recovery mechanisms. In EOR processes, mass transfer and compositional effects are deemed essential to model accurately as they may become the driving mechanisms for the displacing process. A typical reservoir fluid may consist of several different chemical pseudo-components. Fully compositional models must be used when the fluid flow depends strongly on component transfer between phases. In fact, many EOR techniques, mainly chemical and miscible gas processes, are specifically designed to take advantage of the phase behavior of multicomponent fluid systems. Because these components may be transferred between phases (and change their composition), the basic conservation laws must be expressed for each component instead for each phase. For a chemical flooding compositional model, the governing differential equations consist of a mass conservation equation for each component, Equation (16), and Darcy's law for each phase [20] [71] [72].
Here z i is the overall concentration of component i calculated by Equation (17) [73]).
In addition to the advective, directional movement of a component described by the Darcy phase velocity, components may also move due to dispersive forces. The simplest movement is molecular diffusion described by the random Brownian motion of molecules. Such motion is usually considered in reservoir simulation of negligible importance compared to other forces acting on the fluid. A more substantial phenomenon is mechanical dispersion. Narrow channel flows experience parabolic diffusion along the fronts (Taylor dispersion) and the irregular pore networks disperse the mass at a microscale ( Figure 6). The tensor of hydrodynamic dispersion considering the mentioned effects is expressed in Equation (18) The dispersion part of the tensor is significantly larger than the molecular diffusion; also, dl j is usually considerably larger than dt j and their relationship can be expressed as a function of the Peclet number [20] [46]. Nevertheless, at low Peclet numbers (Pe < 10) both dispersitivies are of the same order of magnitude [75]- [77].
Dispersion can represent small scale movements not captured by the REV used in the mathematical model, but according to Heimsund [28], taking this into account in the numerical model may be troublesome [78] [79]. It is worth mentioning that some numerical schemes, especially first order methods (see point 3.2), add artificial diffusion which is most of the times far greater than the physical dispersion discussed here. It is then advisable that when using naturally dissipative schemes the physical dispersion to be neglected [28] [54] [80].
If a compositional model is formed by a system of N comp components and N p phases, there are a total of ( ) , , , , .

Energy Equation
In case the flow cannot be considered isothermal, or the recovery process involves the addition of considerable amounts of energy to the reservoir, an extra condition and variable must be introduced to the system. The conservation of energy, Equation (22), and its dependent variable, the temperature, are then added to the system. The major difference with respect to the other equations is that the energy is also conducted by the rock formation, and not only between the phases. If the local thermal equilibrium concept is applied, the temperature in the REV for all the phases and the porousmedium is considered to be the same and the energy equation is as follows [20], where U j is the specific internal energy, C s the specific heat capacity of the rock, H j the enthalpy of the phase, t k the thermal conductivity tensor, q H the enthalpy source term, and q L the heat loss. In most of the chemical EOR operations, the heat transfer is considered negligible and therefore an isothermal assumption is valid.

Well Models
A production/injection well is a vertical (or vertical/horizontal in case of horizontal wells), open hole through which fluid can flow in and out of the reservoir, according to the strategies or its degree of maturity. These are cemented and then perforated along specific intervals (multi-zone wells). The primary function of production wells is to extract hydrocarbons and later on, the water/chemical products injected as part of EOR processes. On the other hand, injection wells can be used for disposal of certain fluid (e.g. CO 2 storage) as well as to inject chemical solutions so as to increase the recovery efficiency, sweeping the oil towards production wells.
where s is the skin factor, r w is the well radius and r 0 is the effective block radius at which the steady-state pressure equals the computed block pressure. The Peaceman model has been also extended to horizontal wells and it was modified to take into account non-square grids, boundary blocks and non-Darcy effects [45] [86]- [89].

Numerical Techniques for Fluid Flow in Porous Media
Reservoir flow problems can be highly complex, consisting of many different physical effects when it comes to EOR processes. The analysis of all these phenomena can be achieved, up to some extent, by laboratory experiments or field tests at small scale, but these tend to be expensive to conduct may not be extrapolated to the whole reservoir. In order to solve this problem, mathematical models became progressively more important. Using these along with analytical solutions, engineers provided basic performance predictions so as to modify production strategies.
Several numerical formulations are employed to solve the non-linear systems of equations. The most stable approach is a fully implicit solution technique in pressure and saturation/concentration, but this generally leads to large, ill-conditioned matrices. Another scheme broadly utilized in compositional formulation consists, in order to reduce the level of implicitness, to solve the pressure equation system implicitly (which can be viewed as an overall volume balance) plus a sequence of (N comp − 1) components conservation equations [20] [90]- [95]. The conservation equations have a strongly hyperbolic character due to the advective term. The overall technique is called an IMPEC (implicit in pressure, explicit in concentration). IMPEC methods are often limited by stability restrictions on the time step size due to their explicit scheme, but solutions do not suffer less smoothing than fully implicit methods, which strongly affect the performance prediction of compositional models [72]. Different implicit techniques are also used and often provide enhanced efficiency, allowing bigger time steps than those employed with IMPEC. One of these methods is the implicit pressure and saturation (IMPSAT) procedure, in which pressures and (N p − 1) saturations (but not compositions) are determined implicitly [82]. Different techniques worth of mentioning are: Adaptive implicit (AIM) [82] [96], bilinear approximation techniques [97], preconditioning schemes [98] [99], parallel computing and adaptive mesh refinement (AMR) in compositional simulation [72] [100] [101].

Numerical Schemes
The aim of this section is the derivation and explanation of the numerical schemes to be used for solving differential equations presented above, as well as also explain the reasons for the occurrence and possible numerical solutions of certain phenomena that affect simulation results [102]. The equation to be used as a model is the one-dimensional advection-diffusion Equation (24), which is a simplification of the continuity equation presented for the compositional model, Using a finite-difference approach the continuous domain is transformed into a discrete representation with a finite number of points in both a spatial (i) and temporal ( n ) grids. Then, the time derivative in previous equation is expressed using a Taylor series expansion around the point n i u yielding Equations (25) and (26), is the Bachmann-Landau notation used to describe the extra error terms in the Taylor approximation to the time derivative. Time and spatial operators may have finite-difference schemes with different orders of accuracy and in this case the overall order of the equation is determined by the differential operator with the largest truncation error. Noteworthy is also that while the truncation error is expressed for the differential operator, the numerical algorithms will not be expressed in terms of the differential operators and will therefore have different (usually smaller) truncation errors. Following a similar procedure, the finite-difference approximations in Equations (27) and (28) are obtained for the space derivative in backwards and centered forms as, Even though the centered scheme has a higher order of precision, it generates unstable results when applied to the advection equation (wave transport equation). Hence, numerical methods using only points that are "upwind" of the wavefront are employed [103].
The inclusion of diffusion phenomena in the description of a fluid flow leads to non-trivial complications in the numerical solution of the mass conservation equations. From an analytical point of view, the resulting equations are no longer purely hyperbolic PDE's but rather mixed hyperbolic-parabolic PDE's. This means that the numerical method used to solve them must necessarily be able to cope with the parabolic part of the equations. For the diffusive term, the second order derivative is usually discretized using a centered scheme yielding Equation (29) In Equation (30) the order of accuracy is expressed as function of both independent variables ( ) , inferring that both discretization schemes for spatial and temporal grids have influence in the total error of the numerical model. The solution at the new time-step 1 n + can be calculated explicitly from the quantities that are already known at the previous step n . This differs with an implicit scheme in which the finite-difference representations of the differential equation are expressed of terms at the new time-level 1 n + . These methods require solving a number of coupled algebraic equations. Table 1 presents several finite-difference schemes indicating their orders of accuracy both in spatial and temporal grids and their representation in a purely , , x t t x ∆ ∆ ∆ ∆   Table 2 summarizes the most common discretization schemes and their orders of accuracy [102]. Higher orders schemes allow increasing the accuracy of the numerical model but at the cost of requiring a higher number of points to make the evaluation of the derivatives. This increases both the computational cost and the difficulty of evaluating derivatives near the boundaries.
A key factor in all numerical schemes is the issue of how treating the solution on the boundaries of the spatial grid as the time evolution proceeds. Two types of con-

Numerical Dissipation and Dispersion
The exact solution of the discretized equation satisfies a PDE different from the one being solved. This difference is represented by the local truncation error (LTE) of The procedure to calculate this error and assess its contribution to the numeric solution is straight-forward. It consists in performing an expansion in a double Taylor series around a single point n i u , both in spatial and temporal grids to obtain a modified PDE. Besides, the high-order time derivatives as well as mixed derivatives must be transformed in terms of space derivatives using this modified PDE. The analysis for the truncation error in the 1D advective-diffusive equation is presented using two numerical schemes: firstly the upwind scheme and subsequently the Lax-Wendroff scheme. For the upwind scheme it yields Equation (32), ( ) Introducing these terms from Equations (33), (34) and (35) Rearranging Equation (36) Now the temporal derivatives in Equation (37) The numerical scheme does not solve the original PDE, but a modified PDE with extra terms of higher order derivatives. The extra term containing the second order derivative is interpreted as a numerical diffusion, additional to the physical coefficient D. As long as the Cr < 1 condition is met, the numerical solution will produce an artificial smearing given by the term ( ) One way to reduce these numerical errors is by means of additional, artificial factors which stabilize or decrease the previously seen effects. As an example, the streamline diffusion method consists in adding a term of artificial diffusion to counteract the added terms by the numerical scheme; the non-oscillatory shockcapturing methods, TVD (Total Variation Diminishing) or flux-limiters [19] can also be listed as improved schemes to overcome these effects. To reduce the influence of these differences a possible solution is also to use schemes of superior orders (see Table 2). As an example of these techniques, the same analysis done for the upwind scheme is performed using the Lax-Wendroff method. The latter due to the fact that the time derivative is expressed as a second order Taylor series expansion [104] [105]. The numerical model to solve is presented in Equation (40), Following a procedure similar to the previous scheme it renders Equation (41), As shown, the first diffusive term in the LTE has disappeared leaving the dispersive term as the main source of error in this method. Since this term is negative, the spurious oscillations occur behind steep fronts. It is worth mentioning that a more accurate numerical scheme is not necessarily a preferable one. As an example, the upwind and the Lax-Friedrichs methods are both dissipative, thelatterisgenericallymoredissipativedespitebeingofhigherorderaccuracyinspace.

Flux Limiters
In the previous section two of the most common numerical schemes utilized were introduced, and the advantages or disadvantages of each were studied and inferred.
While the upwind scheme can handle steep gradients, it is very diffusive and moreover a first order scheme; on the other hand, Lax-Wendroff is a second order scheme, less diffusive but presents serious problems when sharp gradients are present in the system. Therefore, new numerical schemes were published coupling low-and high-resolution methods, taking advantage of the mentioned characteristics [104] [106]- [115]. For this analysis, it is considered the 1D advection-diffusion Equation (42) in term of fluxes, with no diffusive terms.
( ) The idea behind this concept is then to write the fluxes as a function of low-and high-resolution numerical schemes as in Equation (43), using a proportionality factor.
( ) ( ) The proportionality factor, also called flux limiter function introduced in Equation (44), depends on the ratio of consecutive gradients in the numerical mesh, this is, Using the FTUS and Lax-Wendroff in Equation (45) as low-and high-resolution schemes respectively, the flux is calculated according to Equation (46), Finally, the discretized advection Equation (47) Equation (48) This is valid for positives values of r, when the following two conditions in In Equation (50) are met ( Figure 9).
Further, more restrictive constraints are applied in order to make the scheme second order in accuracy ( Figure 10). Several high-order flux limiter functions were developed within this region. These are characterized by a low numerical dispersion in high gradient fields as well as being less diffusive than traditional schemes (FTUS) in low gradient advection phenomena (Table 3 and Figure 11).

Consistency and Stability
This review concludes with the study of three concepts related with numerical simulation. The first to be defined is the consistency of a numerical scheme: Given a PDE in its operator form ( )    This concept may seem trivial, but numerical errors introduced during simulation Figure 11. Second-order TVD functions presented in Table 3.
Van Leer [122] can grow boundless, amplifying the errors until the system eventually collapses. This is ensured introducing the concept of numerical stability. To understand this, a system of equations written in the form of modified PDE Equation (31) where G is the amplification matrix of the numerical perturbations in the system. Using Equation (51), a relationship between the perturbations at the time step n with the initial perturbation can be established.
Equation (52) These three concepts can be summarized in the following theorem: Given a properly posed initial-value problem and a finite difference approximation to it, that satisfies the consistency criterion, stability is the necessary and sufficient condition for convergence. This theorem, known as the "Lax-Richtmyer equivalence theorem" or Fundamental Theorem of Numerical Analysis (Figure 12), is quite important since it shows that consistency, stability and convergence are strictly Figure 12. Scheme of the Lax-Richtmyer equivalence theorem.
related. In general, proving that the numerical scheme adopted is stable will validate that the discretized equations represent the PDE as well as the numerical errors in the simulation are bounded at all times [102].

Conclusions
Reservoir simulation is a branch of engineering that emerged in recent years since oil companies needed to justify and evaluate E&P investments. The former is not related to only one discipline, but includes the assistance and collaboration of various specialists so as to characterize a reservoir, estimate its profitability and give the "green light" to the project development phase. Oil reservoirs are geological traps where oil migrated and remained for long periods of time. An accurate determination of their physical characteristics has not yet been developed and statistical tools are used along with complex field tests to extrapolate these properties. In addition, oil is not a homogeneous, pure fluid but is composed by a large group of different components which may alter its properties. Hence, a set of variables must be previously evaluated and studied in order to get feasible results. The research and development of new exploration technologies to reduce the model uncertainties are deemed essential. These, along with production studies at reservoir scale on pilot wells will allow performing accurate history matching analysis. Due to the current oil reserve conditions and future production estimates, these new technologies should also consider their applicability also in non-conventional oil reservoirs (e.g.
oil sands, tight oil, shale oil), or in geographical areas with harsh conditions (e.g. off-shore platforms).
Fluid flow simulations are performed using two different approaches: the Navier-Stokes equations throughout a complex network of microchannels in the porous medium; or the assumption of a continuum with averaged properties using Darcy's equation, rendering a system independent of the geometry at a microscopic scale.
The first is only circumscribed to specific laboratory tests and has limited application in field studies. This is due to several factors, among them the uncertainties associated with the poral geometry in the field as well as high computational costs required to solve the system of equations. However, this approach might be useful in the design of new chemicals while being evaluated at a microscopic level. These studies should then be supplemented with scale reservoir simulations and field tests using Darcy's equation.
Subsequently, the mathematical tools for reservoir simulation using Finite Difference Methods (FDM) were presented. The errors introduced by these schemes as well as possible solutions to tackle these problems have been addressed. The numerical convergence of a system of PDE's is a critical aspect that must be taken into