Simulation of Surfactant Oil Recovery Processes and the Role of Phase Behaviour Parameters

Chemical Enhanced Oil Recovery (cEOR) processes comprise a number of techniques which modify the rock/fluid properties in order to mobilize the remaining oil. Among these, surfactant flooding is one of the most used and well-known processes; it is mainly used to decrease the interfacial energy between the phases and thus lowering the residual oil saturation. A novel two-dimensional flooding simulator is presented for a four-component (water, petroleum, surfactant, salt), two-phase (aqueous, oleous) model in porous media. The system is then solved using a second-order finite difference method with the IMPEC (IMplicit Pressure and Explicit Concentration) scheme. The oil recovery efficiency evidenced a strong dependency on the chemical component properties and its phase behaviour. In order to accurately model the latter, the simulator uses and improves a simplified ternary diagram, introducing the dependence of the partition coefficient on the salt concentration. Results showed that the surfactant partitioning between the phases is the most important parameter during the EOR process. Moreover, the presence of salt affects this partitioning coefficient, modifying considerably the sweeping efficiency. Therefore, the control of the salinity in the injection water is deemed fundamental for the success of EOR operations with surfactants.


Introduction
Oil is and has been for the last century the main energy source, and the economy on a global scale still depends largely on it [1][2][3][4]. After the traditional exploitation techniques, which include natural driven mechanisms (primary recovery) and waterflooding (secondary recovery), up to 50% of the Original Oil in Place (OOIP) might be extracted from the oilfield [3,[5][6][7]. Additional production techniques, or Enhanced Oil Recovery processes, target the remaining oil still trapped [8]. The most commonly used EOR methods are [6]: thermal (fire flooding-continuous/cyclic steam injection); chemical (polymers, surfactants and/or caustic); miscible (CO 2 or miscible solvent); and others (e.g., microbial EOR). The cEOR techniques are one of the most used for low and medium viscosity crude oils and they can be employed in a wide range of rock formations. Another important advantage of cEOR processes is that the chemicals to be used as sweeping agents can be synthesized or improved based on the desired characteristics for determined crude oils, rock formations and porous media conditions in general, optimizing the overall process and increasing the operational life of the reservoir. The objective in this paper is to present the numerical simulation of surfactant EOR flooding. Thus, a novel numerical simulator is presented to analyze tertiary recovery processes in a 2D oilfield [9][10][11], simulating a two-phase, four-component flow.
in an accurate, although relatively simple and robust way, solving the inconvenient described in previous simulators with respect to the phase behaviour. The former is assumed to be a function of the total dissolved solids (TDS) present in the water phase, using a piece-wise function based on previous literature. The phase behaviour modeling was done increasing the number of factors used to describe the chemical behaviour, but avoiding creating a system of non-linear partial differential equations which would complicate the simulation process. This leads to a new set of optimum design parameters which can be used during the synthesis of future surfactants.
The combination of these concepts has resulted in a novel simulator, which can be used for the design and screening of new surface active agents to be used in EOR. Nevertheless, the 2-dimensional modelling does not allow to consider the volumetric sweeping efficiency in EOR agents and how the difference in the densities between crude oil and the EOR agents could affect the recovery efficiency, specifically when the vertical permeability cannot be neglected with respect to horizontal ones. Furthermore, the influence of the temperature and energy transfer was not considered in the present model, which mainly affects the phases' rheological behaviour as well as the influence of the salinity in the surfactant flooding. Besides, the model proposed is limited to two phases, aqueous and oleous, neglecting the presence of the gas phase, which is usually present in many reservoir (e.g., gas cap). The mentioned phenomena should be the objective of further developments, tackling out one of its current disadvantages with respect with commercial reservoir simulators. Simulation of multiphase and multicomponent flow in porous media requires solving a number of coupled, highly non-linear equations dealing with temporal and spatial derivatives of pressure and mass concentrations. It is adopted in this model a second-order accuracy discretization scheme, along with a Total Variation Diminishing (TVD) flux limiter, which reduces the occurrence of numerical diffusion and dispersion mechanisms [11]. The compositional flow is a suitable approach to study the different cEOR techniques, which can be described as the mass transfer of a number of components (e.g., polymers, surfactants, salts, etc.) in a two-or three-phase systems. Generally speaking, commercial simulators present more advantages that academic ones. However, the objective of this paper is to present a novel cEOR simulator. The second-order accuracy in the discretization allows obtaining better results than standard academic simulators and it is comparable to commercial ones. Chemically speaking, the adaptation of the ternary diagram to model the phase behaviour (and the influence of the salinity in the calculation of the partition coefficient) and the phases' rheological behaviour (e.g., the calculation of aqueous viscosity in a step-wise approach considering the influence of all the components) were never considered in academic simulators and some of them are not considered by commercial ones. Thus, mathematically speaking it represents a novelty with respect to academic simulators and chemically speaking, it represents an advance when compared to academic simulators and presents some minor advantages when compared to commercial simulators. A major disadvantage with respect to these is the number of dimensions analyzed, although it is considered that a 2D approach is appropriate for a first version of the simulator in order to study the behaviour of chemical agents in EOR. Moreover, this approach allows to study the chemical performance in different oilfields in order to determine a set of desired features when synthesizing new products.

Physical Model
The physical model represents a reservoir (Ω) of known geometric dimensions, which has an absolute permeability (K) and porosity (φ). These can be constant or represented by a normally distributed fields (Figure 1). The flow is considered isothermal, Newtonian, incompressible and 2-dimensional (it is assumed as in previous models that the vertical permeability is negligible when compared to horizontal permeabilities); it is also considered that the system is in local thermodynamic (phase) equilibrium. The Darcy's law is valid and gravitational forces are neglected [38]. Surfactant EOR flooding, such as other chemical techniques, involves the flow of fluids in two phases (aqueous and oleous), and several components (water, salt, surfactant and petroleum). It is noteworthy that these components can be mixtures of a number of pure components, since for instance, petroleum is a mixture of many hydrocarbons, water may contain dissolved monovalent and divalent salts, and the surfactant is composed by a number of different molecules (e.g., cosurfactants) [3,8,16]. The recovery process consists in injecting during a first stage an aqueous solution with the surfactant, followed by a water bank (water or brine) in order to drive the chemical slug. This model is represented by a system of strongly non-linear differential equations, complemented by a set of algebraic relationships representing physical properties of the fluid and the rock, which are: component partitioning as a function of the salinity, interfacial tension, residual phase saturations, relative permeabilities, rock wettability, phase viscosities, capillary pressure, adsorption on the formation, and dispersion. It is clear that most of these properties are strongly dependent on the surfactant concentration, which alters the flowability of the aqueous phase. The interfacial tension depends on the former and it affects the residual saturations, which are dependent on the capillary number, a dimensionless parameter based on the IFT and the viscosity.

Mathematical Model
The compositional approach offers the versatility to model a number of components in porous media affecting the properties of the phases involved. In a previous paper three different components were considered, whilst in this new model the salt is included as a fourth component. This does not modify the equations, but only the number of auxiliary relationships to determine the system.
The mathematical relationship, defined by N comp · N phases − 1 , can be derived from the compositional approach in order to determine numerically the system.
The numerical simulator presented in this paper is based on previous ones developed also for EOR flooding [38,53]. The 2D model was validated against academic and commercial simulators in a previous paper, discussing also its order of convergence ( Figure 2) [9,11]. Since this model it is an extension of a previous one, it is considered that the validation process previously done is valid for this model [54][55][56]. Thus, the Darcy, mass conservation and aqueous pressure equations describing the system are presented in Equations (1)-(3), φc r ∂p a ∂t

Discretization of the Partial Differential Equations
The differential form of the Darcy's and mass transport renders a system of non-linear parabolic partial differential equations, discretized using a Finite Difference Method scheme. The aqueous phase pressure is implicitly solved using a centered scheme and a second-order Taylor approximation for the time derivatives. Moreover, the total and aqueous Darcy velocities are also discretized and explicitly solved using a centered difference stencil. Therefore, Equations (1) and (3) Finally, the discretization of the mass conservation is presented. Equation (2) is the advection-diffusion PDE used in hydraulic, mass transport and flow processes in porous media. The advective terms have an hyperbolic nature and the traditional first-order Upwind schemes cause numerical diffusion. In this model a complete second-order explicit discretization scheme is employed based on flux limiting techniques. This increases the numerical accuracy of the simulator and decreases the influence of numerical phenomena on the results [11,25,[57][58][59]. Thus, the discretized mass conservation equation is, x,m,n ∆t 2φ

Boundary Conditions
At the beginning of the EOR operations, the oil saturation in the oilfield is assumed to be equal to values after a primary recovery, or the saturation after a waterflooding which reached a threshold fractional flow at the producing well, rendering the exploitation no longer profitable. The pressure is assumed to be constant and moreover there is no chemical agents present.
The EOR process starts with the injection of the chemical agent (surfactant) with constant concentration. After this period, a water-bank follows in order to sweep the remaining oil. The boundary condition for the oilfield consists in a 'no flow' on the contour (Γ), since it is assumed that the porous formation is surrounded by an impermeable rock layer. Regarding the advection, this boundary condition is satisfied by imposing the mobilities to be zero.

Chemical Component Partition
The most important part of a numerical simulation using a compositional model is to understand how the components distribute into the phase, what it is called phase behaviour of the system. A surfactant flooding can be reasonably well modeled in a ternary phase diagram, wherein the chemical component is located in the apex, and water and petroleum occupy the lower vertices. The composition of a mixture is then determined by any point inside the triangle, which is known with two of the total concentrations of the components [33,37,60]. The numerical simulation of this model comprises a two-phase, four pseudo-components system, which can be modeled using this ternary phase diagram. If the concentration of the surfactant increases, oil and water become miscible, hence the triangle can be divided into two zones: a miscible in the upper part and the immiscible in the bottom. The curves delimiting these regions are determined by volumetric concentration ratios calculated as, Depending on the value of k c , two different behaviours can be observed: Type II(−) (for k c < 1), and Type II(+) systems (for k c > 1). These parameters affect greatly the efficiency of the EOR process. They should be correctly measured to predict recovery efficiencies with existing surfactants, or be the cardinal values when the objective is to design new chemical products. The partition coefficient depends on the composition of the surfactant as well as on the water characteristics, such as temperature and salinity. Sheng [51] described this two-phase approximation to surfactant flooding and the dependence on the salinity of the partition coefficient. He presented an empirical, piece-wise function (Equation (15) and Figure 3) to describe this, which also depends on the optimum volumetric salt concentration in the aqueous phase (V a s,opt ) [49], Figure 3. Partition coefficient dependence on the salinity. In this example, according to Equation (15), the optimum volumetric concentration of salt in water (V a s,opt ) was adopted from Sheng [51] and it is equal to 2%.
In this case is assumed that all the salinity present in the reservoir is composed by monovalent cations and that the salt component is only present in the aqueous phase (V o s = 0). With these extra relationships the physical system becomes, mathematically speaking, determined with a unique solution, and the parameters previously introduced can now be calculated for each representative elementary volume (REV). Figure 4 shows the original and simplified ternary diagrams used for this simulation.

Interfacial Tension
The IFT of the system depends on the presence and the concentration of the surfactant as well as the emulsion type present in the reservoir. In this simulator, a simplified IFT correlation is used [32,37,38,52,53]. For Type II(−) systems (oil/water emulsion): For Type II(+) systems (oil emulsion/water): where G 1 and G 2 are input parameters, while the term F is the correction factor introduced by Hirasaki [32], In cEOR processes, the surfactant decreases the interfacial energy, allowing the mobilization of the oil trapped in the reservoir. Thus, the residual saturations of the phases depend on the IFT ( Figure 5). The IFT of the water-oil system (with no surfactant present) is assumed to be constant throughout the simulation.

Residual Saturation
The residual saturations represent an important value in conventional and enhanced oil recovery processes. They establish a limit to how much oil can be mobilized by the sweeping agents. If these values can be reduced, this will increase the efficiency of the whole process. As explained before, they depend on the IFT in the water-oil system. The presence of the surfactant in an EOR flooding can modify the residual saturations in porous media. This relationship is ruled by a dimensionless group, the capillary number, The functionality between this group and the residual saturation is described by the following model [38]: The piecewise function is defined by a set of constant parameters which depend on the fluids and the rock formation. The relationship between the residual saturation after EOR and waterflooding processes is known as normalized residual saturation of phase j. The Equation (20) determines what is known as capillary desaturation curves ( Figure 6). At low capillary numbers, the behaviour is similar to a secondary recovery and the normalized residual saturation does not decrease. As the IFT decreases and/or the viscosity increases, the capillary number raises to values several orders of magnitude higher than those of secondary recovery. It is for this reason that in areas of high speeds (i.e., nearby the wells) oil saturations lower than those attained in a waterflooding can be achieved. As can be seen, the aqueous phase requires much higher values of N vc to achieve a full desaturation [8].

Relative Permeabilities
The relative permeabilities influence the phase velocities in the Darcy's equation, and thus the efficiency of oil recovery depends greatly on them. These are functions of the residual saturations calculated in the previous point. The model adopted to calculate these permeabilities is taken from Camilleri [41,42], which is generally used for most chemical flooding processes. Knowing beforehand the phase saturations, the relative permeabilities are calculated according to, where k j0 r and e j represent the end point and the curvature of the function k j r (S j ). These values are calculated by the following equations: where k j0H r and e jH are the endpoint values of curvature and relative permeability function system for water-oil without the presence of chemical agents, respectively.

Phase Viscosities
The viscosity of each phase depends on its composition as a function of the volumetric concentration of each component, according to the following function [38,41,42]: where α k are constants, and µ aH and µ oH are the viscosities in the water-oil system without surfactant, respectively. Nevertheless, in this work a modification is proposed to the formulation presented in simulators (e.g., UTCHEM) and the model developed by Porcelli [37], introducing the dependence on the salt in the pure water/brine viscosity and replacing the latter in Equation (24), according to the following expression [61], where A sal and B sal are constants based on rheology experiments. For the oil phase, as did before, a Newtonian behaviour is adopted. According to the literature [62], light and medium oil cuts exhibit Newtonian behaviour while heavy oil might present a slight shear-thinning region.

Introduction
The tertiary recovery process starts with an initial period in which an aqueous surfactant solution is injected, followed by a water bank used to sweep the surfactant slug and oil bank towards the producing wells. The chemical concentration profile has a shape similar to a wave propagation in a two-dimensional field, travelling at a speed equal to Darcy velocities in each block. For the first part of simulations, the influence of the salt will be neglected in order to study the phase behaviour (three-component system). Subsequently, a set of simulations is analyzed comprising the influence of the fourth component (salt) on the oil recovery mechanisms during surfactant flooding, establishing guidelines for the synthesis of new surfactants.
The initial conditions are set to consider that the surfactant EOR process starts after a waterflooding (or secondary recovery). This is aimed at studying how the surfactant can mobilize the remaining oil which is trapped and can no longer be swept by conventional recovery techniques. Thus, this initial oil saturation allows to study how the surfactant modifies the oil mobility.

Data
The simulation conditions and physical properties are established aimed at emulating an existing oilfield apt for a EOR recovery process with surfactants (Tables 1 and 2). These parameters were used during the simulator validation process as well as in previous papers in which different oil recovery processes were analyzed [9,11,37,38].

Three-Component System
This part of the study was done considering the salt is not present in the system, with no influence on the partitioning component and phase behaviour. The goal is to determine the influence of these parameters on the oil recovery process. Subsequently, this fourth component was added to the system in order to determine its influence on the abovementioned parameters. The latter will be used as reference case for the analysis when salt is considered in the mathematical model.

Influence of the Component Partitioning Parameters
The first part of the analysis, and the one deemed as the most important in surfactant flooding, is the study of the influence of the phase behaviour parameters on the three-component system. This importance has already been discussed previously [33,41,42,44,52,[63][64][65][66][67][68], and it plays an important role in the recovery process, because it determines the fractions of chemical, oil and water in each of the phases. It is noteworthy that, while these parameters are dependent on the salt content present in porous media, in this section this dependence will not be considered since this component is not taken into account. The objective in following this approach is to study firstly the influence of these parameters and determine, if possible, a trend and relationship between the values of these and the oil recovery. Subsequently, the salt is included as fourth component and also its functionality with the partition coefficient.
The partition coefficient was considered by several authors as the most important parameter to study the phase behaviour and hence it was chosen as the first to be evaluated. This parameter determines the type of emulsion present: lower values (k c < 1) render oil-in-water emulsions, known as Type II(−), whilst for higher values it renders water-in-oil emulsions, also known as Type II(+). Several recovery processes were simulated keeping constant the other two parameters (L a pc = L o wc = 1), obtaining the results presented in Table 3 and Figure 7.  Table 3 and Figure 7 show a clear relationship between oil recovered and the partition coefficient, a trend that agrees with the results published by Porcelli [52] and Bidner [38]. Numerical simulations were stopped after the chemical flowrate decrease to negligible values, which virtually stops any further oil recovery. For partition coefficients larger than the unity (k c > 1), complete desaturation profiles are almost achieved in the area around the injection well, which also coincides with the results of the 1D model published by Porcelli (Figures 8 and 9). Due to increased sweeping efficiency, the time employed for the chemical breakthrough to take place is also proportional to the partition coefficient. As more chemical component goes into the oleous phase, the chemical bank is less dispersed and it rises above its injecting concentration [38]. This was also observed during the simulations, with maximum concentration values of surfactant (obtained at the injection node) proportional to k c . The oleous phase swells up with water and chemical, and a lower amount of oil is then a left behind after the chemical slug, leading to higher oil recovery efficiencies is as k c increases. In turn, Figure 7c shows the pressure drop observed during the whole simulation. Firstly, it is important to note that the maximum values obtained with surfactant are significantly lower than those with normal polymer EOR techniques due to its rheological characteristics. In addition, at the end of flooding process, there is no difference between waterflooding and surfactant EOR flooding because in this the case the smaller size of its molecules do not cause the disproportionate permeability reduction (DPR) phenomenon, commonly found with polymers. Furthermore, there is no significant difference between the different partition coefficients with respect to the final pressure. This is due to the fact that, when no chemical is present, the oil relative permeability becomes negligible and all the pressure drop comes from the aqueous phase, which is nearly at the saturation when its relative permeability reaches its maximum. This renders that pressure drops are approximately the same at the end of the flooding. With respect to the time evolution, there is an initial decrease during the injection period related to the fact that surfactant in the injection well region starts mobilizing the trapped oil. Afterwards, as the oil bank is formed and displaced towards the production well, there is a gradual increase of the pressure drop. When the oil bank (along with the surfactant) finally reaches the production well and the total saturation in the reservoir decreases, there is a sudden decrease in the pressure drop. The minimum in the latter is also a function of the partition coefficient: as the k c is increased, the oil recovery is improved and therefore the oil saturation in the reservoir decreases, thus the pressure difference between injector and producer is also diminished. Finally, where the surfactant concentration is lower than the critical micelle concentration (CMC), the capillary number decreases due to an increase in the IFT, and thus the relative permeability values become similar to those at the beginning of the EOR process. Moreover, the pressure drop returns to the initial values. Figure 7b presents the flows for each component as a function of time. As mentioned above, the chemical bank suffers a delaying as k c increases due to increased sweep efficiency. It is also noteworthy in this graph that the total injected chemical differs in all the cases. When the simulations were performed, the injection time was selected so as to have approximately the same maximum concentration of the chemical slug in the injection well region in the reservoir. Nonetheless, this trend was verified by simulating these scenarios injecting the same total amount of chemical and the results, although differ to those presented in this paper, show the same behaviour: higher partition coefficients yield higher recovery efficiencies. Figures 8 and 9 depict the oil saturation profile, surfactant concentration, IFT and pressure during simulations for different partition coefficient cases.   Therefore, the viscosity does not vary significantly due to partition coefficient, nor it reached similar values to those obtained in polymer flooding processes. It is also due to this reason that surfactants must be used in petroleum reservoirs with low/medium viscosity, so as not to diminish the sweep efficiency because of increased mobility ratios. Finally, top right (b) plots depict the IFT and pressure profiles, respectively. In the absence of surfactant, the system yields the water-oil IFT, while ultra-low values are achieved in the chemical bank area, increasing the Capillary number and reducing the residual oil saturation. The sweeping efficiency in the vertexes with no wells is reduced due to the negligible pressure drop in that region, which renders the lowest velocities in porous media. Following the study of the influence of the partition coefficient, the attention is set on the other two parameters: the solubilization and swelling coefficients. The partial miscibility of both aqueous and oleous phases due to the presence of the surfactant is quantified by these parameters. The former is generated either by the solubilization of the petroleum component into the aqueous phase or by the swelling of the oleous phase with water and chemical components [38,52]. Taking into consideration this goal, the following simulations were performed for each of the previous partition coefficients (Table 4), varying the remaining two parameters in order to estimate their degree of influence on the oil recovery process (Figure 10).  The summary of influence of the phase behaviour parameters is presented in Table 4. As mentioned previously, the partition coefficient is the most important factor in terms of recovery efficiency. This determines the type of emulsion (water-in-oil or oil-in-water), affecting the phases' flow behavior and the amount of oil that is mobilized by the presence of the surfactant. The solubilization and swelling coefficients also affect the amount of oil recovered, though to a lesser extent. From these two, the results showed a bigger sensitivity to the solubilization coefficient, whilst the importance of the swelling coefficient was more appreciable in Type II(+) emulsions.
The results of the sensitivity to the swelling and solubilization parameters show similar trends to those reported by Bidner [38]. Both define the area and shape of the interior triangle (two-phase zone) in the ternary diagrams ( Figure 4). The position of the apex of the triangle, its area, and the distance between the injection point and the binodal interior triangle along the extended tie-line which passes through the injection determine the sweeping effectiveness. Furthermore, the values of these parameters affect the IFT: low values of L a pc and L o wc increase the IFT, which reduces the capillary number and therefore increases the relative desaturation of the system.
As in the study of the partition coefficient, the oil recovery resulted proportional to the value of these parameters (Table 4 and Figure 10). Nevertheless, in the case of the swelling coefficient, this influence was less significant. In the case of Type II(−) emulsions, the values of oil recovery in the case of L o wc were substantially similar, and showed greater sensitivity in Type II(+) emulsions due to the reasons mentioned above in this paragraph (Figures 11 and 12). In the case of the solubilization coefficient (L a pc ), the values obtained showed greater sensitivity. This is because the distance between the injection point and the binodal interior triangle along the extended tie-line, which passes through the injection, varies with the values of L a pc , which does not occur with L o wc [38]. This decreases the influence of the latter in the recovery process. As conclusion of this analysis, the trends and sensitivities shown were confirmed for the three parameters defining the phase behaviour. The values found in the 2D model are consistent with those reported [9,38,53,69]. As mentioned in the introduction, the phase behaviour plays a fundamental role in surfactant flooding and therefore the parameters must be carefully studied and chosen according to the chemicals being used. However, it is important to mention that these values depend on the salt content (mono-and divalent) present in the reservoir. This was not considered in previous simulations using this simplified model.
All in all, based on previous results [9,38,69], it is considered that the main difference between them and this simulator is the treatment of the last two parameters: the simplified model treats them as constants whilst it is evident from Figure 4 that their values change as function of the concentration. This is because in the model presented there is a direct proportionality between the recovery efficiency and these coefficients. In order to match the results of the simulations with the observed values, it is necessary to modify this and establish the functionality with the concentration and the partition coefficient. Then a set of functions would be obtained defining all the phase behaviour parameters based on the chemical species' properties, thus they can be adjusted in order to achieve a matching between the simulation results and the laboratory tests. This is deemed a critical point in the model to be addressed in future research, in order to achieve a full dependability of the parameters on the system conditions. The partition coefficient will take the form presented by Sheng [51], and the solubilization and swelling coefficients will be represented by functions of the concentrations and the partition coefficient. Such a system would allow, for example, obtain maximum oil recovery values for salt content near the critical value, as it was published by several authors [49].

Influence of the Partition Coefficient in Random Permeability Fields
The second part of the study verifies and extends the previous simulations to a case with a heterogeneous absolute permeability field, as found in real oil reservoirs. Based on a previous work [11], a Gaussian random permeability field was adopted to model the heterogeneities in the oilfield. However, it is well known that permeability field usually follow a log-normal distribution. The objective of this simulator was not focused on the geostatistical reservoir-modeling but on analyze the influence of heterogeneities on the EOR process. The upgrade and modification of the permeability field distribution to account for more real oilfields will be considered in future studies. The influence of this in the performance of the EOR agent can be studied using different distribution laws, as presented in this paper. For this study it was generated an artificial permeability field according to a Gaussian random distribution ( Figure 13).
The operating conditions and physical properties of the phases involved were kept constant from the previous simulations. An additional difference, besides the random permeability field, is the number of elements used to discretize the domain, since the simulation in this study was done in a field with refined grid (40 × 40 elements). Table 5 and Figures 14-16 present the results of this analysis. Absolute Permeability (mD) Figure 13. Permeability random fields [mD] with normal distribution in X (left) and Y (right) axes.
Oil Recovered (-) Figure 14. Oil recovered in a random permeability field as a function of time.
As expected, the heterogeneity in the absolute permeability affected severely both the recovery and the sweeping efficiency, rendering a major part of the oilfield to not be contacted by the EOR agent ( Figure 15). It is well-known the effect of such heterogeneity in secondary and tertiary recovery techniques, causing that large areas of the oilfield may remain unaffected by the recovery process. Figure 15a-d show the mentioned phenomenon in the results from the simulator. Nevertheless, the field permeabilities did not alter the proportionality between the partition coefficient and the oil recovery efficiency. The variations in permeability primarily affect the surfactant concentration reaching every cell of the reservoir, which can significantly alter the results with respect to a constant permeability field. The heterogeneity in the permeability field did not affect however the trends regarding pressure drop, chemical breakthrough or residual saturations, which are similar to those presented above, although it did affected the absolute values of the variables mentioned. Although the permeability does not depend on the characteristics of the surfactants used, it should be taken into account during previous operations to determine an optimum well-pattern distribution, because a non-optimal arrangement in the wells location can cause a major by-pass by the chemical slug of an important part of the OOIP, reducing the effectiveness of the EOR method employed.

Four-Component System
The second part of this study comprises the analysis and discussion of an EOR process considering a fourth component, namely the salt. On this new simulator no distinctions are made between monoor divalent salt ions. The salt changes the phase behaviour in systems with surfactants [49,[64][65][66]70], so this is crucial to be taken into account during the simulation and laboratory or field experiments, in order to understand the actual behaviour of the system. The partition coefficient is affected by the latter according to Figure 3 [51]. To model this in the simulator, the last curve was adapted to a function (Equation (15)), dependent on the salt content (and the critical salt content, when k c = 1). During this paper an initial constant salt content was assumed in porous media.
The scope of the study is then to analyze the influence of salt content on the two phenomena flow, specially in the phase behaviour and the adsorption of chemical species onto the rock formation. The simulations performed in this section considered porous media with several initial salt contents, both below and above the critical salt concentration. This describes the possible two-phase systems in porous media, Type II(+) and II(−), together with different contents of salt in the injection water (brine), in order to study the influence of salinity with the injection process.

Influence of the Salt Component on the Phase Behaviour
This analysis and its subsequent discussion is divided into two different streams: how the salt present in the reservoir affects the recovery (with no salt content in the water being injected), and how the salinity in the water injection influences the sweeping process. In the mathematical model the salt is only present in the aqueous phase, so the transport phenomenon of this component depends only on the properties of the abovementioned phase. A number of simulations were performed with different initial salt contents (constant) in the reservoir, which, according to Equation (15), generate different phase behaviours and types of emulsions (Table 6 and Figure 17). The injection water contains only surfactant and the water itself. Oil Recovered (-) Figure 17. Oil recovered as a function of time for different initial salt concentrations in the reservoir.
The results demonstrate that, although the salt content has influenced the recovery efficiency, the difference in oil recovered between the different initial salt contents is not as significant as in the previous cases, for different partition coefficients. This is mainly due to what was exposed above: the salt component is present in the aqueous phase and travels faster than the oleous one, so that the salt content in the surfactant bank decreases rapidly, and thus it does the partition coefficient. However, the salt content increases the sweeping efficiency, which agrees with previous results [64][65][66]. This is shown in Figures 18 and 19, where the partition coefficient and the oil saturation are plotted. These simulations show that not only the content of salt in the reservoir is important for the effectiveness of the EOR process, but also the salinity in the injection water influences the performance of the whole process. Therefore, the second part of this four-component study is devoted to present and discuss the results of simulations when brine is being injected and how that can significantly impact in the oil recovery efficiency.
Partition Coefficient (-) The influence of salinity in the water injection has proved to be more significant than the initial one in the reservoir. This is because a constant input of salt in the water phase, along with the salt already present in the formation, modifies the properties of phase behaviour throughout the reservoir. Different configurations and concentrations of brine being injected were simulated in order to study the influence on the process efficiency. The results are presented in Table 7 and Figure 20.  Regarding the influence of salt in the reservoir, for the same concentration in the injection water, the same trend is observed as at the beginning of this section: the salt content influences the recovery efficiency, but it is not as significantly as the other parameters. This is valid for both Type II(−) and II(+) emulsions. However, the most important difference occurs with the injection water. Small variations in the salinity affect the oil recovered, and this is accentuated for Type II(+) emulsions, where the partition coefficient is more sensitive to the salt content (Equation (15)). Slightly increasing the salt content causes a substantial increase in the partition coefficient and the oil recovered (Figures 21 and 22). This is clearly appreciated in Table 7: the case with a salt concentration in the injection of 1.5% incremented substantially the recovery efficiency due to the formation of a type II(+) emulsion as well as the time required to sweep the reservoir (compared to the other cases). The reason is that the injected salt keeps the partition coefficient in non-zero values. In the previous case the salt, being only present in the aqueous phase and moving faster than the oil, was transported faster than the surfactant bank and that generated low or negligible partition coefficient values. This was avoided with continuous injection of low salt concentration brine. The effect of the salt injected is clearly appreciated in cases with initial salt content equal to 1.5%. A 3-fold increase in the salt content caused an oil recovery growth of 2-fold, due to the significant increase in the partition coefficient (Table 7). Therefore, it is advisable in oil recovery processes with surfactants the injection of low salt concentration brine in order to improve sweeping efficiency, which agrees the results reported by several authors [47,[64][65][66]71].

Conclusions
This study aimed at designing a new two-dimensional simulator for a four pseudo-components, two phase flow, applied to a surfactant enhanced oil recovery process. The problems involved, both physical and economic, were formulated and discussed, concluding that new methods of exploitation are necessary to improve the current oil production, until new energy sources are developed. The physical model is described by a system of non-linear partial differential equations, solved by a Finite Difference Method, using an algorithm implemented in MathWorks MATLAB. The higher-order numerical scheme coupled with the TVD flux limiter functions allowed reducing the numerical errors. The compositional approach proved to be a versatile tool performing both secondary and tertiary recovery processes without complications. The phase behaviour is represented by three parameters: the partition, solubilization and swelling coefficients. Along with these three parameters, there is a set of ancillary values which must be carefully defined so as to match the experimental results with the simulation. These values depend on the surfactant chosen (e.g., CMC), so experiments must be carried out beforehand to determine these properties accurately. Besides, further research is deemed necessary to determine a set of mathematical functions describing the phase behaviour parameters in order to match the theoretical results with the laboratory or field tests.
The simplification of the components ternary diagram proved to be an useful and valid approach to model the surfactant properties and phase behaviour, which evidenced to play a significant role in the EOR method. The partition coefficient showed to be the most important parameter in terms of oil recovered, followed by swelling and solubilization ratios. When the surfactant formed Type II(+) (water-in-oil) emulsions, the highest recovery efficiencies were obtained. Furthermore, the influence of the salt only proved to be relevant when it was injected (as brine) in the oilfield, since only in these cases its value remained above the critical salinity in the area where the surfactant slug was present, increasing the sweeping efficiency. Surfactant EOR simulations presented the potential of cEOR methods to sweep the residual oil by means of reducing the IFT. The simplification of the ternary diagram may provide a simple way to design chemical surfactants with the desired properties before carrying out experimental or field tests. Nevertheless, even though the numerical validation has been done, future work is deemed necessary before a full field-scale application of this model. It is considered that lab-scale experiments should be done in 2D or 3D models with surfactants in order to demonstrate the applicability of this numerical simulator. Core flooding is also recommended, though it is advisable the usage of models in which the area and volumetric sweeping efficiencies can be also studied. This will allow determining the parameters in an accurate way before a field-scale application is pursued in order to design novel EOR agents.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: