Regular Article
Experimental and CFD studies on determination of injection and production wells location considering reservoir heterogeneity and capillary number
^{1}
Enhanced Oil Recovery (EOR) Research Centre, IOR/EOR Research Institute, Shiraz University, Shiraz, Iran
^{2}
Institut für Verfahrenstechnik und Umwelttechnik, Technische Universität Dresden, Germany
^{*} Corresponding author: malayeri@shirazu.ac.ir
Received:
13
September
2018
Accepted:
18
October
2018
The indepth knowledge of reservoir heterogeneity is imperative for identifying the location of production and injection wells. The present study aimed at experimentally investigating the process of water flooding in the viscous oilsaturated glass micromodels, which contain layers with different permeability where the fractures were placed in different locations. Computational Fluid Dynamics (CFD) simulations of flooding were also conducted to study the impact of different water flow rates and wettability states. The results showed that the fractures, which have a deviation with the trend of maximum pressure gradient line, would widen the water path and vice versa. The performance of injection wells would increase the recovery factor by 18% if these would be located in the zones with high permeability for low flow rates of water. With changes in wettability state from water to oil wet conditions, the oil production will increase by 11%. Computational Fluid Dynamics results also indicated that an increase in the capillary number from 0.8 × 10^{−6} to 1.6 × 10^{−5}, would cause the recovery factor to decrease as much as 14.34% while further increase from 1.6 × 10^{−5} to 2.24 × 10^{−5}, the oil production will increase by 9.5%. Comparison between the obtained oil recoveries indicates that the maximum oil recoveries will happen when the injector well is located in the zone where ascending permeability, capillary number greater than 4.81 × 10^{−6} and also fracture with the most deviation with pressure gradient line (i.e. angular pattern) are gathered in an area between the injection and production wells.
© P. Ahmadi et al., published by IFP Energies nouvelles, 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Highlights

Significant effects of ascending and descending permeabilities between injection and production wells.

Importance of fracture location and deviation respect to maximum pressure gradient line.

Effects of capillary number on the formed water injection pattern and recovered oil.

Impact of wettability state on the trapped oil in the porous media.
Nomenclature
V_{dr,k} : Drift velocity of kth phase
V_{dr,p} : Drift velocity of a secondary phase
V_{f} : Velocity of the primary phase
V_{p} : Velocity of a secondary phase
ΔP_{1}: Pressure difference between the inlet and outlet with finer mesh
ΔP_{2}: Pressure difference between the inlet and outlet with coarser mesh
ρ_{p} : Density of a secondary phase
φ : Volume fraction of nanoparticles
φ_{p} : Volume fraction of a secondary phase
F : Force that acts on the whole body
ε : A capillary width that scales with thickness of the interface
E_{pf} : Parameter controlling interface thickness
n : Unit vector normal to the wall
б : Surface tension coefficient
γ : Reinitialization parameter
∂f/∂∅: Phiderivative of external free energy
TPF: Half of the maximum mesh element size in the model
Subscripts
1 Introduction
Energy supply still strongly depends on fossil fuels which still account for almost 80% of the required energy worldwide (Ellabban et al., 2014; Zou et al., 2016). This, in turn, increases demands for more efficient Enhanced Oil Recovery (EOR) approaches as oil reservoirs would be depleted over time. Viscous oil in the fractured reservoirs accounts approximately for 75% of the reserved oil worldwide (Alvarado and Manrique, 2010). In the past few years, there has been a growing interest in waterflooding as an EOR method. Nonetheless, drawbacks such as high oil viscosity as well as fractures and heterogeneity have led to low oil recovery of waterflooding process (Palizdan et al., 2018). During waterflooding in heterogeneous oil reservoirs especially those which contain heavy oil and/or oilwet pores, the injection of water cannot adequately push the oil towards the production well thus a significant percentage of oil would remain intact. Among numerous reasons, viscous fingering should be accounted for the reduced efficiency of waterflooding processes. Water has a higher mobility in the reservoir due to its low viscosity compared to that of oil hence by selecting the region of low resistance against the invading, it would bypass the low permeable regions (Henderson et al., 2015; Valavanides, 2012, 2018). In general, the knowledge of reservoir heterogeneity would help to appropriate approach for control of phase diversion towards low permeability channels to avoid bypass phenomena (Quennouz et al., 2014). Therefore, indepth knowledge of parameters that would impact and also control the waterflooding process is imperative to develop an efficient EOR method.
The heterogeneity of reservoir rocks is a key feature that should be given special attention when water injection is considered. Furthermore, of many decisive factors, the physical properties of injected water, geological layering and the presence or absence of fractures in the reservoirs, would substantially affect the performance of the injection process. Oil reservoirs are not usually homogeneous and that is why the capillary pressure is not the same through the porous media and consequently, the injected fluid will not be uniformly distributed. Moreover, pressure gradient which affects the oil production plays an important role in fluid movement towards lower pressure region. Wettability state of the porous medium would, in turn, strongly affect how the water/oil displacement occurs in an oil reservoir. Similarly, the reservoir heterogeneity takes a large range of spatial scales and has an imperative role in fluid flow channeling (Noetinger et al., 2004; Noetinger and Zargar, 2004). Accordingly, the characterization of parameters such as heterogeneity effect (MortonThompson and Woods, 1993), site of injection and production wells (Alvarado and Manrique, 2010; Henderson et al., 2015; Monteagudo and Firoozabadi, 2007; MortonThompson and Woods, 1993) that lead to the control of the injection process could potentially be decisive.
In the last few years in particular, the viscous fingering instability has attracted much attention since the underlying mechanisms are related to fluid flow within the porous media (Gupta et al., 1988; Hewett, 1986; Monteagudo and Firoozabadi, 2007). In both miscible and immiscible displacement processes, instabilities could happen exactly in the area where the phases are in contact (as the interface between oil and water) as a result of the interaction of different forces as well as heterogeneities in the porous medium. The forces that would influence the fingering phenomenon include viscous forces which are proportional to inverse phases’ viscosities ratio, the gravity forces due to phases’ density differences, capillary forces due to immiscible phases’ interfacial tension and dispersive forces because of miscible fluids concentration gradients (Araktingi and Orr Jr., 1993; Baveye et al., 1998; Blackwell et al., 1959; Cobb and Marek, 2003; Jha et al., 2011).
During water injection, the size of the formed flowing water path within the porous medium will be proportional to the inverse of flow rate so that for lower flow rates, the formed patterns of flowing water within the porous medium would be wider (Blackwell et al., 1959). The experimental results have shown that the fluid viscosity, interfacial tension (Homsy, 1987), and displacing phase rate (Jha et al., 2011; Pei et al., 2012; Sheorey and Muralidhar, 2003) are effective on the residual oil and the amount of oil production. The effect of viscosity should not be considered only along with the capillary number or viscosity ratio, because both will simultaneously affect the fingering pattern in the porous media. Under the oil wet conditions, for higher water injection rates, the size of the formed fingering patterns would be decreased while at the water wet condition, it would result in wider water paths when imbibition occurs. For low water injection rates, the fluids in the porous medium have much more time to be in capillary balance and as a result of that, the fingering patterns will be wider and oil production will be increased (Wijeratne and Halvorsen, 2015; Yadali Jamaloei et al., 2010). Visual observations also revealed that in highpermeable layers, the early formation and widening of multiple fingering patterns occurred, whereas fingering was postponed and more pistonlike displacement was observed in low permeable layers (Adegbite et al., 2017; Ahmadi et al., 2017; DiCarlo, 2013; Tsuji et al., 2016; Von Rosenberg, 1956; Wijeratne and Halvorsen, 2015). Investigation of the fingering phenomenon in sand stone cores, considering the variation of viscosity ratio, velocity, and diameter, shows that the developed viscous finger number does not work in an oilwet porous medium. To deal with this challenge, the correlating dimensionless number is modified for the oilwet porous medium (Worawutthichanyakul and Mohanty, 2017). At lab scale, fingering is dominant due to small heterogeneity and consequently causes oil bypassing. An effectivefingering model has been developed to scale up fingering effects. Its model parameters from history matches of a set of laboratory experiments show clear power law correlations with a dimensionless viscous finger number, a function of viscosity ratio, fluid velocity, permeability, interfacial tension, and core crosssectional area (Luo et al., 2018).
In the present study, the effects of ascending and descending permeability conditions and the impact of different water injection rates on the direction of displacing fluid movement, patterns of the formed fingers and ultimately waterflooding efficiency in glass micromodels have been experimentally investigated. Thereafter, a number of ascending permeability experimental results have been validated with simulation results of Computational Fluid Dynamics (CFD) where consistency between the simulation and experimental results are examined. The effect of wettability alteration and flow rates in the nonfractured micromodel has also been investigated using CFD.
2 Experimental setup and procedure
2.1 Micromodel and materials
2.1.1 Micromodel design
In the present study, three 60 mm × 60 mm micromodels have been designed with CorelDraw software that each one consists of three equal size regions with different permeability and octagonal shape grains (see Fig. 1). The throats size of low, intermediate and high permeable zones are 200 μm, 250 μm and 350 μm, respectively, which have been shown in different blocks as shown in Figure 1. The first micromodel is without fracture (simple pattern) and the second micromodel has angular fractures with 450 μm in width, and the fracture that is located between the high and intermediate permeable zones, creates a 45° angle with the connector line of outlet and inlet and the fracture which is located between intermediate and low permeable zones, makes an angle of 90° with the same line (angular fractures pattern). The third one has parallel fractures with 450 μm in width, that each of the fractures makes an angle of 45° with the line that connects the inlet and outlet points of the micromodel (parallel fractures pattern). Hereafter, as shown in Figure 1, in the fractured patterns, the fracture that is located between the low and intermediate permeable zones is known as fracture NO #1 and the one that is placed between the high and intermediate permeable zones as fracture NO #2. Figure 1 also indicates that the fractures are designed in such a way which the half of them is located in its upper permeable region and the other half is in the location of its lower permeable region. Flow direction has been shown with black arrows in all figures.
Fig. 1 (a) Simple pattern, (b) angular fracture pattern, (c) parallel fracture pattern, A and B are the ports that the fluid can be injected from one of them and produced from the other. 
2.1.2 Micromodel construction
To construct the glass micromodels, float glass (for more transparency) was used. By the utilization of laser technology, three porous medium patterns as shown in Figure 1 were engraved on the glass with power of 120 W, and the flow paths were carefully cleaned to achieve a perfect etching process. Thereafter, another glass with the same size was placed on the prepared etched glass and then both glasses were placed in the furnace. The furnace temperature was increased up to 750 °C with a specific procedure at which the two glasses were firmly adhered together.
2.1.3 Injecting fluid and crude oil properties
In the majority of the injection processes with the goal of Enhancing Oil Recovery (EOR), brines of various types are used. In the present study, the designed micromodels with new pattern contain regions with various throat sizes (i.e. various capillary pressure). Furthermore, in order to avoid the impact of the interaction between the polar molecules of crude oil and the brine’s ions on the oil/water interfacial tension, distilled water is used to investigate the viscous fingering phenomenon. The density and viscosity of distilled water are 0.989 gr cm^{−3} and 1.07 cP, respectively. The crude oil used was obtained from one of the southern Iranian fractured carbonate oil fields with an API (American Petroleum Institute) and viscosity values of 27.9 and 96 cP, respectively. The oil wet condition was applied to the micromodels by maintaining the oil saturated models at 60 °C for 120 h as the aging process. More details of crude oil properties are provided in Table 1.
The physical properties and SARA analysis of the crude oil.
2.1.4 Determination of wettability and interfacial tension
To determine the wettability state of the micromodels after the aging process (Ahmadi et al., 2017), contact angle measurement was conducted using cover glass, which had been aged at the abovementioned condition (Fig. 2A). The average value of Contact Angles (CA) was obtained 110° using ImageJ software. Interfacial Tension (IFT) measurement of crude oil and distilled water system was also carried out using pendant drop method, and the value of the interfacial tension of oil and water was measured to be 32 mN/m (Fig. 2B).
Fig. 2 (A) Shaped oil drop on cover glass, (B) IFT of crude oil and distilled water system. 
2.1.4.1 Capillary pressure calculation
Capillary pressure (Pc) calculation results are shown in Figure 3. Capillary pressure is defined as Pc [mN cm^{−2}] = 2σ [mN cm^{−1}] cos (θ [Degree])/r [cm], where σ is IFT, θ is contact angle and r is the pore radius. Having measured σ, θ, and range of variation in pore radius one can calculate Pc. The calculated capillary pressures with negative values of −32, −25.5, −18.3, and −14.20 (mN cm^{−2}) are related to low (200 μm), intermediate (250 μm), high (350 μm) permeable regions and fracture (450 μm), respectively. The negative sign of the capillary pressure represents the oilwet condition of the aged micromodel (θ = 110°) which withstand the invading of injected. Figure 3 indicates the capillary pressure diagram for different throat size and fracture of the parallel fracture pattern. With decreasing (more negative) of Pc, the micromodel resistance increases, since that the injected water has more tendency to enter the larger throats and fractures (with less negative capillary pressure).
Fig. 3 Capillary pressure variation in throats of the designed micromodels. 
2.1.5 Experimental procedure
The experimental procedure was as follows:

The prepared oil wet micromodels were saturated with the crude oil using a LA30 syringe pump. A schematic of experimental setup is shown in Figure 4.
Fig. 4 Schematic of micromodel setup.
Distilled water was flooded into the micromodels with three flow rates of 0.5, 1.5 and 3 mL h^{−1}, and the injection was continued until the front of the water and oil was reached the outlet of micromodel (the outlet pressure was kept constant at 1 atm during the runs and waterflooding was stopped at the moment of breakthrough (BT)).
After the occurrence of fingering, the taken photos at the breakthrough time were analyzed using an inhouse built code (based on counting of pixel numbers) to calculate oil recovery in MATLAB.
In the course of experiments, for investigating the effect of ascending permeability, distilled water was injected into port A (see Fig. 1) and the produced fluid was taken from port B, and for investigating the descending permeability vice versa.
2.1.5.1 Capillary Number (CN)
Capillary Number (CN) as a dimensionless number server as an indicator of capillary forces or viscous forces which is defined as the multiplication of viscosity and velocity values, divided by interfacial tension (Abrams, 1975; Moore and Slobod, 1955). The value of 10^{−4} is the determinant boundary to check that the viscous forces or the capillary forces will prevail in a porous medium. That means that the values which are less than it (<10^{−4}) represent the capillary dominant condition and those which are greater (>10^{−4}), indicate a viscous dominant condition.
The calculated capillary number for the flow rates of 0.5 mL h^{−1}, 1.5 mL h^{−1} and 3 mL h^{−1} are 1.6 × 10^{−6}, 4.8 × 10^{−6}, and 9.6 × 10^{−6}, respectively. The obtained values indicate that all of them are at capillary dominant condition but with increasing capillary number (1.6 × 10^{−6} < 4.8 × 10^{−6} < 9.6 × 10^{−6}) the viscous forces will be intensified and the capillary forces will be weaker. Due to the small size of the micromodels, the forces which govern the flow in porous medium should remain at capillary forces dominant range, since for higher viscous forces, a powerful pressure gradient will affect the formation of water paths.
2.2 Experimental results
2.2.1 Ascending permeability patterns tests
2.2.1.1 Simple pattern
In the nonfractured micromodel, several points were observed, which are discussed as follow. Figures 5A–5C show that immediately after the injection of distilled water, water sweeps low permeable zone and moves to the right side of the micromodel. The water path direction (Fig. 5, blue arrows) is due to increased viscous forces so that in 3 mL h^{−1} flow rate, water has washed horizontally the low permeable zone more than that of 0.5 mL hr^{1} and 1.5 mL h^{−1}. Water injection tends to move towards lower capillary pressure (Fig. 5, red arrows). As the flow rate increases further, the energy of the injected water would be greater which can overcome the generated driving force by the maximum pressure gradient. This would, in turn, converge the water paths towards the outlet. The maximum pressure gradient is also applied through the pores that are around the line with minimum length which connects the inlet and outlet of the micromodel (i.e. the diagonal of the square). Ultimately, the injected water tends to move along the paths which are affected by the pressure gradient which is basically the difference in pressures of inlet and outlet of the micromodel. As it can be seen in all flow rates, near the outlet, the pressure gradient impact has been intensified and the flow path in the high permeable zone has converged towards the outlet. In addition, in all formed patterns, the diagonal flow is quite noticeable (Fig. 5, green arrows).
Fig. 5 Breakthrough moment, simple pattern, ascending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 
For the rate of 0.5 mL h^{−1}, due to the capillary dominant condition, there is usually enough time for the capillary balance to occur at the interface between water and oil and grains, and consequently the formed pattern is fatter (wider) than those that are formed at the condition which the viscous forces are dominant (i.e. 1.5 mL h^{−1} and 3 mL h^{−1}).
2.2.1.2 Parallel fractures pattern
In the micromodel with parallel fractures pattern, the fractures #1, and 2 are located among the grains that are along the line of maximum pressure gradient. Thus, due to the low capillary pressure of the fractures, the convergence of the water paths would occur towards the micromodel outlet (Figs. 6A–6C, blue arrows). Red arrows in Figure 6A, show the water path would be getting wider similar to that shown in Figure 5A of the simple pattern. Any increase in flow rate would intensify the energy of injected water hence for higher flow rates, viscous forces would be able to overcome the capillary pressure. This would, in turn, cause water to invade in the different directions with the maximum pressure gradient (Figs. 6A–6C, green arrows).
Fig. 6 BT moment, parallel fractures pattern, ascending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 
2.2.1.3 Angular fractures pattern
In these runs, horizontal water invading into low permeable zone would increase at such flow rates of 1.5 mL h^{−1} and 3 mL h^{−1} showed more zones that are washed by the injected water (Figs. 7A–7C, red arrows). Fracture #1 has an angle of 90° along the maximum pressure gradient line. This would instead causes water to flow through the grains for getting the outlet of the micromodel (Figs. 7B and 7C, blue arrows). Accordingly, for the flow rate of 3 mL h^{−1}, fracture #2 has not been filled completely with water (Fig. 7C). At a lower flow rate of 0.5 mL h^{−1}, due to the more laminarity of the water flowing regime, its path is influenced by the pressure gradient and consequently after passing through the grains that are along the maximum gradient line, the injected water would then reach the outlet.
Fig. 7 BT moment, angular fractures pattern, ascending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 
2.2.2 Experiments on descending permeability patterns
2.2.2.1 Simple pattern
As shown in Figure 8, because of minimum capillary pressure in the high permeable region, the injected water would invade the micromodel much preferably than medium and low permeable zones. Due to the capillary dominant condition, during the injection of water with the flow rate of 0.5 mL h^{−1}, the capillary forces would prevail as it is shown in Figure 8A, the higher permeable region was washed more than the same regions of Figures 8B and 8C (red arrows). This can also be shown in Figures 8B–8C (purple ovals) as there are more trapped oil droplets for a flow rate of 3 mL h^{−1} flow in the same high permeable regions. During the invasion of the injected water in the high permeable zone (Figs. 8A–8C, green arrow), the impact of pressure drop impact will be more noticeable as the injected water enters the intermediate permeable zone from different positions (at the boundary of high and intermediate permeable regions) and will deviate its path to reach the outlet of the micromodels due to increasing the pressure gradient so that the water bypasses the oil in some areas (Figs. 8A–8C, red square dot). The low permeable regions (above the brown lines) in the oil wet micromodels have no tendency to allow the water to enter. This lies in the fact that in the highpressure gradient near the outlet, the viscous forces overcome the capillary pressure causing water to reach the production port.
Fig. 8 BT moment, simple pattern, descending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 
2.2.2.2 Parallel fractures pattern
During water injection, the fractures with the low capillary pressure regions in an oilwet porous medium, the injected water tends to move through the fractures in order to reach the outlet. The results also indicate that the existence of the fractures in the same direction of the pressure gradient lead the shrinkage of the formed pattern by water in comparison with the simple pattern results. The blue arrows in Figures 9A–9C outline the direction of maximum pressure drop that the water selects to pass through the fractures and grains to reach the micromodel’s outlet.
Fig. 9 BT moment, parallel fractures pattern, descending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 
2.2.2.3 Angular fractures pattern
The direction of water flow in fracture #1 is perpendicular to the direction of that which is flowed towards the outlet (Figs. 10A–10C, blue arrows). In this case, the formed patterns demonstrate that the near the outlet fracture due to the 90° angle which makes with the maximum pressure gradient direction, would weaken the impact of pressure drop. Therefore, even with the presence of a fracture (#2) near the inlet, the injected water still invades through the high permeable region (Figs. 10A–10C, green arrows). The orientation of the fracture plays a significant role in determining the amount of injected water that would move the fractures or passes through the grains. If one compares all the results obtained in this study, the angular fractures pattern represents the maximum oil recovery for similar flow rates and fracture location (Fig. 10).
Fig. 10 BT moment, angular fractures pattern, descending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 
3 CFD simulation
To validate the accuracy of analyses described in the preceding section which were based on the underlying importance of the capillary or viscous forces, a simulation was also conducted with the help of CFD (Ghaffarkhah et al., 2017; Jafari et al., 2008; Mohammadmoradi and Kantzas, 2016a, 2016b; Moraveji et al., 2017; Nadimi et al., 2016; Xia et al., 2017). To do so, three micromodels with three different permeability layers, as shown in Figure 1, have been considered to investigate the impact of micromodel heterogeneity on the paths that the fluid chooses to flow. Mass transport and momentum equations were solved based on the finite volume method using CAMSOL 5.1 software to investigate the water and oil displacement in the porous medium. The governing equations that can model n phases (liquids and solid) through solving the continuity, momentum for the mixtures model are as follows:
Conservation of momentum:$$\frac{\partial}{\partial t}\left({\rho}_{\mathrm{m}}\right)+\nabla \bullet \left({\rho}_{\mathrm{m}}\stackrel{\u20d7}{{V}_{\mathrm{m}}}\right)=0.$$(1)
Conservation of mass:$$\frac{\mathit{\partial}}{\mathit{\partial}\mathit{t}}\left({\mathit{\rho}}_{\mathbf{m}}\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{m}}}\right)+\nabla \bullet \left({\mathit{\rho}}_{\mathbf{m}}\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{m}}}\stackrel{\u20d7}{{\mathit{v}}_{\mathbf{m}}}\right)=\nabla {\mathit{P}}_{\mathbf{m}}+\nabla \bullet \left({\mathit{\mu}}_{\mathbf{m}}\nabla \stackrel{\u20d7}{{\mathit{V}}_{\mathbf{m}}}\right)+\nabla \bullet \left(\sum _{\mathit{k}=\mathbf{1}}^{\mathit{n}}{\mathit{\varnothing}}_{\mathit{k}}{\mathit{\rho}}_{\mathit{k}}\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{dr},\mathit{k}}}\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{dr},\mathit{k}}}\right).$$(2)
kth phase drift velocity is:$$\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{dr},\mathit{k}}}=\stackrel{\u20d7}{{\mathit{V}}_{\mathit{k}}}\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{m}}}.$$(3)
The second phase volume fraction equation can be obtained by:$$\frac{\mathit{\partial}}{\mathit{\partial}\mathit{t}}\left({\mathit{\varnothing}}_{\mathbf{p}}{\mathit{\rho}}_{\mathbf{p}}\right)+\nabla \bullet \left({\mathit{\varnothing}}_{\mathbf{p}}{\mathit{\rho}}_{\mathbf{p}}\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{m}}}\right)=\nabla \bullet \left({\mathit{\varnothing}}_{\mathbf{p}}{\mathit{\rho}}_{\mathbf{p}}\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{dr},\mathbf{p}}}\right).$$(4)
The relative velocity (the velocity of the secondary phase relative to the velocity of the primary phase) and also the drift velocity related to the relative velocity are as follow:$$\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{pf}}}=\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{p}}}\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{f}}},$$(5)$$\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{dr},\mathbf{p}}}=\stackrel{\u20d7}{{\mathit{V}}_{\mathbf{pf}}}\sum _{\mathit{k}=\mathbf{1}}^{\mathit{n}}\frac{{\mathit{\varnothing}}_{\mathit{k}}{\mathit{\rho}}_{\mathit{k}}}{{\mathit{\rho}}_{\mathbf{eff}}}\stackrel{\u20d7}{{\mathit{V}}_{\mathit{fk}}}.$$(6)
For the CFD simulation, distilled water and crude oil were considered as displacing and displaced fluids, respectively. Their physical properties are reported in the previous section. Waterflooding was simulated for the same flow rates of ascending permeability to validate the accuracy of the experimental results. Thereafter, the effect of flow rate and wettability state for different conditions have been examined. In all figures of CFD simulations, the oil and water phases are shown in red and blue colours, respectively. The discretization approaches and boundary conditions for the present simulation were as follows:
Mesh: As expected, the finest mesh size would lead to more accuracy but at the expense of longer computing time. In the present CFD simulation, to boost the computing accuracy, from building mesh module, the physicscontrolled mesh has been selected from sequence type and extremely fine was chosen as mesh size.
Flow: The equations of laminar twophase flow were used.
Interfacetracking: In this study, the phase field flow has been used where the interface of two fluids was calculated indirectly by means of a phase field variable. The NavierStokes equations are manipulated by adding the surface tension force, as a body force by multiplying the chemical potential of the system by the gradient of the phase field variable as per following equations (Provatas and Elder, 2011; Wörner, 2012):$$\rho \frac{\partial u}{\partial t}=\nabla \xb7[\mathrm{pl}+\mu \left(\nabla u+{\left(\nabla u\right)}^{T}\right]+\rho g+{F}_{\mathrm{st}}+F,$$(7)$$\nabla \xb7u=0,$$(8)$$\frac{\partial \phi}{\partial t}+u\xb7\nabla \varnothing =\nabla \xb7\frac{\gamma \lambda}{{\u03f5}_{\mathrm{pf}}^{2}}\nabla \Psi ,\varnothing =\varnothing \mathrm{pf}$$(9)$$\Psi =\nabla \xb7{\u03f5}_{\mathrm{pf}}^{2}\nabla \varnothing +\left({\varnothing}^{2}1\right)\varnothing +\frac{{\u03f5}_{\mathrm{pf}}^{2}}{\lambda}{\partial}_{f}\partial \varnothing ,\Psi =\mathrm{psi}$$(10)$$\lambda =\frac{3{\u03f5}_{\mathrm{pf}}^{2}\delta}{\sqrt{8}},\gamma =\chi {\u03f5}_{\mathrm{pf}}^{2}$$(11)
Boundary conditions:
For injection port the Inlet boundary condition was considered constant rate during each test.
For production port the outlet boundary condition was considered constant atmospheric pressure.
The wettedwall node boundary condition was considered with contact angle of 110° along with zero fluid velocity (U = 0).
The phase field would satisfy the following boundary condition:$$n\bullet \frac{\gamma \lambda}{{\u03f5}_{\mathrm{pf}}^{2}}\nabla \Psi =0.$$(12)
The flowing boundary condition defines a contact angle between wall and the fluid:$$n\bullet {\u03f5}_{\mathrm{pf}}^{2}\nabla \varnothing ={\u03f5}_{\mathrm{pf}}^{2}\mathrm{cos}\left({\theta}_{w}\right)\left\nabla \varnothing \right.$$(13)
5Initial value: For conducting each CFD simulation run, the required values of pressure (=1 atm), velocity (=0 at x and ydirections), and initial saturation (=100 for oil) are selected.
3.1. CFD validation of experimental runs of ascending permeability patterns
In this study, to simulate the injected water in the oil wet porous medium, the fluids are considered incompressible. The oil and water densities were defined to be 0.948 gr cc^{−1} and 0.988 gr cc^{−1} and viscosities to be 96 cP and 1 cP, respectively. The contact angle of slightly oil wet glass with oil droplet was measured to be approximately 110° for all grid points of the patterns using contact angle measurement method (Fig. 2). A comparison of the simulation results (Fig. 11) with the experimental results (Figs. 5–7) confirms, to some extent, that they are qualitatively comparable. Following, the impact of wettability state and different flow rates are examined and discussed.
Fig. 11 BT moment (A) simple pattern, (B) parallel fracture pattern, (C) angular fracture pattern: (1) 0.5 mL h^{−1}, (2) 1.5 mL h^{−1}, (3) 3 mL h^{−1}. 
The comparison of the simulation results presented in Figures 11A1–11A3 and experimental results at similar conditions (see Fig. 5) confirmed that for the flow rate of 0.5 mL h^{−1}, the expansion of injected water pattern is evident. This lies in dominance of capillary forces compared to those for higher flow rates. For flow rates of 1.5 mL h^{−1} and 3 mL h^{−1} intensifying the viscose forces has caused the deviation of water/oil front towards the high permeable region. The presence of water/oil front would lead to invasion of low permeable porous medium by the injected water as shown with green arrows in Figures 6B and 6C (experimental results) and Figures 11B2, and 11B3 (CFD results). As Figures 7B and 7C (experimental results) and Figures 11C2 and 11C3 demonstrate, failure to completely enter the flowing water to the fracture #2 (in the high permeable region) are related to the impact of pressure gradients boost near the outlet. This is because when the flow rate increases then in the case of 3 mL h^{−1} the tendency of flowing water to enter the fracture would significantly be decreased and in the CFD results (Fig. 11C3) the near outlet fracture has not been filled with water.
3.2 Impact of flow rate on water flow paths in the simple pattern
For investigation of viscous fingering, the flow rates of 0.25 mL h^{−1} and 0.5 mL h^{−1} (CN = 0.8 × 10^{−6} and 1.60 × 10^{−6}, respectively) in the capillary dominant range and 5 mL h^{−1} and 7 mL h^{−1} (CN = 1.6 × 10^{−5} and 2.24 × 10^{−5}, respectively) in the close range to viscous dominant condition, all with a contact angle of 110 degrees, are considered.
As it can be seen from Figures 12A to 12D, in the low permeable zone, with higher flow rates, the invaded zone has been increased as viscous forces would be becoming more dominant. Under such circumstances, the capillary forces and pressure gradient are not able to resist, thus a more washed area is expected. When the rate of injected water increases, the viscous forces will overcome the pressure gradient thus the injected water will deviate its way and instead invades the low permeable zone (Fig. 12, black arrows).
Fig. 12 BT moment, slightly oil wet state (110°), invaded zones in the low permeability layer due to viscous forces: (A) flow rate = 0.25 mL h^{−1}, (B) flow rate = 0.5 mL h^{−1}, (C) flow rate = 5 mL h^{−1}, (D) flow rate = 7 mL h^{−1}. 
The result for the flow rate of 0.25 mL h^{−1} (Fig. 12A) also proves that for very low viscous forces, there is a very lowpressure gradient. In addition, there is enough time to achieve a capillary balance between oil and water and grains in the porous medium. As a result, in the region near the inlet the washed zone is wider than that of the higher flow rates in the same regions. When the water is getting close to the outlet, it somehow feels the pressure gradient more than before and will converge towards the outlet. As it is shown in Figure 12, in the oil wet micromodel, the high permeable zones overcome the low capillary pressure region causing the injected water to move towards the outlet (Fig. 12, red arrows). This means that the low capillary pressure overcomes the pressure gradient and therefore the water chooses the high permeable region instead of passing through the low and intermediate permeable throats to reach the outlet.
3.3 Impact of wettability state on water paths in the simple pattern
In this part, the focus was to investigate the effect of wettability state on the condition that no viscous forces could be able to seriously influence the wettability role. Therefore, the CA of 30°, 110°, 125° and 140° were chosen to run CFD with flow rate of 1 mL h^{−1}.
As it can be seen from Figure 13, in the water wet condition (CA = 30°) the tendency of water to remain more pistonlike was more than that of the other wettability states. Furthermore, when the oil/water front has reached near the outlet, the pressure gradient has been amplified so that the water reaches the outlet directly along the line of maximum gradient pressure.
Fig. 13 BT moment, different wettability states in the simple ascending permeability pattern: (A) CA = 30°, (B) CA = 110°, (C) CA = 125°, (D) CA = 140°. 
Adhesion of water to the grains causes no trapped oil between the throats (Pei et al., 2012) and pores but with increasing the contact angle towards oilwetness, the amount of trapped oil has been increased so that in the contact angle of 140° has reached the highest trapping compared to the other cases. In fact, as the oil wetness of the pores and throats increases, the pores and throats would have less tendency to drain oil thus the amount of trapped oil would be maximum for a contact angle of 140°.
4 Oil recoveries
4.1 Experimental oil recovery of descending and ascending permeability patterns
Figure 14 shows that for the low flow rate of 0.5 mL h^{−1}, due to its slow movement and also low permeable region near the outlet, the tendency of water flow towards it, will be reduced and consequently water will enter into the high permeable regions that are around the inlet. Thus, for higher flow rates of 1.5 mL h^{−1} and 3 mL h^{−1}, the injected water would be able to pass through the low permeable throats that are located around the outlet.
Fig. 14 Experimental recovery factor of descending permeability patterns oilwet state, flow rates of 0.5 mL h^{−1}, 1.5 mL h^{−1}, and 3 mL h^{−1}. 
The calculated RF of the same flow rates of 3 mL h^{−1} in the angular pattern has shown a superiority of about 18% in the case of injection into an ascending permeability region. Comparison between the obtained experimental oil recoveries (Figs. 14 and 15) indicates that the maximum oil recoveries will happen when the injector well is located in the zone where ascending permeability, CN greater than 4.81 × 10^{−6} (1.5 mL h^{−1}) and also fracture with the most deviation with pressure gradient line (i.e. angular pattern) are gathered in an area between the injection and production wells.
Fig. 15 Experimental recovery factor of ascending permeability patterns, oilwet state, flow rates of 0.5 mL h^{−1}, 1.5 mL h^{−1}, and 3 mL h^{−1}. 
4.2 Oil recovery comparison of experimental and CFD simulation for ascending permeability patterns
As Figure 11 shows, the findings of CFD runs of ascending permeability patterns for the flow rates of 0.5 mL h^{−1}, 1.5 mL h^{−1}, and 3 mL h^{−1} imply that there are a good accuracy with those of experimental results. Figure 17 illustrates a better comparison of the normalized RF of experimental runs and their respective CFD simulations. Furthermore, the similarity of oil recovery trends as shown in Figures 15 and 16 confirms that for different flow rates and patterns for both experimental and CFD results, the oil recoveries are comparable.
Fig. 16 CFD runs recovery factor of ascending permeability patterns, oilwet state, 0.5 mL h^{−1}, 1.5 mL h^{−1}, and 3 mL h^{−1}. 
Fig. 17 Comparison of experimental and CFD results at BT moment for pattern with ascending permeability. 
4.3 Oil recovery estimation for different flow rates using CFD
The oil recovery (Fig. 18) for the runs of Section 3.2 demonstrate a decreasing and increasing trends due to the roles of capillary and viscous forces. As stated earlier, the dominancy of each of these forces will determine the flowing water paths and consequently the amount of oil recovery.
Fig. 18 CFD simulation recovery factor of descending permeability (simple pattern), flow rate investigation, oilwet state. 
Comparison of Figures 15 and 18 shows that RF reduction happens in flow rate of 1.5 mL h^{−1} while the CFD simulation of the same porous medium pattern demonstrates this RF reduction in the flow rate of 5 mL h^{−1}. This difference is likely due to the governing equation which controls the oil/water front velocity in CFD simulation so that in the 5 mL h^{−1}, the velocity (V_{dr,p} – Eq. (6)) of the segment of the front which its direction is towards the outlet is greater than the other front segments and consequently the washed area in Figure 12C (which is shown with black arrow) is less than that of Figures 12B and 12C. Although in experimental results the role of viscous force on the amount of oil recovery is observed at lower flow rate (1.5 mL h^{−1}) which is probably due to changes in uncontrollable parameters such as wettability and its impacts on the phases velocities.
4.4 Oil recovery of the investigated different wettability state using CFD simulation
As it can be seen from Figure 19, the change of wettability towards oil wet condition, would lead to the reduction of oil recovery due to the increased adhesion of oil to the grains and reduction of oil permeability.
Fig. 19 CFD simulation recovery factor of descending permeability (simple pattern), contact angle investigation. 
5 Conclusion
The experimental and simulation results obtained in this study, showed that:
Generally speaking, the lower injection rate in the patterns that the high permeable regions are close to the injection well will cause more oil recovery than those that are close to the outlet. This is likely because the high permeable region and increasing the pressure gradient will facilitate the reduction of breakthrough time.
The presence of the fractures which are at the same direction of maximum pressure gradient, would result in lower oil recovery while the fractures with the proper angle would help the distribution of the injected water in the porous medium.
For the high flow rates, due to ability of the injected water to invade through the low permeable zones, the amount of oil recovery will increase.
For low flow rates, the comparison of descending and ascending permeability shows that the superiority of descending permeability condition.
In water wet condition, the flow paths would move more pistonlike in the porous medium and causes more oil production. On the contrary, if for more oil wetness micromodels, the amount of trapped oil will increase.
Acknowledgments
The authors would like to acknowledge the support and assistance of EOR Research Center of Shiraz University during this study.
References
 Abrams A. (1975) The influence of fluid viscosity, interfacial tension, and flow velocity on residual oil saturation left by waterflood, Soc. Pet. Eng. J. 15, 437–447. [CrossRef] [Google Scholar]
 Adegbite J.O., AlShalabi E.W., Ghosh B. (2017) Modeling the effect of engineered water injection on oil recovery from carbonate cores, SPE International Conference on Oilfield Chemistry, Society of Petroleum Engineers. [Google Scholar]
 Ahmadi P., Riazi M., Malayeri M.R. (2017) Investigation of co and counter current flow behaviour in carbonate rock cores, International Symposium of the Society of Core Analysts, held in Vienna, Austria, 27–30 August. [Google Scholar]
 Alvarado V., Manrique E. (2010) Enhanced oil recovery: an update review, Energies 3, 1529–1575. [CrossRef] [Google Scholar]
 Araktingi U.G., Orr F.M. Jr. (1993) Viscous fingering in heterogeneous porous media, SPE Adv. Technol. Ser. 1, 71–80. [CrossRef] [Google Scholar]
 Baveye P., Vandevivere P., Hoyle B.L., DeLeo P.C., de Lozada D.S. (1998) Environmental impact and mechanisms of the biological clogging of saturated soils and aquifer materials, Crit. Rev. Environ. Sci. Technol. 28, 123–191. [Google Scholar]
 Blackwell R.J., Rayne J.R., Terry W.M. (1959) Factors influencing the efficiency of miscible displacement, Society of Petroleum Engineers, SPE1131G. [Google Scholar]
 Cobb W.M., Marek F.J. (2003) Determination of volumetric sweep efficiency in mature waterfloods using production data, SPE Repr. Ser. 182–190. [Google Scholar]
 DiCarlo D.A. (2013) Stability of gravitydriven multiphase flow in porous media: 40 Years of advancements, Water Resour. Res. 49, 4531–4544. [Google Scholar]
 Ellabban O., AbuRub H., Blaabjerg F. (2014) Renewable energy resources: Current status, future prospects and their enabling technology, Renew. Sustain. Energy Rev. 39, 748–764. [Google Scholar]
 Ghaffarkhah A., Shahrabi M.A., Moraveji M.K., Eslami H. (2017) Application of CFD for designing conventional three phase oilfield separator, Egypt. J. Pet. 26, 413–420. [CrossRef] [Google Scholar]
 Gupta A.D., Pope G.A., Sepehrnoori K., Shook M. (1988) Effects of reservoir heterogeneity on chemically enhanced oil recovery, SPE Reserv. Eng. 3, 479–488. [CrossRef] [Google Scholar]
 Henderson N., Brêttas J.C., Sacco W.F. (2015) Applicability of the threeparameter KozenyCarman generalized equation to the description of viscous fingering in simulations of waterflood in heterogeneous porous media, Adv. Eng. Softw. 85, 73–80. [Google Scholar]
 Hewett T.A. (1986) Fractal distributions of reservoir heterogeneity and their influence on fluid transport, Society of Petroleum Engineers. [Google Scholar]
 Homsy G.M. (1987) Viscous fingering in porous media, Annu. Rev. Fluid Mech. 19, 271–311. [Google Scholar]
 Jafari A., Zamankhan P., Mousavi S.M., Pietarinen K. (2008) Modeling and CFD simulation of flow behavior and dispersivity through randomly packed bed reactors, Chem. Eng. J. 144, 476–482. [Google Scholar]
 Jha B., CuetoFelgueroso L., Juanes R. (2011) Fluid mixing from viscous fingering, Phys. Rev. Lett. 106, 194502. [CrossRef] [PubMed] [Google Scholar]
 Luo H., Delshad M., Pope G.A., Mohanty K.K. (2018) Scaling up the interplay of fingering and channeling for unstable water/polymer floods in viscousoil reservoirs, J. Pet. Sci. Eng. 165, 332–346. [Google Scholar]
 Mohammadmoradi P., Kantzas A. (2016a) Petrophysical characterization of porous media starting from microtomographic images, Adv. Water Resour. 94, 200–216. [Google Scholar]
 Mohammadmoradi P., Kantzas A. (2016b) Porescale permeability calculation using CFD and DSMC techniques, J. Pet. Sci. Eng. 146, 515–525. [Google Scholar]
 Monteagudo J.E.P., Firoozabadi A. (2007) Controlvolume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects, SPE J. 12, 355–366. [CrossRef] [Google Scholar]
 Moore T.F., Slobod R.L. (1955) Displacement of oil by watereffect of wettability, rate, and viscosity on recovery, Fall Meeting of the Petroleum Branch of AIME, Society of Petroleum Engineers. [Google Scholar]
 Moraveji M.K., Sabah M., Shahryari A., Ghaffarkhah A. (2017) Investigation of drill pipe rotation effect on cutting transport with aerated mud using CFD approach, Adv. Powder Technol. 28, 1141–1153. [Google Scholar]
 MortonThompson D., Woods A.M. (1993) Development geology reference manual: AAPG methods in exploration series, no. 10, AAPG. [Google Scholar]
 Nadimi S., Miscovic I., McLennan J. (2016) A 3D peridynamic simulation of hydraulic fracture process in a heterogeneous medium, J. Pet. Sci. Eng. 145, 444–452. [Google Scholar]
 Noetinger B., Artus V., Ricard L. (2004) Dynamics of the water–oil front for twophase, immiscible flow in heterogeneous porous media. 2isotropic media, Transp. Porous Media 56, 305–328. [Google Scholar]
 Noetinger B., Zargar G. (2004) Multiscale description and upscaling of fluid flow in subsurface reservoirs, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 59, 119–139. [CrossRef] [Google Scholar]
 Palizdan S., Ahmadi P., Riazi M., Malayeri M.R. (2018) Effect of formed waterinoil emulsions by cobalt oxide nanoparticles on the oil recovery, an experimental approach, 80th EAGE Conference and Exhibition 2018, Copenhagen, Denmark. [Google Scholar]
 Pei H., Zhang G., Ge J., Tang M., Zheng Y. (2012) Comparative effectiveness of alkaline flooding and alkalinesurfactant flooding for improved heavyoil recovery, Energy Fuel. 26, 2911–2919. [CrossRef] [Google Scholar]
 Provatas N., Elder K. (2011) Phasefield methods in materials science and engineering, John Wiley & Sons. [Google Scholar]
 Quennouz N., Ryba M., Argillier J.F., Herzhaft B., Peysson Y., Pannacci N. (2014) Microfluidic study of foams flow for enhanced oil recovery (EOR), Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 69, 457–466. [CrossRef] [Google Scholar]
 Sheorey T., Muralidhar K. (2003) Isothermal and nonisothermal oilwater flow and viscous fingering in a porous medium, Int. J. Therm. Sci. 42, 665–676. [Google Scholar]
 Tsuji T., Jiang F., Christensen K.T. (2016) Characterization of immiscible fluid displacement processes with various capillary numbers and viscosity ratios in 3D natural sandstone, Adv. Water Resour. 95, 3–15. [Google Scholar]
 Valavanides M.S. (2012) Steadystate twophase flow in porous media: Review of progress in the development of the DeProF theory bridging pore to statistical thermodynamics scales, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 67, 787–804. [CrossRef] [Google Scholar]
 Valavanides M.S. (2018) Oil fragmentation, interfacial surface transport and flow structure maps for twophase flow in model pore networks. Predictions based on extensive, DeProF model simulations, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 73, 6. [Google Scholar]
 Von Rosenberg D.U. (1956) Mechanics of steady state singlephase fluid displacement from porous media, AIChE J. 2, 55–58. [Google Scholar]
 Wijeratne D.I.E.N., Halvorsen B.M. (2015) Computational study of fingering phenomenon in heavy oil reservoir with water drive, Fuel 158, 306–314. [CrossRef] [Google Scholar]
 Worawutthichanyakul T., Mohanty K.K. (2017) Unstable immiscible displacements in oilwet rocks, Transp. Porous Media 119, 205–223. [Google Scholar]
 Wörner M. (2012) Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications, Microfluid. Nanofluidics 12, 841–886. [Google Scholar]
 Xia Y., Goral J., Huang H., Miskovic I., Meakin P., Deo M. (2017) Manybody dissipative particle dynamics modeling of fluid flow in finegrained nanoporous shales, Phys. Fluids 29, 56601. [CrossRef] [Google Scholar]
 Yadali Jamaloei B., Kharrat R., Torabi F. (2010) Analysis and correlations of viscous fingering in lowtension polymer flooding in heavy oil reservoirs, Energy Fuel. 24, 6384–6392. [CrossRef] [Google Scholar]
 Zou C., Zhao Q., Zhang G., Xiong B. (2016) Energy revolution: From a fossil energy era to a new energy era, Nat. Gas Ind. B 3, 1–11. [Google Scholar]
All Tables
All Figures
Fig. 1 (a) Simple pattern, (b) angular fracture pattern, (c) parallel fracture pattern, A and B are the ports that the fluid can be injected from one of them and produced from the other. 

In the text 
Fig. 2 (A) Shaped oil drop on cover glass, (B) IFT of crude oil and distilled water system. 

In the text 
Fig. 3 Capillary pressure variation in throats of the designed micromodels. 

In the text 
Fig. 4 Schematic of micromodel setup. 

In the text 
Fig. 5 Breakthrough moment, simple pattern, ascending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 

In the text 
Fig. 6 BT moment, parallel fractures pattern, ascending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 

In the text 
Fig. 7 BT moment, angular fractures pattern, ascending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 

In the text 
Fig. 8 BT moment, simple pattern, descending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 

In the text 
Fig. 9 BT moment, parallel fractures pattern, descending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 

In the text 
Fig. 10 BT moment, angular fractures pattern, descending permeability: (A) 0.5 mL h^{−1}, (B) 1.5 mL h^{−1}, (C) 3 mL h^{−1}. 

In the text 
Fig. 11 BT moment (A) simple pattern, (B) parallel fracture pattern, (C) angular fracture pattern: (1) 0.5 mL h^{−1}, (2) 1.5 mL h^{−1}, (3) 3 mL h^{−1}. 

In the text 
Fig. 12 BT moment, slightly oil wet state (110°), invaded zones in the low permeability layer due to viscous forces: (A) flow rate = 0.25 mL h^{−1}, (B) flow rate = 0.5 mL h^{−1}, (C) flow rate = 5 mL h^{−1}, (D) flow rate = 7 mL h^{−1}. 

In the text 
Fig. 13 BT moment, different wettability states in the simple ascending permeability pattern: (A) CA = 30°, (B) CA = 110°, (C) CA = 125°, (D) CA = 140°. 

In the text 
Fig. 14 Experimental recovery factor of descending permeability patterns oilwet state, flow rates of 0.5 mL h^{−1}, 1.5 mL h^{−1}, and 3 mL h^{−1}. 

In the text 
Fig. 15 Experimental recovery factor of ascending permeability patterns, oilwet state, flow rates of 0.5 mL h^{−1}, 1.5 mL h^{−1}, and 3 mL h^{−1}. 

In the text 
Fig. 16 CFD runs recovery factor of ascending permeability patterns, oilwet state, 0.5 mL h^{−1}, 1.5 mL h^{−1}, and 3 mL h^{−1}. 

In the text 
Fig. 17 Comparison of experimental and CFD results at BT moment for pattern with ascending permeability. 

In the text 
Fig. 18 CFD simulation recovery factor of descending permeability (simple pattern), flow rate investigation, oilwet state. 

In the text 
Fig. 19 CFD simulation recovery factor of descending permeability (simple pattern), contact angle investigation. 

In the text 