Numerical effects on the simulation of surfactant flooding for enhanced oil recovery

Numerical simulation of surfactant flooding using conventional reservoir simulation models can lead to unreliable forecasts and bad decisions due to the appearance of numerical effects. The simulations give approximate solutions to systems of nonlinear partial differential equations describing the physical behavior of surfactant flooding by combining multiphase flow in porous media with surfactant transport. The approximations are made by discretization of time and space which can lead to spurious pulses or deviations in the model outcome. In this work, the black oil model was simulated using the decoupled implicit method for various conditions of reservoir scale models to investigate behavior in comparison with the analytical solution obtained from fractional flow theory. We investigated changes to cell size and time step as well as the properties of the surfactant and how it affects miscibility and flow. The main aim of this study was to understand pulse like behavior in the water bank, which we report for the first time, Our aim was to identify their cause and associated conditions. The pulses are induced by a sharp change in relative permeability as the interfacial tension changes. Pulses are diminished when adsorption is modeled, and ranged from 0.0002 kg/kg to 0.0005 kg/kg. The pulses are absent for high-resolution model of 5000 cells in x direction with a typical cell size as used in well-scale models. The growth or dampening of these pulses may vary from case to case but in this instance was a result of the combined impact of relative mobility, numerical dispersion, interfacial tension and miscibility. Oil recovery under the numerical problems reduced the performance of the flood, due to large amounts of pulses produced. Thus, it is important to improve existing models and use appropriate guidelines to stop oscillations and remove errors.


Introduction
An oil field development may be subdivided into three stages: appraisal, mature and late field development for oil recovery. Unfortunately, most of the world's oil fields have reached, or are approaching, the late field development phase that is a few years of the total recoverable field oil production before abandonment. At this point, secondary recovery has been exhausted and production may be considered as being no longer economically viable. Therefore, companies have explored alternative methods which are also referred to as tertiary, whereby oil is recovered by the injection of substances not naturally found in the petroleum reservoir. Chemical flooding is an important example of such a method. Chemical flooding techniques are often applied to alter fluid mobility and enhances the efficiency of the flood. This process is considered because conventional or secondary flood techniques have the potential to recover a third of the initial oil-in-place, leaving residual or by-passed oil in the reservoir. Therefore, it is crucial for companies to explore tools that are capable of modeling the physical phenomena as well as solving the resulting sharply changing fluid interfaces associated with such techniques to improve the field economic viability [1,2].
Experimental and analytical modeling approaches have been used to formulate laws and correlations for chemical flooding techniques, which are used to predict fluid flow and frontal advances through a porous medium. Numerical modeling approaches use a gridded format that can accommodate the reservoir description. This is integrated with nonlinear, multiphase partial differential flow equations which cannot be solved by an analytical approach, except in limiting circumstances. Numerical models have been widely used as an inexpensive alternative to experimental or physical modeling to assist with decision making through reservoir simulators. Simulations help predict and optimize recovery from petroleum reservoirs. Over the years, commercial simulators have utilized proven computational advancements to combine systems of mathematical equations to describe the process of enhanced oil recovery with surfactant flooding, which is the focus of this study. Numerical stability and accuracy of the surfactant flood model are very important to accurately compare recovery techniques for oil production. Today simulations are used on a daily basis, with a degree of dependability of the solution.
The simulators solve systems of nonlinear partial differential equations used to define the physical behaviors of surfactant flooding and recovery process. The simulations approximate the solutions by discretization of time and space which can lead to spurious oscillations, instabilities, or deviations in the model outcome [3]. The results may be used to predict or optimize oil production, so inaccuracies in the model could lead to false predictions, which can harm the economic valuation. Similar potential numerical issues may be encountered using the commercial simulator, which was documented in the simulator technical manual [4] with local grid refinements leading to numerical differences between the coarse and highly refined cases. Nevertheless, modeling and simulations are necessary tools in evaluating the practical application of surfactant flooding and potential implementation along with other chemicals.
Analytical methods have been used to compliment numerical simulation. Welge [5] built on fractional-flow theory and presented an analytical method for calculating average saturation and oil recovery performance by gas flooding in a onedimensional, three-component, two-phase system. Such analytical solutions can be used to study linear displacement for accuracy purposes. These approaches have been extended to waterfloods as well as chemical floods including polymer and surfactant methods [6]. Moreover, stability analysis can be used to ascertain if the numerical simulation result is reliable or requires further study to identify any underlying problems. Therefore, we define a solution as being unstable in the usual way such that a perturbation causes errors or inaccuracies that grow in time.
In a recent study by Paula et al. [7], instabilities of the physical process were reported for miscible displacement at high mobility ratios which the effects were shown to depend on dispersion and viscosity. The paper focused on viscous fingering patterns in two dimensional flow and discretization of the pressure equation to improve numerical solution, which differs to our approach. However, the existence of such complicated behavior requires that we are aware of potential numerical effects. Sadegh et al. [8] reported some discrepancies between the analytical and numerical results for surfactant compositional flooding which was caused by changes in the capillary number when partition coefficient of the CMG simulator was greater than unity. Moreover, a validation of the solution with a black oil simulator and high-resolution grid is important, as the study used a maximum of 100 grid-blocks which can be considered as a coarse model. This is the topic of our study. To control inconsistencies caused by numerical dispersion for miscible and near miscible flow, AlSofi and Blunt [9] suggested a physical dispersion-reduction scheme based on segregation-in-flow within a grid block as an alternative to implementing a higher-order discretization method. They proposed 50 grid blocks would be a sufficient alternative to predict oil production accurately. On the other hand, Connolly and Johns [10] used a high-resolution numerical scheme for miscible flow in 1D models to show the adverse effect of numerical dispersion, that caused mixing associated with permeability variation. The authors researched miscible flow but did not consider surfactant flood and its properties effect on the numerical scheme.
Mohammad et al. [11] formulated a black-oil model capable of modeling surfactant microemulsion phase behavior in the University of Texas Chemical Flood Simulator (UTCHEM) to improve the numerical accuracy of nonlinear chemical processes. The new model was validated using low-tension gas flooding with a surfactant-alternating gas process. This study ignored the potential for perturbations resulting in erroneous instabilities which should be considered for the black oil model.
In low salinity water flooding, the change between oil-wet (high salinity) and water-wet (low salinity) relative permeability curves have been reported to cause inconsistencies between the numerical and analytical solutions [12 -15] which took the form of pulses. These unphysical effects were observed in one-dimensional models of homogeneous reservoirs. In this study, we extend the method of Al-Ibadi et al. [13] to model a surfactant miscible flood. We also show how the magnitude of the pulses can be examined using the fractional-flow theory as a function of adsorption which was neglected in the paper.
In this paper, we discuss physical and numerical models of the miscible surfactant flood enhanced oil recovery (EOR) processes. We will examine the linear stability and the varying adsorption effect on the numerical scheme is also displayed. Then, we focus on the causes of perturbations, and report observations determined using the analytical approach.

Surfactant mechanism
Surfactant flooding has been studied in this paper and several mechanisms of recovery have been reported in the literature for over 40 years [16 -23]. The surfactant molecule has one end that is oleophilic and the other hydrophilic. The key mechanism therefore involves the fluid-fluid interface (oil and the displacing fluid), with the aim of reducing interfacial tension (IFT) and altering the capillary forces which influences the mobility of bypassed oil in the porous media. Oil and water are initially immiscible but the reduced IFT makes the fluids miscible. This helps displace the residual oil, normally trapped at the pore scale. These effects are expressed through the capillary number (N c ), which is the dimensionless ratio of viscous to capillary forces to improve sweep and displacement efficiency. The relationship between the capillary number and oilwater IFT is given as follows [23]: where v is the interstitial velocity of water, μ w is the aqueous phase viscosity and IFT ow is the oil-water interfacial tension also referred to as IFT.
In water flooding, the N c is usually 10 −9 to 10 −7 . The surfactant chemical interaction leads to a decrease in IFT and hence an increase in N c . If N c reaches 10 −3 then nearly 100% of oil initial in place may be recovered. This occurs when the surfactant reduces the oil-water IFT from typical levels of 20-50 mN/m, down to 10 −2 to 10 −3 mN/m.
The decrease in IFT has been shown to alter relative permeability curves in core floods. The fluids are immiscible for waterflooding but at infinite capillary number, they become miscible and the relative permeability curves are linear. Improved oil recovery is usually observed when the capillary number exceeds 10 −5.5 would not improve recovery [24]. These effects are captured in the simulation model.

Modeling of surfactant flooding
Changes to the relative permeability curve as a function of IFT and N c are modeled by interpolating from a set of input relative permeability curves for waterflooding under immiscible conditions to straight line curves for fully miscible flow [4]. In Fig. 1, the waterflood relative permeability curves are K ro (orange) for oil and K rw (blue) for water. At infinite N c , the straight-line curves are used (K rw_surf , gray and K ro_surf , yellow, for water and oil respectively). Since the relationship between surfactant concentration and IFT calculation was unveiled [25,26], surfactant transport models have been widely used in both analytical and numerical models for surfactant flooding to quantitatively analyze the mechanisms [6,27].
The surfactant concentration C o influences the change in oil-water IFT (Fig. 2), which causes a switch in the relative permeability curves. The immiscible set of relative permeabilities curves, denoted by K rw, and K ro , are applied in the absence of surfactant. With the influence of surfactant, the miscible relative permeabilities curves, denoted by K rw_surf , and K ro_surf are applied. Along with the switch from immiscible to miscible relative permeability curves, the surfactant contributes to the combined effect of the IFT, and aqueous phase viscosity. Firstly, the IFT is calculated using a function equivalent to that shown in Fig. 2. The IFT ow is initially at 0.05 N/m when there was no surfactant concentration in the solution. The critical micelle concentration (CMC) is achieved at the lowest IFT ows value with surfactant concentration of 1 kg/sm 3 . Once the CMC is achieved, the IFT remains constant with an increase in surfactant concentration. Meanwhile, the aqueous phase viscosity increases linearly with the surfactant concentration to the highest value of 5 cP.  The relative permeability curves reflect the miscibility as a function of capillary number interpolation is controlled by the parameter F kr . This term is calculated from the logarithm (base 10) of the capillary number and can be used to define an instantaneous switch between the set of relative permeability curves. However, this instantaneous switch may lead to convergence problems, slowing the calculations and even produces numerical effects in form of oscillations and pulses in the solution. Consequently, a lower threshold of F kr equivalent to zero, results in the immiscible flow while the upper threshold of F kr , equal to one, results in miscible displacement. In between these limits, the curves in Fig. 1 are interpolated using F kr as a linear weighting in relation to the two limits.
It is very important for the prediction of surfactant flooding that the fronts are resolved accurately, which can be studied through the fractional flow theory [28,29]. We illustrate the fractional flow function F w (S w , C, F Kr ) for certain combinations of surfactant concentration, C, and F Kr written as (also plotted in Fig. 3): where C o is the maximum concentration of surfactant and is equivalent to that of the injected water. F w (S w , 0, 0) is the fractional flow function for a normal waterflood, F w (S w , 1, 0) is the fractional flow as if the surfactant concentration has reached 1 kg/Sm 3 without any change from the immiscible waterflood relative permeability behavior. Therefore, F w (S w , 0, 0), F w (S w , 1, 0) and F w (S w , C o , 1) are the fractional flow functions given by The surfactant solution is modeled as a single solute dissolved in water and alters the viscosity, μ w , to stabilize the flood front and improve sweep efficiency. The influence of the surfactant concentration may be examined via the fractional flow curves (See Fig. 3). The value of the surfactant adsorption can be observed through the constant D i , which is achieved through the use of the Welge tangent on the fractional flow curve as seen in Fig. 3 [6].
Fractional flow theory can be used to represent flow in a 1-D homogenous reservoir, neglecting the effects of gravity, capillary pressure, and physical dispersion [29]. This can provide a mathematical prediction of the physical transport of fluids and concentration through the porous medium, including the chemical reaction on the rock surface. We acknowledge this analytical solution will be acceptable for coarsescale models in which the pressure difference from the viscous forces is usually much larger than from capillary forces, owing to high flow rates.
Researchers [27][28][29][30][31][32][33] have examined the dispersion of the surfactant solute front, involving numerical and physical dispersion at the core scale. However, owing to the fact that physical dispersion is ever so often offset by numerical dispersion; the first is ignored. We use the Peclet number to characterize dispersion.
The value of N pe ≪ 1 indicates dispersion dominates and advection is negligible. When N pe > 10 advection dominates on the other hand and spreading is almost non-existent. If N pe~1 then neither dominates over the other and no approximation to the equation can be justified.
D in Eq. 6 is given by:  We assume molecular diffusion to be insignificant in comparison with the dispersion. The total dispersivity α total represents a combination of physical and numerical dispersion [29].
Where α Physical is derived from continuous experimentally or field test measurement and to be zero for this research. While α Numerical can be estimated by: For grid block sizes Δx, time step Δt and ∅ is porosity and the arithmetic subtraction operator (−) applied for Implicit Pressure-Explicit Saturation (IMPES) procedure and addition operator (+) for the fully implicit procedure [30]. The decoupled implicit scheme was used in this study such that physical dispersion was neglected. Eq. 10 can be further simplified for the fully implicit Peclet number to obtain [34]: N pe is used to evaluate numerical problems that may arise from the transport calculations due to overly large time steps and cell sizes. The numerical effects lead to smeared spatial gradients of saturation and concentrations as well as grid orientation effects. The problem results from the truncation error in calculating the movement of saturation fronts as if additional physical dispersion were present. We will also show that other effects include pulses in fractional flow and saturations. Unfortunately, in explicit schemes, instabilities are apparent when the combination of Δx and Δt violates the numerical stability criterion [30]. These instabilities can lead to behavior not unlike that which we will report. However, it is not our contention that we are observing instabilities.
In the literature [30] for a sufficiently small time step in which the level of numerical dispersion is due to the spatial term, the level of dispersivity inbuilt in the simulator ought to be about half of the grid block size. The Courant-Friedrichs-Levy found instability issues due to the time step effect, and to mitigate such problems, the Courant condition was derived to suggest a formulation for a time step size of each grid cell at a known fluid velocity and cell size. We would like to highlight that this condition was derived for an explicit scheme whereas this study used a decoupled implicit scheme that appears to have some numerical issues. To examine the stability of the numerical solution in the surfactant flood model, we assume Eq. 11 as the limit, and to ensure our model is an advection dominated system. The time step limit can be given as [32]: where the v is the advection velocity. As mentioned earlier, the key mechanism of the surfactant is alteration of interfacial tension, influenced by the level of adsorption. The adsorption level makes the process less efficient and also alters the Welge tangent points across the fractional flow curve, which influences change in the fluid flow behavior and directly affects the velocities of surfactant induced and formation waterfronts. This is incorporated in the numerical modeling approach to impact the changes in the concentration distribution and mobilization of residual oil in the reservoir. Therefore, it is important to study and test the role of adsorption on the sharpening of the front and numerical solutions.
Modeling the impact of adsorption on the surfactant flood is captured by using a Langmuir-type isotherm. Adsorption is linear for small values of surfactant concentration but approaches a constant at high values. The analytical adsorption function allows for dependencies of adsorption on the rock permeability, salinity, and surfactant concentration [4]: Where C ads is the rock adsorbed concentration, C is the concentration in the solution surrounding the rock, m is an exponent for concentration dependence, C se is the effective salinity for surfactant concentration, K is the grid block permeability, K ref is the reference permeability, n is an exponent for permeability dependence, while a 1 , a 2 , b are adsorption coefficient input parameters for the simulation. The surfactant model range of adsorption can be evaluated using the retardation equation: where R ci is the retardation factor, v is the advection velocity of water and v ci is the concentration velocity. Using the equations above we could evaluate a relationship between adsorption and numerical effects. We can use the analytical solution as a baseline to study adsorption effects and estimate the numerical effects on 2D and 3D models.
An algorithm for numerical stability analysis, based on the perturbation method, was developed by Druetta et al. [3] to ensure convergence in a two-phase, three-component system. The approach is appropriate for two-dimensional models; however, one must account for the impact of the surfactant mechanisms such as how adsorption may influence the precision of the solution.
Using numerical analysis, Al-Ibadi et al., [12][13][14] reported the propagation of pulses within the cells as numerical errors and highlighted the importance of precision for low salinity waterflooding. A similar approach is used in those models to switch relative permeability curves as salinity is reduced. This can induce a sharp increase of the flow of oil from a cell with reduction in water and creates a pulse. Simulation of surfactant flooding is dependent on the logarithm of the capillary number, influenced by the ratio of viscous to capillary forces and a sharp change to the relative permeability can also be induced.
Nevertheless, experiments have demonstrated that the synthesis of low surfactant concentration with new chemicals can achieve ultra-low IFT values and acceptable adsorption levels for EOR. Therefore, it is necessary to study how this type of flow might cause numerical problems in form of pulses that influence incorrect results. The overall goal is to present design consideration for surfactants aimed to attain low interfacial tension at low surfactant concentration, as well as low adsorption levels on the rock surface to recover bypassed oil saturation [15,35].

Description of the simulation model
We began with a simple homogeneous 1D model solved through the black oil commercial simulator [4] to evaluate the flow behavior of the surfactant flood. The dimensions of the reservoir setup are 500 m long, 500 m wide, and 1.89 m thick. The model was chosen to provide the simplest possible illustration of surfactant flooding, to streamline the interpretation of the flow behavior in a reservoir. One injection well was used with a production well at the extreme ends of the grid block in the x-axis, and both are completed vertically and fully perforated in all the interconnecting grid cells.
The base case model was represented with 100 cells, each defined in the x-direction. Models were tested with a range of time steps to consider the influence of numerical dispersion based on the time step alone. We further compared the impact of grid refinement on the numerical solution profile against the analytical models which predict the shock front behavior [29,30]. The displacing fluid consisted of an aqueous phase, a mixture of water and the surfactant solution, where the injected surfactant increased the viscosity of the displacing fluid and reduced IFT between the oil and aqueous phases. Some features of the reservoir characteristics for the surfactant flood model were considered constant and due to the large dimensions of the model. Gravity, and capillary pressure effects were considered negligible compared to the viscous forces. These assumptions simplified the model and removed their impact on the numerical dispersion. The model used isotropic permeability to achieve the desired pressure limit. A commercial simulator [4] was used to solve the aqueous phase solutions simultaneously with surfactant solute transport and all related equations using the decoupled implicit finite difference scheme.
The simulation models were studied using the rock and fluid properties along with the reservoir conditions as tabulated in Table 1. To focus on the surfactant flood solution, the water, and surfactant slug was injected from the beginning as a secondary process, through the injection well sweeping the mobilized oil to the production well. The injection rate was constant at 400 rm 3 /day. Time step and grid size were our primary objective to evaluate the appearance of numerical effects and impact on solution accuracy.

Simulations of surfactant flooding
In this work, we evaluated the stability, precision, and efficiency of the numerical model and the components that influence the flow of surfactant flooding. The precision is termed as the minimal deviation of the solution affected by the numerical problem when compared to the exact partial differential equation that represents flow. Meanwhile, instability occurs when numerical effects result in a divergence of the numerical solution as a function of small perturbations. The numerical effects, also known as numerical dispersion, may also result in pulses and convergence problems [12][13][14]. They may be imprecise but stable, nonetheless. The efficiency of the numerical model refers to the simulator's ability to solve a wide range of physically distinct problems. The computational efficiency of the solution scheme accounts for the number of numerical values computed after the required "Newton iterations" per time step and the time it took to solve the linearized system on each iteration. Firstly, we evaluated the numerical solution of four models which consist of coarse (10 × 1 × 1), base case (100 × 1 × 1), refined (1000 × 1 × 1) and high-resolution (5000 × 1 × 1) block setup with time step defined according to the Courant condition, and compared to the solution obtained by fractional flow theory for surfactant flooding. The advection velocities of the surfactant-induced and formation waterfronts are 1.9 m/ day and 2.7 m/day, respectively. We obtained solutions for various time steps and obtained the Numerical dispersion and Peclet number (See Table 2).
The example showed the appearance of numerical effects in the form of pulses ( Fig. 4 and 5) and a shifting of the surfactant concentration (Fig. 6) and viscosity behavior (Fig. 7). The high-resolution solution (N pe = 2000) tends to follow that of the fractional flow theory as the decrease in the cell size to 0.1 m and time step of 0.04 days produced a sharp front in the water cut (See Fig. 4). The injected waterfront traveled at S w = 1 and f w = 1. The formation waterfront moved ahead of that at lower saturation and fractional flow values given by the Welge tangent which passes through the blue curve in Fig. 3 to -D i on the saturation axis.
An increase of the cell size to 0.5 m and time step to 0.18 days in the refined solution (N pe = 400), introduced numerical effects in the form of small pulses in the plateau (oil bank) with low amplitude reducing the precision of the result compared to the water cut solution obtained from fractional flow theory (Fig. 5). The base case (N pe = 40) with a cell size of 5 m and a time step of 1.82 days was observed to show higher pulse amplitude and longer wavelength, which caused an even greater deviation in the water cut compared to the refined and high-resolution water cut profile. A further increment of the cell size to 50 m and time step to 18.24 days as observed in the coarse model (N pe = 4) is shown in Fig. 5. The flow behavior tended to lose numerical precision due to the pulses which had significantly higher amplitude on the formation water plateau and failed to match the analytical solution. The coarse model also produced formation water breakthrough slightly earlier and the surfactant showed a more dispersed front due to larger cells.
The displacement profile of the surfactant concentration ( Fig. 6) was obtained numerically from the cell closest to the production well, with a comparison of the corresponding effect on the effective aqueous phase viscosity (See Fig. 7) and  the oil-water IFT (See Fig. 8). The analytical solution of the surfactant concentration was very similar to that obtained from refined (N pe = 400) and high-resolution (N pe = 2000) but with a slight deviation in the base case (N pe = 40). We observed that the numerical solution for the coarsest model was highly imprecise due to the combined effect of oil-water IFT and solution viscosity combining with numerical dispersion from the discretization errors. This effect may be seen through the earlier IFT change at 0.8 PVI for the coarse model compared to other models at approximately 1.1 PVI (See Fig.  8). Increasing the cell size had a large effect on dispersion of surfactant. However, this mixed within the formation water in this case. While viscosity increased and the miscibility changed, the general speed of the injected waterfront was largely unaffected by comparison although the second waterfront did rise more slowly.
To show how the variation in effective surfactant concentration may lead to unfavorable flow, the coarse model at 1.5 PVI could be seen to reach viscosity of 5 cP with far lower concentration of 27kg/sm 3 in comparison to the accurate input value of 30kg/sm 3 as achieved in the base case, refined and high-resolution models. The effective surfactant concentration also triggered the change in miscibility, which is a function of the logarithm of the capillary number. Therefore, imprecise values of the effective surfactant concentration will lead to false values of the oil and water relative permeabilities, impacting the precision and stability of the model.
The relative advance of the surfactant profile as the cells were increased in size depends on the balance between dispersion and adsorption. With less adsorption we anticipate that it would advance faster. If it moves as fast as the head of the formation waterfront then more complex behavior could be seen, particularly if there is little adsorption or more dispersion. Besides precision, any significant improvement in the timestep domain would require an adaptive step solver, which would involve a larger computational time. The precision must be weighed against computational efficiency. The error in the coarsest model renders it useless as a tool for predicting behavior. The origin of this numerical effect on the fluid flow behavior is further evaluated through the mobility analysis in the grid cells and parameters influencing the numerical computation of the surfactant flood model.

Numerical computation of the mobility and weighting function
We examined the origin of the numerical effects from the coarse model (which consisted of 10 cells). We evaluated the dynamic behavior of relative permeability, saturation, log(N c ) and also the fractional flow (into the downstream cell) on a cell by cell basis (Fig. 9). The simplified fractional flow was obtained by the ratio of This analysis follows that presented in [13,14] for low salinity waterflooding with a similar hypothesis water flow rate to the total flow rate with a value between 0 and 1. The relative permeability of water and oil and the fractional flow ( Fig.9) was observed to have pulses which grew progressively as each cell transitioned from immiscible to miscible behavior and also as the viscosity changed. Similar effects have been reported previously albeit for low salinity waterflooding [12][13][14] where the change in relative permeability (wettability) as a function of salinity appeared to generate equivalent pulses. We adopt a similar analysis and interpretation here.
Water saturation in the first cell (i = 1) increased quite quickly reaching a maximum value at 0.2 PVI (Fig. 9a). The surfactant concentration also increased (Fig. 9b) and IFT was reduced. There was an initial change in K rw (Fig. 9c) and K ro (Fig. 9d) at around 0.02 PVI when we see the log(N c ) passed through the transition, switching flow from immiscible to miscible behavior (Fig. 9e). The change to fractional flow (Fig. 9f) occurred slightly later as the viscosity in the cell changed.
The rise in saturation in the second cell occurred slightly later as the front moved with finite speed. Formation water was mobilized ahead of the front. There is a slight "kink" in the water saturation at around 0.1 PVI. Numerical dispersion caused the surfactant concentration to start to spread out. We see that the relative permeability to water shows a slight plateau around this time. This was coincidental with the rise in log(N c ) as IFT dropped and flow behavior became immiscible. The viscosity then rose (as it did in i = 1) creating a broader plateau in the fractional flow.
This process continued as the fluid fronts moved through the reservoir. The formation waterfront moved faster than the later high viscosity surfactant front creating a plateau in saturation which spread out. Each cell saw a minor increase in log(N c ) ahead of the main surfactant front, which was caused by formation water before a sharp switch in as the fluids become miscible. The initial rise in log(N c ) was due to the mobilization of the formation water and a change to the interstitial velocity. The water relative permeability jumped up with the switch to miscibility as did oil relative permeability jumped up. Then the water phase became less mobile and relative permeability declined somewhat. Fractional flow was largely unaffected by the switch in miscibility as the increase in oil and water relative permeabilities seemed to cancel each other out.
Qualitatively similar behavior was observed in the base case (100 cell) model in cells 91 to 99 (See Fig. 10a and  Fig. 10c). The cells are closer together so the changes in behavior occur more frequently and we magnify the PVI on the x-axis of the plots. The initial short rise in log(N c ) was observed in each cell followed by a sharp change as fluids became miscible. Then after that there was a change in relative permeabilities and fractional flow as viscosity changed. In this case, we see pulses in saturation and relative permeabilities as well as fractional flow. On top of that there seems to be a dual frequency of pulses. In cell 99, for example, there were some high frequency pulses just before the cell became miscible by the arrival of the surfactant. These high frequency pulses were also seen and correlated with behavior in the preceding cells. The high frequency pulses appeared to correlate with the change in log(N c ) in Fig. 10c.
As mentioned above, Al-Ibadi et al. [12][13][14] have shown similar such pulses for low salinity waterflooding. In that case relative permeability is changed rapidly due to the reduction of salinity as injected water arrives. This process emulates a change in wettability. The change of relative permeability increases oil flow and reduces water flow. The resulting pulse in a cell can then propagate downstream and appear in saturation and fractional flow. The pulses were shown to grow in magnitude as the low salinity front arrived. In the numerical case, dispersion slowed down the increase of water saturation that is normally associated with a sharp low salinity water front in the analytical solution.
In the latter case, saturation and salinity change together preventing the pulses. We hypothesise that the same process occurs here for surfactant flooding as oil becomes more mobile when miscibility changes. Water mobility was initially increased by the change of relative permeability but then decreased when the viscosity was reduced with increasing surfactant concentration. We observe the effects in the 100 cell model and in the water cut of the 10 cell model by the presence of the double bump before the surfactant waterfront arrived. Al-Ibadi et al. [12][13][14] showed that the pulses grew in number as the formation and low salinity waterfronts separated out. We observe the same here in longer models.
In earlier time, longer wavelength pulses were also observed in the fractional flow curve (but not in the water relative permeability curve). These appear as troughs which become deeper in time and have a periodicity of about 0.08 PVI (roughly the time for 10 cells to change from immiscible to miscible upstream). The troughs also correlate across cells. This effect was not reported for low salinity waterflooding. The grid refinement reduced the effect of numerical dispersion but there remain some numerical effects. We can conclude that the pulses were a result of the presence of multiple phases at a given point, where the upstream weighting of the relative permeability values produced unphysical pulses causing a ripple effect that was transmitted through the cells to the production well. Grid refinement actually introduced more pulses, but they were compressed and began to converge to the analytical solution. Al-Ibadi et al. [12][13][14] showed that the number of ripples increased when cells were considered further from the injector because the formation and injected waterfronts separate out. Refinement also saw pulses diminish in magnitude.

Time step reduction and effect on the appearance of numerical effects
We examined the impact of time step reduction on the surfactant flood base case model concerning the appearance of numerical effect and computational efficiency. The base case grid block setup (100 × 1 × 1) was simulated with various time steps, corresponding Peclet number, and computational time (see Table 3).
We show in Fig. 11 that the time step contributes heavily to the numerical dispersion. The reference case (N pe = 40) produced the noisiest, and with significant pulses spanning across the plateau in comparison with the other case studies.  Fig. 12).
The computational efficiency of the solution profiles also was considered. The reference time step was seen to lack computational efficiency when compared with a time step of 1 day (N pe = 52) and 0.45 days (N pe = 64) (See Table 3). We did observe convergence issues with the base case. The increase of numerical dispersion seems to have exacerbated the issues, hence additional iterations were required as the problem became harder to solve.
Closer inspection of Fig. 12 reveals that the apparent long wavelength pulses were clearly an artifact of the dispersion as they were not present in any of the models with N pe > 40. High frequency pulses were observed as the injected waterfront arrived, but the plateau of the formation waterfront was mostly flat and equivalent to the analytical solution. An exception to this was the case of the time step of 0.2 days (N pe = 72) which produced a solution profile with pulses but there was also a jump in water cut above the value of the analytical solution. Further investigation was required to underpin why the initial bump appeared in the solution profiles that were contrary to the analytical solution. We can conclude that to enhance the precision for the numerical modeling of the transport of surfactant flood weighed against the appearance of nonphysical solution and computational efficiency, the Courant condition considered should be as follows: 8 The effects of adsorption: Impact of zero adsorption on oil recovery We also examined how the variation of adsorption in the surfactant flood model affected recovery factor as a function of the different grid cell sizes. We ran the base case model that was used above with zero adsorption (Fig. 13). The analytical result was almost perfectly matched by the most refined model (N pe = 2000) and oscillations appeared as the model was coarsened. We observed the appearance of an obvious dip in the water cut solution profile just ahead of the surfactant induced waterfront for N pe = 40, the pulse had a different magnitude compared to the base case, which modeled adsorption. As described above, the dispersion of the surfactant initially changed the mobility of the water which changed relative permeability and caused the dip. Then the IFT decreased when the threshold of the capillary number was reached. The dip was largest for the coarsest grid which also saw a much earlier breakthrough of the formation waterfront while other cases correctly predicted this. Overall the main effect of enlarging the cell size was to change the water cut and to predict earlier breakthrough of surfactant. The dispersed surfactant moved faster than the formation waterfront in the absence of adsorption (Fig. 13). This disturbed the formation waterfront which no longer moved at the same steps in saturation and fractional  flow. Fig. 14 shows that, with the exception of the coarsest model, the cumulative volume of oil produced matches the analytical solution very well. The deviations are small compared to the cumulative volumes. We conclude that in the absence of adsorption, a more refined model is required to obtain sufficient precision.

Varying adsorptions
Lastly, we examined the role of the adsorption on the sharpening of the front and appearance of pulses in the numerical solutions. We used the base case model and made some modifications to the adsorption values for this study (See Table 4).
The increase of the adsorption rates resulted in the acceleration of the formation waterfront and a retardation of the injected waterfront. Adsorption removed surfactant from the injected water so that it behaved like formation water. We can observe the impact of adsorption on the interpreted D i achieved through the use of the Welge tangent on the fractional flow curve (See Fig. 3). We also saw a change in the numerical effects.
It was observed in the comparison between the water cut profile with zero adsorption through the appearance of an obvious dip (Fig. 15) that can be directly related to the pulse observed in the water relative permeability curves (Fig. 16), which was caused by the change of viscosity as the surfactant concentration rose to 1 kg Sm −3 . The relative permeability of water appeared to drop as the surfactant arrived and then increased as the oil was displaced. The dip in water cut gradually diminished when adsorption was increased to 0.0003 kg/kg; this was then followed by the recurrence of the pulses as the adsorption was further increased to 0.0005 kg/kg. The interpretation of the second pulse on the solution profile may be caused by the appearance of the numerical effect observed during the mobilization of the formation waterfront (Fig. 17), which was initially less mobile. As the adsorption rate was increased further, the surfactant concentration decreased. The decrease in the concentration lagged in the water solution displacement due to the pore volume that can be accessed by the surfactant; this caused the surfactant to move slowly and thus the initially immobile oil was released and recovered more quickly.
We observed from the surfactant flood fractional flow curves from the simulator data which were computed as the producer cell first saw formation water invade then a change from the immiscible curves towards the surfactant induced miscible curves. The input fractional flow curves are shown as reference (Fig. 18). The waterflood initially appeared to correspond to the predicted fractional flow theory, while there also appeared to be complexity all along the curve prior to the transition to the miscible surfactant flood curve. This might be the reason for the appearance of the dip observed in the zero adsorption and 0.0001 kg/kg value; which showed a jump and deviation from the waterflood curve in value of saturation and fractional flow due to the arrival of the surfactant concentration. From the theoretical analysis, the Welge tangent was used to understand the water saturation shock fronts of both the water formation and surfactant induced fronts. The two cases (zero adsorption and 0.0001 kg/kg) appeared to have two possible tangent points, one along the waterflood curve and the deviated point that corresponds to the tangent point of the adsorption value of 0.0003 kg/kg with less match of the water saturation. The dispersion had a clear impact on changing the shape of the curves making them less precise.

Discussion
In this paper, we examined the numerical errors that can arise during the simulation of miscible surfactant flooding in coarse grid models. We used a 1D homogenous model at the interwell scale and found significant deviations when compared with analytical and high-resolution solutions. Grid and time step refinement studies were used to access the impact of the numerical error on their solution which is a common method used by many researchers for model selection [12,28,30,32,34]. On the other hand, based on the results of this study, we suggest a modified Courant condition for maximum timestep which is distinct from the previous reports to address the numerical instability. This approach for timestep constraint may be specific for the application of the decoupled implicit scheme that is available in the commercial reservoir simulator that we used. This method has been reported to give stability at the expense of accuracy and computational efficiency [4]. Nevertheless, we observed that the errors occurred in the surfactant concentration front at dispersion-dominated flow, which was particularly observed when the Peclet number was less than 40. Other researchers have shown that for low salinity water flooding, no instabilities were reported for Peclet number less than 100 [12]. The relationship between dispersion and propagation of the pulses helps distinguish and eliminate errors. Some features that contribute to the mobility and weighting in the simulation were likely to contribute to instability. Key mechanisms influencing the surfactant flood model may be used to further understand the origin and the spread of the pulses. The switch in relative permeability seems to have been important as has been reported for low salinity waterflooding by Al-Ibadi et al. [12][13][14]. These effects are important to understand, particularly if experiencing instabilities such as those reported by Paula et al. [7]. On the other hand, Abbas et al. [36] reported the importance of numerical analysis of surfactant adsorption to understand the flow behavior using CMG-STARS software. The authors carried out adsorption tests to accurately develop new chemical flood processes with improved oil recovery. Our approach showed how the adsorption/retardation factor observed from the fractional-flow theory can be used to understand and compare the perturbation associated with the coarse model.  Fig. 15 The effect of adsorption sensitivity on the water cut profile from when adsorption was neglected (orange line) to 0.0008 kg/kg (red line)  Fig. 16 The water relative permeability for the adsorption sensitivity versus pore volume injected Furthermore, our approach can be used as an option to examine the magnitude of the pulse for accuracy purposes. The numerical solution of the base case model with the modified Courant condition tested according to various adsorption values produced two separate types of numerical effects. The first occurred when dispersion led to an earlier change of viscosity and miscibility, triggered by the lack of adsorption in the surfactant flood model. Whereas the second effect was formed by the mobilization of the formation water that predominantly caused the solution inaccuracies.
The finding of this study would help in better understanding the numerical effects of surfactant miscible flood, maximizing grid and time step by selecting effective limits, interblock inspection and then considering the analytical solution with adsorption tests to categorize the pulses. Applications of the analytical solution and varying adsorption parameters could be useful to modify the fractional-flow model and upscaling of surfactant flood as a function of the effect of dispersion, miscibility, and solute concentration. In further work, we will use the new analysis to address the key question of how to recognize and reduce these errors for a layered and heterogeneous reservoir for more accurate models. The study showed that particular attention is to be paid to the properties of the miscible flow and that due to the appearance of nonphysical numerical artifacts, the simulations are not certain to correspond to the physical models. Therefore, rigorous studies may be required for both the numerical and physical models to resolve the perturbation and accurately predict the performance of the surfactant flood on oil production.

Conclusions
We conclude the following from this study: & Numerical accuracy problems were identified for models with a typical cell size as used in well-scale models. & Pulses were observed in the oil bank with multiple frequencies or duration. & A dip in water saturation and fractional flow was observed for low levels of adsorption so that more refined models were needed. & The Courant condition should be observed as a minimum to avoid these effects.
Also, we suggest further work is required to examine the relevance of these findings in the 2D and 3D layered and heterogeneous reservoir scale representation of the surfactant flood model. Acknowledgments The financial support from Nigerian Petroleum Technology Development Fund is greatly appreciated. We would like to thank Schlumberger for the application of Eclipse 100 simulator. We also thank Hasan Al-Ibadi for his contributions to discussions on pulses as seen previously in low salinity waterflooding.
Code availability Eclipse 100 simulator.
Author contribution Not applicable.
Data availability Not applicable.

Declarations
Ethics approval Not applicable.
Consent to participate Not applicable.
Consent for publication Not applicable.

Conflict of interest Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ .