Pore-Scale Evaluation of Physicochemical Interactions by Engineered Water Injections

Engineered water injections have gained a lot of interest as an economic and effective method of improving the oil recovery. However, the complexity of the physicochemical interactions between the brines of various compositions, oil and rock has led researchers to provide multiple ways to explain this phenomenon. In this work, we evaluate the previously suggested mechanisms, namely wettability alteration and emulsification, against high-resolution micro-CT coreflood observations in a limestone sample. This is achieved by integrating the effects of above-mentioned mechanisms into a volume-of-fluid simulation by using geochemical modelling and experimental measurements. This has allowed us to explain the effect of capillary force affecting mechanisms, whereby we were able to achieve 6% increase in recovery factor. We have also observed that these mechanisms have limitation in improving recovery due to fingering and subsequent formation of the stagnation zones inside the core samples. When viscous effect is considered in numerical study, 22% increase in recovery is achieved by reorientation of the main flow paths and mobilisation of the previously unconnected oil clusters. This result is closer to 24% increase in recovery factor which was observed in experimental study and signifies that viscosity increase due to emulsification is an important mechanism of engineered water injections.


Introduction
Engineered water injections have become a much-discussed topic in the petroleum industry as an economic and sustainable method of increasing the oil recovery (Meng et al. 2015). Multiple reviews (Al-Shalabi and Sepehrnoori 2017; Awolayo et al. 2018;Hao et al. 2019) suggest that recovery can be enhanced by modifying both the overall salinity of the seawater and the composition of the potential determining ions (PDI) individually. For carbonate rocks, these PDI include Ca 2+ , Mg 2+ , SO 4 2− (Ding et al. 2019). Based on the spontaneous imbibition and sessile drop contact angle studies, researchers have shown that wettability alteration by low-salinity or ion-tuned waters is a mechanism that significantly contributes to increase in oil recovery in ex situ and in situ studies (Selem et al. 2021;Shirazi et al. 2020;Zhang et al. 2007). However, as the number of studies on the topic grew , inconsistencies in the evaluation of wettability change and oil recovery have been identified. This led the scientific community to develop additional methodologies to evaluate other potential mechanisms contributing to the efficiency of engineered waterfloods. Specifically, Masalmeh et al. (2019) have shown that the formation of in situ emulsion at the interface of oil and water is dependent on the oil composition and salinity of the injected water. Other studies Tetteh et al. 2021;Unsal et al. 2019) have shown that the volume of generated emulsion is also controlled by the ionic composition of the water and can be observed and quantified using microtomography. Therefore, in situ formation of emulsion is another mechanism which provides an additional effect of engineered waterfloods. This is due primarily to the fact that emulsions increase the viscosity and interfacial tension at the interface as well as redirect water by blocking large pores, therefore reducing the fingering effects (McAuliffe 1973a;Unsal et al. 2016).
These foregoing complexities have made it difficult to model the governing mechanisms and predict the recovery. Majority of early modelling studies on the topic were focused on wettability change with two distinct geochemical approaches. The first approach is based on Bond Product Sum (BPS) (Brady et al. 2015), which is a summation of products of surface concentrations of opposite ions on oil and mineral. If all the species are negatively charged, BPS would be zero, species would not be electrostatically attracted to one another, and therefore, the surface is considered strong water wet. If both positively and negatively charged surface species are present, BPS would be high, and the wettability will move towards more oil-wet. The second approach, which has received a greater attention, is DLVO (Derjaguin-Landau-Verwey-Overbeek theory of colloidal stability)-based calculation (Derjaguin et al. 1987). In this approach, the attractive-repulsive microscopic forces between a surface and colloid (in this context rock surface and oil particle) are described by a disjoining pressure. Hirasaki (1991) has provided a significant extension to this theory in the context of oil recovery by developing a model for contact angle based on the disjoining pressure. As such, this approach allows to precisely predict specific contact angle values by knowing the chemistry of the injected brine, oil and rock mineralogy (Sanaei et al. 2019). At this moment, there is a lack of clear approach for modelling the fluid-fluid interactions with exception of some early works considering the effect in simple geometries with no experimental validation (Alizadeh et al. 2021).
There are two main direct numerical (DNS) approaches that are frequently applied to simulate multiphase flow on a pore scale: volume of fluid (VOF) and a number of multiphase extensions of lattice Boltzmann method (LBM) (e.g. Rothman-Keller, Shan-Chen, free energy models) (Huang et al. 2007;Lanetc et al. 2022;Vogel et al. 2005). In LBM, such simulation parameters such as viscosity, adhesion and cohesion forces depend on a voxel size and time-step values (Krüger et al. 2017). Thus, the usage of realistic parameters may lead to unfeasibly large number of time steps or computationally expensive mesh refinement. In turn, flow parameters can be set arbitrary for the VOF method which is based on the Navier-Stokes equation (Deshpande et al. 2012). However, the employability of the VOF method is still restricted by high computational costs when a large 3D geometry is considered. This issue can be addressed by the employment of relevant image-based modelling (IBM) tools including grid coarsening and by parallel computing Raeini et al. 2014). Shams et al. and Frank et al. have shown that it is possible to precisely predict localised contact angles and capillary snapping within the scope of abovediscussed limitations (Frank et al. 2018;Shams et al. 2021).
The engineered water injections have been extensively studied by both lattice Boltzmann (LB) and VOF methods (Namaee-Ghasemi et al. 2021). For instance, the wettability alteration due to the injection of low salinity water in 3D micro-CT images of sandstones was modelled by (Akai et al. 2020) with the usage of colour gradient LBM. Likewise, (Aziz et al. 2019) used a VOF method to analyse flow dynamics under wettability alteration on 2D synthetic images. Nevertheless, both studies only considered the wettability alteration by the linear functional dependence between the contact angle and ions local concentration while other mechanisms (e.g. interfacial tension change, viscoelastic effects and emulsion generation) were neglected.
The goal of this study is to provide an exhaustive modelling approach for numerically predicting the results of the engineered waterfloods. To achieve this, we combine the observations from pore-scale imaging of the ion-tuned injections with numerical models of individual mechanisms and then evaluate the significance of both rock-fluid and fluid interactions on oil recovery from carbonate samples. To the best of our knowledge, this is the first work that considers the modelling of multiple mechanisms and at the same time compares the results of these models with pore-scale experimental data.

Rock and Fluids
Indiana Limestone was used as a representative carbonate sample, which is a well-characterised rock with relatively homogeneous mineralogy. A micro-plug with diameter of 5.9 mm was prepared to allow for high resolution of the CT image. Full list of its properties is provided in Table 1.
Three types of brines were used in this study. Firstly, a formation water (FW) represents a carbonate reservoir formation brine. For the secondary recovery, synthetic sea water (SW) was used (Ding and Rahman 2018;Yousef et al. 2011). For the tertiary recovery, ion-tuned water (ITW) was injected. Ionic modification was done by decreasing the Ca 2+ concentration 4 times when compared to SW. Full list of brine properties is presented in Table 2. 1 3 Dead crude oil was used in this study with its original properties listed in Table 3. For the purpose of CT imaging, the oil was mixed at 1:1 proportion with 1-Iododecane. To assure the oil-wetting conditions after the ageing, the total acid number (TAN) of the mixture was adjusted to 1.5 mg KOH/g using the stearic acid (Song et al. 2020;Standnes and Austad 2000).

Sample Preparation and Waterflooding
Samples were first saturated for 10 days with FW and then injected with the oil mixture. Next, they were aged for 2 weeks at 80 °C and 200 psi to establish the initial oil-wet conditions. Finally, the 20 PV of SW was injected at secondary recovery stage, followed by 20 PV of ITW. All the injections were conducted at 0.5 cc/min flowrate (Shapoval et al. 2021). This flowrate was estimated based on the scaling coefficient proposed by Rapoport and Leas, (1953) to minimise capillary end effects: where L is the core length (cm), ν is the flow velocity (cc/min), μ w is the water viscosity (cP).

Image Processing
Images were acquired at the following stages: dry, after ageing, after seawater injection and after ion-tuned water injection. Images were first registered to each other and then segmented using the watershed algorithm . At last, contact angles were (1) S c = L * ν * μ w > 1

Contact Angle Calculation
Rock-fluid interactions, quantitatively expressed as contact angle, are governed by the surface forces at the water-rock and water-oil interfaces. The summation of these surface forces gives rise to disjoining pressure (see Eq. (2)), which governs the ability of rock surface to adsorb or desorb oil.
where Π-disjoining pressure, Π VdW -van der Waals force, Π EDL -electrical double-layer force, Π S -structural force, h-water film thickness. Hirasaki (1991) provided an analytical solution, which allows to use the calculated disjoining pressure for the prediction of the contact angle: where θ-contact angle, σ-interfacial tension (IFT), h 0 -is a limiting film thickness.
(2) It was previously discussed that disjoining pressure curve for calcitic rocks shows both positive and negative values (Jalili and Tabrizy 2014;Khurshid and Al-Shalabi 2022;Mahani et al. 2017;Maskari et al. 2020). At our conditions, disjoining pressure was estimated to be negative (see validation in Supporting Information), an assumption proposed by Hirasaki (1991), which states that the limiting film thickness is a thickness of a water molecule monolayer (0.3 nm), was used.

Surface Complexation Modelling
Out of the three forces that makeup the disjoining pressure, the electrical double-layer (EDL) force is the one that is controlled by the chemistry of the brine. Specifically, this force is based on the electrochemical interactions on the water-oil and water-rock interfaces. Numerically, EDL force is calculated as a linear solution of Poisson-Boltzmann equation (PBS) at a specific boundary condition. There are various conditions proposed in literature for fitting the experimental data; however, two mostly used methods include constant charge (CC) and constant potential (CP). In this work, we have used the constant potential boundary condition, which has been previously used to successfully model the EDL force for carbonates (Bordeaux-Rego et al. 2021;Sanaei et al. 2019): where n b -the ion density of the bulk solution, k B -Boltzmann constant, 1.38*10 -23 J/K, T-Temperature, 298 K, ri -reduced surface potential (see Eq. (5)), k-reciprocal Debye-Huckel length.
where e is elementary charge (1.6 × 10 -19 C), i is the zeta potential.
As such, zeta potential is a parameter, which is required for estimating the contact angle and which is dependent on the ionic composition of the brine. Zeta potential can either be measured experimentally or modelled numerically using surface complexation modelling. Because the goal of this work is to maximise the modelling capabilities, we have used the surface complexation modelling (SCM) to predict the zeta potentials for various brines. Specifically, two SCMs were used to model the reactions between water and calcite and between water and oil, as described in Table 4.
For calcite, density of surface binding sides is assumed to be 4.95 sites/nm 2 ; specific surface area is assumed to be 1 m 2 /g (Ding and Rahman 2018). For oil phase, site densities of carboxylate and nitrogen base groups were calculated using the following relationships: where 0.602 × 10 6 is a conversion from mol/m 2 to sites/m 2 , a oil is the specific surface area of oil (3.7 m 2 /g), MW KOH is the molecular weight of potassium hydroxide in g/mol (Sanaei et al. 2019).
The results of the modelling and the contact angle calculations are presented in Fig. 2. The resulting contact angle distribution fits well to the average values identified by in situ contact angle measurements, as discussed in Sect. 3.1 Experimental data analysis.

Fluid-Fluid Interactions Modelling
The effect of fluid-fluid interactions is related to the formation of microemulsion at the interface between oil and water. Specifically, emulsification affects the flow of fluid by two mechanisms and was previously shown to be a mechanism of EOR by modified brine injections at pore scale (Masalmeh et al. 2019;Mohammadi et al. 2022).
Firstly, emulsion acts as an IFT decreasing agent. It is shown that modification of the water composition does not significantly affect the IFT between water and oil. However, the emulsions formed due to engineered water injections do have a low IFT towards both oil and water. As such, these emulsions can act as a low-IFT bridge, which allows to decrease the resisting capillary forces and improve the recovery. While it is not possible to completely reproduce the properties of in situ emulsions in laboratory environment, we have measured the IFT of various compositions of emulsions against seawater, as shown in Fig. 3. Experimental procedure is based on the works of (Mahani et al. 2015;McAuliffe 1973b) as follows: first a vial was filled with fluids of the ratio of interest, then emulsified and stabilised in the sonic bath for two hours. Then, the emulsified interface was extracted and IFT was measured using the pendant drop method (Drelich et al. 2002). Measurements were provided in the range of stable compositions. In this work, we assume that as the concentration of ion-tuned water reaches the water-oil interface, the IFT decreases. At 100% ion-tuned water, the concentration reaches the lowest measurable value (9.59 dyne/cm). Secondly, emulsions have very high viscosity compared to oil. This allows to significantly improve the microscopic displacement efficiency due to improved mobility ratio. At laboratory conditions, the viscosity of oil-water emulsion (mixed at 1:1 ratio) was measured to be 573cP using the Brookfield DVII + Cone/Plate viscometer (at 0.5 rpm to allow for the approximate representation of simulated shear rate). However, previous research has shown that the viscosity of emulsions can be as high as 2000 cP or even higher (Kokal and Al-Dokhi 2008;Umar et al. 2018). Therefore, in this study, the viscosity of the oil-ITW interface was used as a tuning parameter to compensate for history matching the microscopic recovery to that of the experimental subset.

Fluid Flow Modelling
In this section, a general mathematical description of the employed volume of fluid (VOF) method along with proposed modifications is provided. This approach was validated previously to provide highly precise description of the fluid-fluid interfaces in multiphase flow (Gamet et al. 2020;Li and Li 2018). The details of VOF methodology as well as its numerical representation and software implementation (InterFoam solver of OpenFoam CFD package) can be found in Deshpande et al. (2012).
While volume of fluid method is applied for two phase flow, it is based on a single phase form of governing CFD equations including the continuity equation and Navier-Stokes equation. The continuity equation can be expressed as follows: where biphasic phenomena are accounted throughout the averaged flow parameters of fluids' mixture such as density ( = w w + o o ) and velocity ( � ⃗ ). Subscript indices w and o denote water and oil phases, correspondingly, while is phase saturation. Navier-Stokes equation is employed in the same way as continuity equation above; however, it has additional factor (F µ ) for viscosity term and additional term responsible for capillarity: where p is the pressure and is the average dynamic viscosity of fluids' mixture ( = w w + o o ). The viscosity factor F and interfacial tension both depend on concentration of ion-tuned water ( C ) for each of the grid blocks: F = F (C) and = (C) . In turn, the equation which introduces the VOF advection scheme and provides the biphasic representation of the fluids' mixture is described with respect to the saturation of a particular phase: where relative velocity ( � ⃗ r ) is expressed as max | | � ⃗ | | � ⃗ n wo , is numerical parameter responsible for interface compression and � ⃗ n wo is an average unit normal vector of the water-oil interface.
Since two kinds of water (sea water and ion-tuned water) are considered in this research, it is convenient to introduce governing equation for the concentration of ion-tuned water just after continuity equation: where D is diffusivity coefficient. This transport equation is solved numerically using source code of conventional OpenFoam solver LaplacianFoam (Jasak et al. 2007).
To decrease the number of computational grid blocks in numerical meshes, a gradual local grid refinement is employed. A segmented micro-CT image is used for the generation of a Stereolithography (STL) image, which is then merged with a regular rectangular coarse mesh using SnappyHexMesh OpenFoam utility (Greenshields 2018). SnappyHexMesh is employed specifically to generate 3D and 2D meshes from segmented micro-CT images using hexahedra and split-hexahedra elements. The process refines and morphs the starting mesh, which should have grid blocks several times larger than the voxels of the initial micro-CT image, to closely approach the triangulated surface geometry through iterative procedures. This technique is particularly suitable for micro-CT based simulations due to its advantages, including smooth wall representation, internal connectivity preservation, reduced computational time and a lower grid block count. The SnappyHexMesh approach ensures smooth walls of porous void spaces, undistorted internal connectivity and automatic refinement of thin throats or bottlenecks in the micro-CT image, while preserving smooth wall transitions and a gradual change in the volumes of adjacent grid blocks. This combination of SnappyHexMesh and auxiliary functionality provides an efficient and suitable technique for micro-CT-based simulations.
Each simulation is carried out in three stages: primary injection of oil in a fully watersaturated sample, secondary injection of seawater and tertiary injection of the ion-tuned water. The following boundary conditions are imposed: velocity (inlet, based on the experimental value), pressure (outlet), no-flow (pore walls) and the contact angle are considered dependent on the concentration of ITW and enforced as a boundary condition with respect to saturation. OpenFoam allows considering boundary conditions as dependent on an arbitrary volume field using an interpolation table. Such interpolation table is used to relate the contact angle dependence to the ITW concentration and, thus, account for wettability alteration. Two types of simulation datasets were used in this work: an Indiana Limestone subsection extracted from micro-CT image (uses experimental parameters as input) and a 2D synthetic heterogeneous porous media designed to conceptually replicate a carbonate rock due to a variety of geometries (uses synthetic set of parameters, see Table 5).

Assumptions
A number of assumptions related to the available experimental data and state of the art of numerical modelling are incorporated in this study. These assumptions are discussed, and an overview of the future directions of research on the topic is provided:

Effects of Microporosity and Surface Roughness
Due to the limitation of micro-CT imaging technology, especially when scanning of a fluid-saturated sample is considered, a specific subset of pores is observed to be below the voxel resolution of the image. As a result, these pores are excluded from the fine scale grid meshing of pore surface and two-phase numerical fluid flow analysis. Further development studies are, therefore, essential to integrate the volumetric and connectivity effects of multiphase flow and geochemical reactivity in these pores with fluids. Furthermore, it is not possible to consider certain features of rock matrix that include spatial connectivity of subresolution pores, surface roughness, etc. as discussed earlier within the workflow of image acquisition, segmentation and fine-scale meshing.

Spatial Distribution of Wettability
It is shown in section "Experimental data analysis" that the 3D image analysis has identified a range of contact angles within sub-section of the samples, importantly within both water-wet and oil-wet regions. An average measurement is used to represent the contact angle at the oil-water-rock interface to allow us model this contact angle change in the presence of different types and concentration of the brine at a specific time step. Further development work is essential for extracting local contact angles at the specific locations of the sample and their alteration due to injection of EOR fluids in the laboratory.

Ionic Species Diffusivity
Our model considers the transport of the ionic species through both the water and oil phases. First, the effect of DLVO forces in enhanced recovery is estimated mainly through the changes in the electrical double-layer force. As per DLVO theory, a constant water film is present along the surface of the rock. The more the initial oil-wetting condition-the thinner the water film. As the modified brine is injected into the sample, ionic transport occurs and brings the change in electrostatics between the three phases, thereby changing the thickness of this film and the corresponding contact angle (Ding and Rahman 2017;Myint and Firoozabadi 2015). The other consideration comes from the numerical procedure standpoint. The usage of diffusion concerning the entire pore volume allows accounting for the change of interfacial tension and contact angle on the walls where the two phases come together. Otherwise, these phenomena are not captured appropriately. Regarding the simulation studies, the flux magnitude due to diffusion is considerably smaller than the flux magnitude due to advection. Hence, it is reasonable to assume that the diffusion takes place in all phases. It is important to note that the aim of this model is to simulate the coreflooding experiment, which is conventionally conducted to extract petrophysical data related to the EOR. In such experiment, advective transport has a higher magnitude of influence over the transport of the ionic species. As a result, all the mechanisms we propose are indeed limited to the zones near proximity of the oil-water interface at a given time during which experiments are conducted.

Kinetics of the Geochemistry
Kinetics of reactions is another important factor that was not fully considered in this study. In order to correctly characterise the timing of each mechanism in relation to the injection time, it requires further development studies. Continuum physics approach is a promising direction to provide the insights into the timescale of the modified water injection effects (Poisson-Boltzmann or Poisson-Nernst-Planck). This approach provides unique insight into ionic species diffusion and the dynamics of electrostatic interfaces (Pourakaberian et al. 2021(Pourakaberian et al. , 2022.

Viscosity Reduction Dynamics
Dynamics of viscosity reduction due to emulsification by modified salinity waterflood is another unknown, which could not be incorporated due to lack of data. It was only recently that both academia and industry have realised that fluid-fluid interactions are important and potentially key contributors to the EOR by modified waterfloods (Masalmeh et al. 2019;Tetteh et al. 2021). As such, one of the goals of this work is to provide clear pathways for more specialised research in future.

Experimental Data Analysis
A subset (300 × 300 × 300 voxels) of the micro-plug was analysed and used to validate the numerical model. Specifically, recovery factor and contact angle distributions were calculated after the secondary injection of the seawater as well as the tertiary injection of the ITW. The results of this analysis are presented in Fig. 4.
Experimental recovery data show that the recovery factor achieved by the seawater injection is 54%, while the recovery factor achieved after the tertiary injection of an ion-tuned water is 78%, with a total improvement of 24%. Wettability analysis confirms a change in wetting properties of the subsample towards more water wet, with the average contact angle across the sample after the seawater injection being 126 degrees (based on 26,800 measurements) and the average contact angle after the tertiary injection being 120 degrees (based on 24,468 measurements). The results show that the contact angle values are representative of the whole sample for which a detailed discussion is presented in previous study, see Shapoval et al. 2021.

Capillary Pressure Mechanisms Effects
In this section, we have first explored the effects of capillary-pressure affecting mechanisms, namely the contact angle reduction (attributed to rock-fluid interactions) and IFT reduction

2D Synthetic Model
To allow for an explicit visualisation of the change of the capillary effect, first we have completed a simulation study using a synthetic model of a 1000 × 1000 voxel size. This model was completely saturated with water; after that, oil was injected to achieve irreducible water saturation. Next, seawater was injected until steady-state conditions were achieved. Finally, ITW was injected in a tertiary mode. In this case, a set of synthetic parameters was used, as presented in Table 5, which clearly demonstrate the effects of various forces related to ITW application on fluid flow and oil recovery.
The results of the tertiary injection in the 2D dataset are presented in Fig. 5. It becomes apparent from the figure that as soon as the shortest available flow path is established, the enhanced recovery becomes limited to the adjacent pores. Further, the additional recovery is assured only as the ion-tuned water spreads to the surface of the rock, at which point stagnated oil is mobilised.
A similar phenomenon was previously observed by Aziz et al., (2019), who have identified stagnant zones due to the established flow path as the main limitation of further recovery increases.

3D Indiana Limestone Model
The results of the modelling of capillary mechanisms in Indiana Limestone are presented in Fig. 6. It can be observed that the oil recovery is noticeably increased from the previously unconnected zones, proximity to the main flow path of the waterflood ( Fig. 6A and B, examples of affected zones are highlighted in circles).
Quantitatively, it can be stated (see Fig. 6C) that recovery factor has increased by 6% from the secondary-(RF is 0.7) to tertiary recovery stage (RF is 0.76). However, this increase in recovery is lower than that observed in the experimental study (24%), in the same subsection of the micro-plug. This can be explained by the fact that this modelling approach accounted for only capillary effect (wettability alteration and IFT reduction) and ignored the viscous effect.

Combined Capillary and Viscous Mechanisms
The effect of the viscosity increase can be attributed to the change in the local capillary number of the porous media during the multiphase flow (see Eq. (13)). This phenomenon allows to redirect the main flow path of the fluid in the sample as the emulsion is being generated. As a result, oil recovery is increased even further by engaging the pores unaffected by the capillary pressure change. Tertiary waterflooding simulation on Indiana Limestone sample using capillary pressure affecting mechanisms. A Saturation map after secondary recovery, B saturation map after tertiary recovery, C recovery factor vs pore volumes injected (green area shows the recovery enhancement measured experimentally from the same section of the sample)

2D Synthetic Model
In this case, we have explored the change in the flow behaviour due to increase in the viscosity (by 35 times to 97.65 cP), as presented in Fig. 7. One can observe that as the ITW front advances, the previously isolated pores are connected, thereby enabling oil to flow. When compared with the capillary mechanisms, these previously isolated pores away from the main flow path are now connected to the main flow path.

3D Indiana Limestone Model
The results of the combined model are presented in Fig. 8. It is observable that the additional oil recovery becomes more evenly distributed across the volume of the sample ( Fig. 8A and B, examples of affected zones are highlighted in circles). Quantitatively (see Fig. 8C), when the capillary decrease and viscosity increase mechanism is combined, the recovery factor is improved from 0.70 to 0.92, with a total increase of 22%. This total recovery factor is reasonably close to the one observed in the experimental study (24%), which was achieved by increasing the viscosity of the ITW water-oil interface (by 35 times to 97.65cP). Notably, history matched viscosity is significantly lower than that indicative viscosity measured in the laboratory (97.65cP vs 573cP). This is related to the complexity of the formation of in situ emulsions during the fluid flow and other potential factors related to the ITW injections, which are outside of the scope of this work (time scale of wettability alteration, capillary number effect, calcite dissolution, pH change effect on IFT, etc.).

Conclusions
In this work, we have provided a set of models to represent the effect of engineered water injections and have compared the results of these models to the experimental data from a typical limestone. The results of this study can be summarised as follows: • A numerical model was built to describe the three mechanisms of enhanced oil recovery by engineered waterfloods and grouped them into capillary force affecting (contact angle and IFT reduction) and viscous forces affecting (viscosity increase). Decrease in contact angle is explained by DLVO theory while IFT reduction and viscosity increase are explained by emulsification. • Capillary force change due to ion-tuned water injection was proved to have an effect on oil recovery. However, when the capillary force is considered in the numerical model, the recovery factor achieved for the Indiana Limestone was significantly lower than that of the experiment. • When the viscous force is considered in the numerical model, it was possible to match the experimental data of the Indiana Limestone, which significantly increased sweep efficiency and connectivity of previously isolated zones of the sample. • Details of the implementation of ITW injection mechanisms can be further extended to include such components as the effect of pH change and dissolution-precipitation, which is important in the sub-porosity regions (Ansari et al. 2015;Kwak et al. 2018;