Influence of physical and rheological properties of sweeping fluids on the residual oil saturation at the micro- and macroscale

Oil recovery processes depend on many factors that can be altered in order to maximize the sweeping efficiency in porous media, and one of these is the rheology behavior of the displacing agent. Furthermore, scales in the recovery process should also be considered: from the macroto microscale systems, in which capillary forces become predominant. It is also well-known the non-Newtonian behavior of polymer solutions used in Enhanced Oil Recovery (EOR) processes. This has been considered before, explaining how the polymer’s viscosifying properties enhance the displacing process. Recently, another property exhibit by polymer solutions started being considered: the viscoelasticity. The interaction between the (macro)molecules in the displacing phase generates a complex stress field which cannot be simply addressed by an increment in the shear viscosity. We present a 2D, multiphase simulation at macroand microscale of a recovery process with different fluid models, showing that viscoelastic fluids increase the recovery performance due to the extra stresses generated by the polymer molecules, up to a 15.4% when compared to traditional waterflooding techniques. The viscosity of the displacing phase affects indeed the recovery efficiency, and moreover, the results also evidenced that not only the bulk viscoelasticity, but also the interfacial forces play a vital role in the microscopic sweeping efficiency in polymer EOR flooding processes. This can be used when determining the properties of future EOR agents to be synthesized.


Introduction
The study of Newtonian and non-Newtonian fluids covers a broad spectrum of engineering fields including: pipeline transmission systems [1,2], the food industry [3,4], flow in porous and/or underground media [5,6], and the study of foams and emulsions [7]. Within the group of non-Newtonian fluids, it is of particularly interesting the research in the so-called viscoelastic fluids, which also comprises several areas of engineering, such as the flow of slurries [8,9], polymer melts [10], viscoelastic jets [11,12], the study of the circulatory system and other body fluids [13][14][15], and the employment of polymer solutions in Enhanced Oil Recovery (EOR) [16]. Polymer flooding is one of the chemical EOR methods, which aims at improving the rheological and physical properties from the displacing fluid so as to increase the efficiency of the recovery process.

Polymer flooding
Water soluble polymers have been used in chemical EOR processes during the last 50 years. The first work involving polymers as a possible which modifies the fluid flow and therefore, the oil recovery. Moreover, viscoelastic fluids exhibit properties which are not observed in Newtonian ones, such as elastic instabilities described by Shaqfeh [39], Pakdel [40], Khomami [41], and lately by Clarke [42], when analyzing the anomalous extra oil recovery by means of polymer solutions. Thus, the numerical modeling of high viscoelastic fluids presents serious challenges, as shown by De [43], Alves [44], and Poole [45]. The development of both physical models describing the relationship between the stress and strain-rate tensors, and numerical techniques to solve these systems is a topic of continuous research. Nevertheless, it is not only important to consider the bulk viscoelastic properties but also the interfacial viscoelasticity exhibit by polymer solutions and polymeric surfactants, which present a great potential as EOR agents [16].
In order to understand how to increase the microscopic efficiency, an analysis of the acting forces and their importance at this scale of magnitude should be performed. Capillary forces impede the mobilization of the remaining oil (a phenomenon known as capillary trapping), while the forces generated by the pressure gradient are responsible for the displacement. A way to establish a mathematical relationship between these two is by means of a dimensionless group known as the Capillary number [46,47]. If the pressure gradient is high enough, the oil droplet can be displaced. This pressure gradient depends mainly on the viscosity of the circulating fluids. It is because of this phenomenon that, during waterflooding processes, in order to mobilize the remaining oil the pressure gradient to be applied should be 4 to 5 orders of magnitude higher that the currently employed, which is mainly limited due to economic reasons. Thus, capillary forces predominate and oil remains trapped [7,19,48,49]. Possible strategies used in chemical EOR aims at lowering the interfacial tension using surfactants, or increasing the viscosity of the displacing fluid and, in doing so, the pressure gradient. However, according to some authors [20,23,46,[49][50][51], the latter, as applied in polymer flooding, is still not enough to mobilize the oil trapped. This statement contradicts the experimental results obtained micro-and macroscales [25,26,31,35,36,[52][53][54].

Aim of this work
The purpose of this paper is to understand the mechanisms affecting the oil recovery both at the macroscopic and microscopic level. Therefore, a new numerical model was developed aimed at investigating the mechanisms and effectiveness of polymer flooding on the residual oil after waterflooding. This was considered necessary since most of the literature found up to now is based on a steady state analysis of monophasic systems, or on a study of the viscoelasticity effects on structures or shapes that simulate or represent the residual oil in its different forms usually found in porous media. In order to tackle the problems described in previous models, we developed and studied a new CFD model, considering a transient recovery process of a multiphase system, and to study the importance of the scales involved, we included both micro-and macroscopic systems, representing different oil recovery scenarios. In doing so, it was possible to simulate more accurately the sweeping process of residual oil in macro-and microscopic models, and to theorize over the possible mechanisms present during the sweeping, using fluids with different rheological and interfacial properties. These go from Newtonian systems to several non-viscoelastic, non-Newtonian fluids, concluding with the simulation of a recovery with a viscoelastic EOR agent. The hypothesis of this study is based on previous literature stating that the extra stress field caused by the interaction among polymer molecules is also relevant for the recovery process. The goal of this paper is to demonstrate that also the interfacial characteristics between the viscoelastic solution and oil play a major role in the sweeping process of the remaining oil trapped in "dead-ends".

Mathematical model
The numerical model is described by the incompressible Navier-Stokes equations, including the surface tension component, and the continuity equation, which together describe the momentum and mass conservation equations. For this study both oil and displacing fluid are considered incompressible. The system of equations yields (neglecting gravity forces) [1,[55][56][57], Furthermore, the interfacial force is introduced, which was not present in previous models since they did not consider multiphase systems. This is calculated according to Eq. (3), where is the surface tension coefficient, is the interface unit normal versor, and is the Dirac Delta function, non-zero only at the fluid interface, and is the viscous/elastic stress tensor. The Level-Set parameter is also used to approximate the Dirac Delta by a smooth function defined by Eq. (5). In order to take into account for the effects of viscoelastic polymer solution, it is introduced in Eq. (1) the extra stress tensor component, Several numerical models allow calculating the extra stress tensor. The simplest ones used for EOR processes are: Upper convected Maxwell (UCM) and Oldroyd-B (Old-B). UCM is a simplification of the Old-B model which does not take into account the contribution of the solvent in the solution viscosity. This model was used by several authors in order to study the effects of viscoelasticity [55,[57][58][59][60][61][62][63][64][65]. In this study the Old-B model was chosen because, even though it does not include the effects of shear-thinning described in the literature [52,63], it allows achieving a good approximation to the flow of polymers in porous media. Moreover, it can also be quantified the solvent contribution to the final viscosity of the solution. The extra stress symmetric tensor is then defined by the following equation, wherěis the upper convected derivative of the stress tensor, is the deformation rate tensor, anďis the upper convected derivative of the deformation tensor. The variable is the relaxation time of the polymer and 2 is the retardation time. These two values define the viscoelastic properties of the polymer solution. In order to solve this system of equations different approaches can be used [59,66]: the most commonly used are elastic-viscous stress splitting (EVSS), explicit diffusion (EDIF), and solvent-polymer stress splitting (SPSS), used in this study, based on evaluate separately the contributions of the Newtonian fluid (solvent) and the non-Newtonian viscoelastic (polymer), and measure their contribution to the overall momentum equation (Eq. (6)).
= + (10) +̌= 2 (11) where is the Weissenberg number. For pure viscous fluids Weissenberg number is zero (no elasticity), whilst an infinite Weissenberg number would correspond to purely elastic response [10]. In practice, when the Weissenberg number is greater than one it is considered P. Druetta and F. Picchioni  as a high Weissenberg number for Oldroyd-B fluids. In addition, the dimensionless system of equations is presented, which it is commonly used in the literature [3]. The viscosity relationship is expressed by Eq. (16), thanks to which can be evaluated the contribution to the viscosity of the solution by the solvent and the polymer. Thus, the final system of equations is,

Geometric model
In order to evaluate the properties of the polymer solutions two models were evaluated in COMSOL Multiphysics ® to simulate the different scales in oil recovery processes. The numerical model and the adaptation to the software was based on a previous paper from Alves [44], who studied numerically (and validated) the monophasic flow of an Oldroyd-B fluid through a cylinder. The macroscopic model ( Fig. 1) was adopted because of the availability of experimental results in order to validate the model. It consists of a 2D cell representing different configurations of "dead-end" pores, a pattern commonly found in porous media. Even though the model was chosen for the study, it is important to note that these settings are not realistic when the width/depth ratio is analyzed with respect to the diameter of the entrance channel of the polymer solution (especially in the last three chambers). Due to the relationship between the thickness and the other two principal dimensions, we considered that a 2D model would represent a good starting point for the research. The extension to a 3D system is recommended in order to achieve a full understanding of the viscoelastic and interfacial mechanisms working in the flow cell.
Before proceeding with the numerical simulations, a detailed analysis of the meshing process was performed in order to maximize the performance of the computational system as well as to improve the accuracy of the model, especially in the so considered critical areas of the system. The refinement of this model is shown in Fig. 2: the P. Druetta and F. Picchioni  original mesh was firstly refined in the whole domain, followed by a specific refinement in the walls of the small top channel as well as in the chambers' walls, in order to capture in a more precise way the movement of the interface.
The microscopic model was imported from the COMSOL library. This model (Fig. 3) is based on the work developed by Sirivithayapakorn [67], who designed their experiments on the basis of scanning electron microscope (SEM) images of thinly sliced rock sections. The geometric pattern from the images was exported to a DXF file, which was ultimately imported into COMSOL.
During the simulations a sensitivity analysis was performed in order to determine is the mesh refinement no longer affect the numerical solution. Thus, the number of elements, nodes and mesh configurations are different for the macro-and micro models presented in this paper.
As it was done with the macroscopic flow-cell, we analyzed the mesh of the microscopic system before the sweeping processes. In this case, the original mesh ( Fig. 4 -left) was refined in specific walls, namely those in the dead-ends and in small throats, in order to keep an accurate track of the interface and a correct velocity and stress profiles at the throats ( Fig. 4 -right).
In both cases, it is deemed that future studies extending the domain to 3D models are necessary to have a full analysis of all the viscous, viscoelastic and interfacial phenomena acting during the recovery process. This is of the utmost importance in the macroscopic model in order to consider, for instance, the plastic cover of the system (see Appendix), and in the microscopic model to consider the different throat diameters and chamber thicknesses which will affect the velocity and stress profiles.

Boundary conditions
The boundary conditions were set also according to literature [59,66] and they are presented in Eqs. (17) to (19). A fixed velocity field (flowrate) and pressure are established for inlet and outlet, respectively.
The assumption of a fully developed flow can be guaranteed in the macroscopic model, although this is quite different in the microscopic system. Indeed, the inlet to this model was assumed to be fully developed flow but in reality this may be hard to achieve, since it is not completely known the entire pore geometry of the rock. Finally, a no-slip condition is applied for the internal walls, along with the orthogonality of the components of the extra stress (applicable only for Oldroyd-B fluids). In the microscopic model, a symmetry condition (i.e., no stress and no normal component of the velocity field) is set for top and bottom boundaries.

Macroscopic model
Numerical simulations were performed and thereafter compared with experimental results obtained with the polymers synthesized for EOR purposes [68]. The base case for comparison and further analysis was made using water, representing a typical secondary recovery process. The rheological and interfacial properties of the two phases affect the dimensionless groups in the Navier-Stokes equations (e.g., Reynolds and Capillary numbers). The viscoelastic properties in the model are controlled mainly by adjusting the Weissenberg number for the displacing phase. The selected strategy on how modeling the viscoelasticity of the polymer solution in this work is based primarily on modifying the relaxation time of the polymer which can be obtained experimentally. For HPAM-based branched polymers, Wever reported the relaxation time between 0.47-10.89 s [68].
The base case (waterflooding) was performed with standard rheological properties. As explained in the previous section, the boundary conditions were established as an inlet velocity field at steady state conditions = × (1 − 2 ) and a fixed pressure ( = 0) at the outlet. The final result of displacement is shown in Fig. 5 (top) in the complete 2D flooding cell (meshed with 27,848 elements and a meshing quality of 0.8741). Moreover, Fig. 6 shows the detailed simulation of the chambers 2 and 3 (meshed with 7144 elements and a meshing quality of 0.8936) compared to the results obtained previously [68,69]. The total recovery factor in the whole model was 61.51% and in the simplified model 31.3% (See Appendix). Due to what was previously explained, the displacement of oil in the last 3 chambers is negligible due to the width/depth ratio of "dead-end" compared to the diameter of the entrance fluid channel. The main inlet and outlet cavities (chambers 1 and 7) do not allow us to make a proper assessment of the effectiveness of the displacing fluid. The agreement between the experimental results obtained previously and the numerical model is satisfactory, especially in terms of oil recovered from chambers 2 and 3. However, the numerical model is based on a 2D field and therefore has its limitations. This is due to the lack of the plastic lid from the experimental model, which adds another no-slip boundary condition in the entire fluid field. Even though this itself cannot alter substantially the results, there is a related factor in the experimental set-up that might affect the process, namely the gap between the metal and the plastic cover. This gap (in theory 0.5 mm) is susceptible to the closure torque of the metal screws. An uneven torquing procedure could originate areas with different gap and this may affect considerably the results. This is something to be checked in future experiments with similar devices. Figs. 5 and 6 show the results of volume fraction of the displacing phase and pressure, respectively, for the three cases studied in the macroscopic model: water, shearthinning non-viscoelastic polymer solution and viscoelastic polymer solution.    7 shows the velocity field for waterflooding case basis. As it was expected, the highest rates are found in the upper channel crossing throughout the entire cell, while the penetration depth in the chambers 2 and 6 corresponds to that found in laboratory experiments (See Appendix). The penetration in the last 3 chambers is negligible when compared to the penetration in the first two "dead-end" channels. The pressure field (Fig. 6 -top) does not provide additional information in the case of waterflooding, with the higher pressure drops occurring in the upper channel, due to the higher velocities. Moreover, due to the difference between water and oil, there is a ten-fold decrease in P. Druetta and F. Picchioni the pressure drop between the initial condition and the end of the waterflooding process.
The first series of simulations allowed us establishing a validation between the numerical model and the laboratory results. A Newtonian fluid was used to simulate a sweeping process in a macroscopic 2D field with different "dead-end" chambers. However, it is evident that the shear-rate is not the same in the whole device, such as it happens in porous media. The next step consists then on investigating what happens when a Non-Newtonian viscous fluid is employed to perform the same sweeping (e.g., a non-viscoelastic polymer solution). A similar process to the previous one is carried out, but using a non-elastic viscous fluid with similar rheological properties to those found in the literature for polymer solutions [23,68]. The model chosen to simulate the rheology is the Carreau fluid, which comprises a part at low shear rates at which the polymer solution behaves like a Newtonian fluid, with the viscosity depending, among others, in the concentration and the polymer molecular weight. Then, at moderate shear-rates, a shearthinning region, similar to the one described by the Power-Law model. Finally, at higher shear rates, the viscosity has the asymptotic limit in the viscosity of the solvent, which in this case is water (Eq. (20)).
where represents the limit viscosity of the solution at a infinite shear-rate; 0 represents the viscosity of the solution at zero shear-rate, is the power law factor, and is the shear rate. The rheological behavior of the solution chosen for the simulation is shown in Fig. 8 [68]. Fig. 5 (middle) shows the final result of the oil recovery process, in terms of the fluid volume fractions. The sweeping efficiency in the whole model was 66.73%, which represents an increase of 5.22% with respect to the waterflooding. As in the previous case, the result obtained in the last three chambers is negligible and cannot be used to judge the performance of the viscous polymer solution. The chambers 2 and 3 show an increased recovery, which results in a greater depth, although not significant.
Figs. 6 and 9 (middle) depict the velocity and pressure fields, respectively. As expected, the pressure drop increased approximately 5-fold due to increment in the viscosity of the displacing fluid. The velocity field, set as a boundary condition at the entrance, has not changed significantly. The velocity profile in the upper channel was slightly modified due to the non-Newtonian behavior of the fluid. The only change worth of mentioning is in the chambers simulating the "dead-ends", where the flow lines achieved slightly deeper penetrations.
Finally, the recovery process with a viscoelastic fluid was simulated in the macroscopic model, according to the Old-B model. The simulation show that the presence of normal stresses generated by the elastic effects improve the sweep fluid oil, and these increase with the elasticity of the fluid, expressed in terms of the Weissenberg number ( ). It is important to note that all the physical properties of the polymer solution depend on the molecular architecture and weight, namely: the zero-shear viscosity, the relaxation time, radius of gyration, and degradation rate [70]. With respect to these, the molecular architecture has become increasingly studied in the development of new polymers for EOR, which exhibit increased viscoelastic and rheological properties when compared to traditional linear ones [16]. In the presented model the process was simulated with = 0.4 for the displacing phase and a dimensionless viscosity ratio * ∕ * = 0.1. As it was performed in the waterflooding analysis, the simulation was carried out with the complete model (Fig. 5 -bottom), followed by the refined model of chambers 2 and 3 (See Appendix). The results were compared with those obtained in laboratory experiments [68], finding a good correlation between simulations and laboratory tests. The total recovery factor in the flow cell was 69.92%, which represents an increment of 8.41% with respect to the base case. Moreover, in the simplified model the recovery factor was 38.90%, which represents a final increment of 7.60% with respect to the waterflooding. This extra recovery obtained by the viscoelastic fluid also shows good correspondence with the average results obtained in coreflood experiments. These showed a great variability and dependence on the polymer's architecture, average molecular weight and concentration, with recovery values ranging from 5 to 25% when compared to waterflooding processes [68]. It is also important to point out that all these polymer solutions exhibited the same rheological properties, in terms of viscosity, and storage and loss moduli, implying that the recovery efficiency is also dependent on other physical/interfacial parameters, as the results of these simulations also suggest.
Such as in the waterflooding, there is an agreement between the experimental results and the viscoelastic model in the flow cell. The penetration depth in chambers 2 and 3, as well as the interface shape show a resemblance with the pictures from the experiments (See Appendix). This also works as a validation for the mathematical description of the viscoelastic problem, complementing the previous validation done in the waterflooding. The velocity field is shown in Fig. 10. It can be seen that the viscoelastic fluid has entered the "dead-ends" at a greater depth than in previous cases, which is related to the viscoelastic properties of the fluid and is directly proportional to the elastic properties, among others. Figs. 6 (bottom) and 11 show the pressure field and the first regular stress difference, respectively. A result to be considered in this model is that the highest elastic stresses are in the upper channel connecting all the chambers, while in the "dead-ends" as well as in the main chambers the values are much lower. This is mainly due to the geometric scale relationship between the mentioned chambers and the circulation ducts and its influence on the velocity field. Therefore, it can be seen in Fig. 11 that the normal stress difference is more than one order of magnitude lower than in the upper channel, which has been also reported by others authors [48]. This occurs at low Weissenberg numbers, and during oil recovery processes in porous media this number should be significantly higher than pressure forces. However, for higher numbers Weissenberg, the effects of non-linearities begin to affect the numerical performance of the model, which has already been previously reported. There currently are lines of research dedicated to develop numerical methods so as to model fluid with high numbers. With respect to the values of the normal stresses 11 and 22 , the results show that the second normal stress is at least one order of magnitude smaller than the first normal stress component, which is in accordance with the results reported previously [61].
One of the conclusions of these simulations was the negligible influence on the oil recovery factor of the Reynolds number, which has P. Druetta and F. Picchioni  already been reported by several authors [48,52,56,[71][72][73][74]. Another factor that has not been considered in previous studies due to the type of simulation performed is the effect of the force originated due to interfacial effects. During the simulation it was observed that this term plays an important role in oil recovery along with the viscoelastic properties, so further research should be undertaken to capitalize both effects, the viscoelastic and interfacial, something that could be achieved by using polymeric surfactants in lieu of other chemical products [75,76].
All in all, three different fluids were tested in an macroscopic 2D flow cell, comparing the results with those obtained in laboratory simulations, showing a good correspondence between physical and numerical model. The waterflooding used as benchmark yielded the P. Druetta and F. Picchioni   lowest results in terms of oil recovered, followed by the non-Newtonian, non-viscoelastic fluid and finally by a viscoelastic Oldroyd-B fluid. Not only the viscoelasticity was important in these results, but also the relationship between interfacial and viscous forces, expressed by the Capillary number (Table 1). Regarding this dimensionless group, it is worthy of analysis what happens with non-Newtonian fluids, in which the viscosity, and thus the , changes throughout the pore geometry as a function of the shear-rate. This number was considered constant during the simulation for a shear thinning fluid, although in real-world applications it changes both as a function of time and the geometry. De [77] studied this phenomenon and the relationship during a sweeping process with the elastic instabilities that appear in viscoelastic fluids, which renders a critical above which the size of residual oil blobs decreases when viscoelastic fluids are employed.
An important point to consider in the simulation of viscoelastic fluids is the influence and presence of the so-called 'elastic turbulences'. In this case, the presence of the polymer molecules, with different architectures and molecular masses, renders the EOR agent a viscoelastic fluid, and part of the energy in the flow can be stored, depending on the amount of the elastic or storage component. The flow of viscoelastic fluids may generate purely elastic instabilities that have a different origin than the inertial ones, created at high Reynolds numbers. These instabilities are generated when the previously mentioned elastic energy stored is higher than the dissipation in the system caused by the polymer relaxation. All the parameters previously introduced are considered in the Weissenberg number; it has been reported that these elastic instabilities appear when > , with this critical a function of the polymer molecules and the process conditions. When this happens, these instabilities lead to what is known as 'elastic turbulence', which is defined as the random flow caused only by nonlinear elastic stresses, and therefore, it can be observed also at low numbers [78,79]. It was reported by Soulies [78] that this phenomenon usually takes place at > 1, and therefore it is considered negligible P. Druetta and F. Picchioni  for the simulations presented here. It is recommended that in future studies, at higher , the elastic turbulence should be considered and its influence, both at macro-and microscopic scales, evaluated. This will ultimately help in optimizing the design of polymer or other viscoelastic agents for EOR.

Microscopic model
The second part of this paper involves the microscopic model (Fig. 3) described previously, performing similar studies as those carried out with the macroscopic one. For this model the meshing was studied carefully because of the grooves and sharp edges present, which required a localized mesh refinement to achieve results independent of the discretization. Thus, the discretized model has 56,602 elements, with a mesh quality of 0.8684. The boundary conditions chosen were similar to the macroscopic model: no-slip condition on the inner walls; a constant outlet pressure and inlet velocity fields. The upper and lower boundaries were modeled as symmetry axes [58][59][60]66].
The first simulation (base case) performed was a waterflooding, using as fluids the same used in the previous model, in terms of physical and rheological properties. The final results of the flooding process are shown in Fig. 12. Water has displaced part of oil trapped in the porous medium but a water circulation channel was created, due to the difference in the viscosity between the two phases, reducing the efficiency of the recovery process. The interfacial properties of the system, as in the previous case, played an important role. Nevertheless, at these microscopic scales the role they play in the capillary forces is even more relevant, especially with ultra-low Reynolds numbers such as those encountered in reservoir simulation. Then, the Capillary number in the Navier-Stokes equation becomes now a preponderant term in the flow equations. The results in "dead-ends" in the model show a similar profile to the macroscopic model, with water unable to penetrate or displace the oil trapped therein. The recovery efficiency of the waterflooding process considered as base case was 67% of the original oil trapped in place.
The velocity field in the final condition of the simulation is shown in Fig. 13. The velocity of the displacing phase depends on the pore size and on the water-oil system interfacial tension, tending to form a parabola-shaped velocity profile that resembles the laminar flow in pipes. The pressure, as in previous cases, depends on the viscosity and velocity of the displacing phase and as it might be expected, higher pressure gradients are found in the pore throats. Due to the small scales and velocities involved, it is deemed not necessary to include in the report the graph showing the pressure gradient throughout the model.
Then, as it was done in the macroscopic model, a flooding using a non-elastic viscous fluid was performed, but with a rheology behavior similar to that found in the literature for polymer solutions at moderate and high shear rates, what is known as the shear-thickening phenomenon [23,68,80,81]. In this case the rheological model chosen, similar to Carreau at low and moderate shear rates, was developed by Delshad [82]. The causes of why the phenomenon of shear thickening occurs at high shear rates in polymer solutions is known and already P. Druetta and F. Picchioni  reported in the literature [23,51]. The polymer molecules flow through a porous medium composed by a series of small pore bodies and throats, and an elongational flow field occurs. When the shear rate is high enough, the polymer molecules do not have sufficient time in comparison to their relaxation time to stretch and re-coil. The result of these elastic strains causes an apparent high viscosity [34,51]. The unified viscosity model (UVM) is composed of two terms: one related to the phenomenon of shear thinning, which occurs at low and moderate shear rates, and an "elastic" term associated with the phenomenon of shear thickening (Eqs. (21) and (22)).
These parameters depend on, among others, the solvent viscosity and the polymer concentration. The rheological behavior is shown in Fig. 14. The parameters that define this behavior are as follows and were based on the literature [51]. The reason why this model was P. Druetta and F. Picchioni  not used in the macroscopic model is that the rates shear encountered during the simulation with water did not justify the use of the UVM model, since shear rates which may lead to phenomena of shear thickening were not registered. The final results of the flooding process are shown in Fig. 15. The recovered oil was increased, with a marked improvement sweep caused by the fluid in the area near the outlet, because, as in the previous case, the rheological properties and also interfacial phenomena between the oil and the polymer solution had a vital influence on the recovery efficiency. The formation of channels as in the case of waterflooding has disappeared and a more uniform sweeping of oil-saturated areas is registered. In this case a slight increase in the recovery efficiency was obtained after the flooding process, with a 71.1% of the original oil in place displaced, a 4.1% increase with respect to the waterflooding. Fig. 16 shows, as in the previous case, the velocity field. Due to the boundary conditions imposed on the system, this value has not changed significantly. The only change worth of mentioning is in the chambers resembling the "dead-ends", where flow lines achieved a deeper penetration, and also a decreasing in the number of channels caused a more uniform velocity field distribution.
The last simulation performed in this paper involves the analysis of the performance of viscoelastic fluids in porous media, considering at a microscopic scale the viscoelasticity exhibited by polymer solutions, as well as interfacial phenomena, which were not considered in previous works. Due to the factors already mentioned above, the fluid simulation with high Weissenberg numbers is still a challenge for existing numerical methods. Furthermore, if it is also considered an intricate geometric (as it is the network of channels in a porous rock formation), and the presence of more than one phase, the problem becomes a strong nonlinear system, with difficulties to find a numerical solution. The Old-B model was used as in the previous case to simulate the effects of elastic polymer molecules. In the presented model the process was simulated also with = 0.4 and a dimensionless viscosity ratio * ∕ * = 0.1. It is worth noting that, although the Old-B model allows us to consider the effects of elasticity and the contributions to the viscosity of the polymer and solvent (so UCM is a particular case of Old-B), it presents no shear  thinning behavior. Thus, the viscosity for any shear rate is constant, contrary to was what simulated with Carreau and UVM in microscopic and macroscopic models, respectively. The oil recovery was increased partly due to the presence of elastic stresses in the "dead-ends" and similar geometries, which increased the sweeping efficiency (Fig. 17). Finally, for the viscoelastic product the highest recovery efficiency was obtained, reaching a 82.4% of the original oil in place, which means a 15.4% increase with respect to waterflooding and an 11.3% with respect to the viscous non-elastic polymer model. The velocity field is shown in Fig. 18. It can be seen that it has entered the "dead-ends" at a greater depth than in previous cases. This is related to the viscoelastic properties of the fluid and it is directly proportional to the elastic properties of the same, as is the case in the macroscopic model, although it was expected that due to the length scale these phenomena to become more evident. Fig. 19 shows the first regular stress difference. As in the case of macroscopic model, the pressure forces are still more than one order of magnitude greater than the forces due to the elasticity of the polymer. This occurs at low Weissenberg numbers, and processes simulating oil recovery this number should be significantly higher. In porous media, the elastic forces can grow faster than the pressure ones and become of the same order of magnitude [29,82,83], and even become the dominant mechanism to recover the residual oil. According to Afsharpoor [48], flows with ∼ (100) may take place in porous media, especially at high speeds found near the wells.
In order to summarize the result of this section, Table 2 presents the oil recovered by different fluid models along with the dimensionless groups which characterize the recovery process. As shown in the macroscopic model, the waterflooding was used as reference case, and rendered the lowest results in terms of oil recovered. Then, the shear-thickening model (UVM) presented an increase in the recovery and finally the viscoelastic Oldroyd-B model yielded the best results. There is a clear relationship between the rheological, viscoelastic and interfacial properties and the oil recovered, as shown in the values of the dimensionless groups in Table 2.

Conclusions
The original theory that polymer flooding only increases macroscopic efficiency without altering the process at a microscopic level could not explain the results, both in field and laboratory tests, showing a decrease of the residual oil. Therefore, other phenomena also play a role in the recovery process, which were not taken into account in the former theories. Viscoelasticity proved to be one of them. The simulations carried out in laboratory experiments in a macroscopic 2D flow cell as well as the numerical simulations performed in the present study demonstrated the role of elastic phenomena, mainly at a microscopic scale. The viscoelastic fluids were able to recovery from 7.6% to 15.4% more oil than traditional secondary processes. This came a cost of an increase in the pressure drop, which is a factor that must be carefully analyzed if these processes are applied on a reservoir scale tests, due to the design of the injection facilities. While the effects of elastic forces are not dominant due to low Weissenberg numbers employed during simulations, they grow at a faster pace than they the pressure forces do. In a porous medium it is possible then that these elastic forces reach values greater than the pressure ones, and thus become the predominant driving force. It also resulted evident from the simulations that the interfacial term in the Navier-Stokes equation played an active and important role during the process of oil recovery. This phenomenon, quantified by the Capillary number in the interfacial force term, is influenced by the IFT and the displacing fluid's viscosity. The importance of this factor, which could not be evaluated in single-phase simulations, shows that the combined action of viscosity, elasticity and low interfacial tensions may be useful for the recovery of low and medium viscosity oils. Regarding the viscoelasticity, due to  the results obtained in this study, it is considered that not only the bulk viscoelasticity, but also the interfacial one should be study in order to assess its importance in the sweeping process with polymer solutions. Based on the aforementioned it is deemed necessary to develop further research in the field of polymeric surfactants, which it is considered to be very promising for EOR applications. With regard to the numerical model and its limitations, it is believed that, as mentioned by other authors, the development of numerical models that can handle instabilities and non-linearities at high Weissenberg numbers are needed as well as the extension to 3D domains, in which all the viscous and viscoelastic phenomena can be captured. Future work should be aimed at using of more complex viscoelastic models which take into account phenomena that the Old-B model does not consider (e.g., shear thinning, finite dumbbell extension) and to study the instabilities that take place at high Wi numbers. different dimensions. The bottom part of this device is made out of Aluminum, and it has a plastic, transparent top cover; the flow-cell has a depth of 0.5 mm. The procedure of these experiments consisted in filling the cell with the oil and then use different agents in order to test their recovery efficiency, especially from the dead-end chambers. Wever [68] reported the results using brine and a series of different linear or hyperbranched polymers, showing the improved efficiency of the last ones in the penetration and sweeping of oil from the dead-end chambers. These results were used to validate the numerical model before proceeding with further simulations. The results of this validation process for two different fluids is shown in Figs. A.1 and A.2.