Two-component fluid front tracking in fault zone and discontinuity with permeability heterogeneity

Fluid front tracking is important in two-phase/component fluid flow in porous media with different heterogeneities, especially in the improved recovery of oil. Three different flow patterns of stable, viscous fingering, and capillary finger ing exist based on the fluids’ viscosity and capillary number (CA). In addition, fluid front and sweep efficiency are af fected by the heterogeneity of the porous medium. In the current study, the heterogeneous porous media are: (1) normal fault zone or cross-bedding with heterogeneity in permeability, and (2) a fracture or discontinuity between two porous media consisting of two homogeneous layers with very low and high permeabilities, in which immiscible water flooding is performed for sweep efficiency and streamlines tracking purposes. By considering the experimental glass micromodel and the simulation results of discontinuity, a crack is the main fluid flow path. After the breakthrough, fluid inclines to penetrate the fine and coarse grains around the crack. Moreover, an increase in flow rate from 1 and 200 (ml/h) in both the experimental and simulation models causes a reduction in the sweep efficiency from 14% to 7.3% and 15.6% to 10% by the moment of breakthrough, respectively. In the fault zone, the sweep efficiency and the streamline of the injected fluid showed a dependency on the interface incident angle, and the layers’ permeability. The presented glass micro model and Lattice Boltzmann Method were consistent with fluid dynamics, and both of them were suitable for a precise evaluation of sweep efficiency and visualization of preferential pathway of fluid flow through cross-bedding and discon tinuity for enhanced oil recovery purposes.


Introduction
Heterogeneity is important in all scales from several kilometres to centimetres, especially for the determination of primary oil saturation and oil distribution of sweep zone by an improved oil recovery process. An increase in heterogeneity and anisotropy depends on different conditions from transportation and sedimentation to diagenesis, lithification, tectonic activity, and dissolution (Bjoerlykke, 1989). Chasing fluid front in definite geometries and heterogeneities is important in enhanced oil recovery. In small scales, because of the lack of visualizations, the interpretation is challenging and vital for the overall evaluation of the main oil recovery mechanism. Fluid flow in fault/cross-bedding structure with permeability contrasts and also, in a crack between two layers with different permeabilities attracted the attention of several researchers (Weber, 1982;Kortekaas, 1985;Ringrose et al., 1993;Huang et al., 1995;Silva et al. had a sequential study on wettability and permeability heterogeneity in some structures in pore-scale as follows: in 1990, the effects of permeability contrast, lens size, viscosity ratio, and lens heterogeneity on fluid flow were numerically considered and the results were experimentally compared with a glass-bead pack (McKean and Dawe, 1990). In 1992 and 1993, miscible and immiscible fluid displacement was experimentally performed on a lens and cross-bedding structure by a glassbead pack. In this study, layer thickness, cross-bedding angle, permeability heterogeneity, mobility ratio, and flow rate were analysed. The results showed that fluid front only depends on the homogeneity of the porous medium, capillary pressure, and fluid saturation (Dawe et al., 1992;Roti and Dawe, 1993). In 1996, the effects of wettability and permeability heterogeneity on the fluid front in cross-bedding, which depicts the sedimentary structure in the fluvial deposits or the fault zones, were investigated (Caruana and Dawe, 1996). In 2010 and 2011, the vertical layers' effects in fluid flow direction with permeability contrast were studied in pore-scale by an etched glass model (Dawe et al., 2010(Dawe et al., , 2011. Chaudhari et al., in immiscible oil/water displacement in three parallel layers with different permeabilities, showed that during the drainage process, the oil prefers to enter the pores with greater permeability. Also, during the imbibition process, the water prefers to enter the pores with a lower permeability because of the wettability and capillary pressure effects (Chaudhari et al., 2011).
Due to heterogeneity in all scales, chasing fluid front in definite geometries and heterogeneities is important. In large scales, clear interpretations show that heterogeneity affects the displacement process, but in small scales, because of the lack of visualizations, the interpretation is challenging. Several studies had visualized the oil/water displacement on a large scale, especially to show water coning (Safari et  . These experimental studies demonstrate the importance of studying two-phase flows adequately, but they were costly, and using a glass micromodel and glass-bead pack are more favourable for their accuracy in small scales and cost-efficiency. A fault/cross-bedding structure with permeability contrasts and also, in a crack between two layers with different permeabilities or a discontinuity between two different layers. These heterogeneous systems geologically exist on both small and large scales. The aim of this research is a precise evaluation of the preferential pathway of fluid flow by imbibition process (water flooding) in these structures. It enables us to visualize and investigate the effects of cross-bedding and discontinuity on fluid pathways and sweep efficiency for enhanced oil recovery purposes.

Glass micromodel fabrication
A micromodel for the evaluation of the fluid front in multiphase flow has been a beneficial tool in fluid mechanics. The main purposes of using a micromodel is the construction of any geometrical and structural patterns, direct observation of the fluid behaviour, and comparing the results with numerical simulations. There are several methods for the fabrication of a micromodel, including an etched glass micromodel (Homsy, 1987 Using a glass micromodel is the most common method for improved oil recovery. The current method, presented by Mohammadi et al., is preferable because of cost savings, fabrication safety, and the preservation of pore morphology and heterogeneity of the porous medium in comparison with other conventional methods (Mohammadi et al., 2013). In this study, a laser with 100 μm engraving resolution was used for etching the pore spaces on silicon glass. Porous media were designed using Corel-Draw software. The glass used is a common flute-type flat glass with high-temperature resistance and with 4 mm thickness. After cleaning, washing, and drying the engraved surface with a metal brush and water, the second flat piece of the glass was put on the etched glass to cover the engraved surface. At the inlet and outlet of the model, a penetration hole was drilled to allow the fluid to pass through the etched porous medium. Then, the injection needles were placed in these holes and were temporarily fixed with glue. The bonding process of two glasses is known as "fusing". It is accompanied by a gradual temperature increase inside the oven to the melting point temperature of the glass. So, the oven temperature was increased to 460°C with an incremental gradient of 5°C. Then the temperature was increased to 670°C with an incremental gradient of 1°C. After that, the oven was switched off, and the micromodel was gradually cooled in the oven for one day. Then, the needles were permanently fixed by glue, and accordingly, the model was ready to setup the experiment as shown in Figure 1. A syringe pump with 10 ml syringe, and K connector was used for injecting the fluid into the model. Injection speed was calculated at the inlet of the glass micromodel according to the cross-sectional area of the syringe.

Lattice Boltzmann Method
LBM can be obtained both by discretizing the Boltzmann equation, and Lattice Gas Automata (He and Luo, 1997). LBM has the capability to simulate the fluid flow in both small and large scales (Nabovati, 2009). It is more favourable due to studying viscous fingering on a small scale, reducing computational time, performing parallel computations, and tracking fluid front continuously (Chen and Doolen, 1998). Mean-field, free energy, colour gradient (Rothman-Keller), and pseudo-potential (Shan-Chen) are different models of LBM that are extended for two-phase/component fluid flow simulation (Gunstensen et al., 1991;Shan and Chen, 1993;He et al., 1998;Gonnella et al., 2007). In the current study, a pseudo-potential (Shan-Chen) model was used for its simplicity, computational time efficiency, and parallelism. In this method, density, viscosity, and surface tension are not independent variables (Sbragaglia et al., 2007). The Shan-Chen model is according to the intermolecular interactions of all fluids, in which the total force is equal to the summation of fluid-fluid intermolecular force, external force (like gravity), and fluid-solid adhesion force (Huang et al., 2015). (1) where: σ and ͞ σ -the symbols of two separate phases, F σ -total force on σ fluid (gr.lu/ts 2 ), F σ int -fluid-fluid intermolecular force on σ fluid (gr. lu/ts 2 ), F σ ext -external force on σ fluid (gr.lu/ts 2 ), F σ ads -fluid-solid adhesion force on σ fluid (gr.lu/ts 2 ), G σ͞ σ -the fluid-fluid particle interaction or cohesion parameter, G σ ads -the fluid-solid interaction parameter, G σ ads depends on the wettability. It is repulsive (nonwet) with a positive value or is attractive (wet) with a negative value. Shan and Chen introduced effective density (ψ σ (x,t)) as follows (Shan and Chen, 1993).

(5)
The total velocity (u σ eq ) and force of each phase (u total ) are applied by: These equations are used for the determination of velocity distribution functions in the equilibrium state (f i σ(eq) ) as follows: (7) where: c s -the velocity of sound speed in the lattice network (lu/ts), w i -the weights of the velocity distribution functions in the equilibrium state, i -node number, t -time (ts), f i σ -velocity components of each phase (dimensionless), f i σ(eq) -the equilibrium state of each phase (dimensionless), ρ σ -the density of each phase (gr/lu 3 ), u σ -the velocity of each phase (lu/ts), τ σ -the relaxation time (dimensionless). The procedure in two phases/components LBM is to consider a layer for each phase and in each node of the lattice, PDFs of phases are updated. After the initialization step the fluid variables are determined in a time marching process. This process will be stopped when the final criterion, such as breakthrough, is achieved.

Interfacial Tension and wettability Estimation
In the Shan-Chen method, the interfacial tension and wettability are determined by a series of simulations. Surface tension is dependent on the cohesion value of G coh and needs simulations for determination. The contact angle is dependent on the adhesion value of G ads .
The procedure of surface tension and contact angle determination is explained by Shiri et al. (Shiri et al., 2018). For the estimation of interfacial tension, droplets with different radii were placed in the centre of fully saturated container with another fluid. After equilibrium, the pressure inside and outside of the droplets were calculated. The slope of pressure differences versus the inverse of radii is interfacial tension. The estimated surface tension by the simulations for fluids with densities of 1 and 0.792, τ=1, and G coh =3 is 0.033 (see Figure 2). Estimations of contact angles with different adhesion parameters (G ads ) are shown in Figure 3. In this simulation, droplets were place on a solid surface, and after reaching equilibrium, contact angles were determined.

Grid Size Independent Test
The grid size independent test is done by investigating four lattice sizes. The displacing fluid (water) with a density of 1 and a relaxation time of 0.5555 is injected into porous media saturated with displaced fluid (oil) with a density of 0.875 and a relaxation time of 0.875. The injection rate is 200 (ml/h). The results for lattice sizes of 400×884, 600×1326, 800×1768, and 1000×2210 are shown in Figure 4. Lattices smaller than 600×1326 had some deviations. In Figure 4(e), the interpolated lines are water saturation (S w ) versus dimensionless time (t D =t/breakthrough), and no significant changes are seen in bigger sizes than 600×1326. So, lattice sizes above it do not affect the fluid flow pattern and sweep efficiency.

Results and Discussion
In immiscible fluid flows, the dominant force being applied to fluids is important. By considering the capillary number and the viscosity ratio, the flow pattern falls  In immiscible displacement with a low flow rate, the capillary force, as the main force, causes blob trapping. By increasing the injection rate, the viscous force will be dominant and consequently, this causes viscous fingering (Silva and Dawe, 2003;Holtzman, 2016). Displacing fluid viscosity over displaced fluid viscosity is known as dynamic viscosity ratio. If this ratio is smaller than one, the displacement is "unfavourable", and if it is bigger than 1, the displacement is "favourable" (Lenormand et al., 1988). Based on the range of dynamic viscosity ratio, three zones were defined. On the left-hand zone with low dynamic viscosity, the viscous force of the displaced fluid is the main force. In the middle zone or the transitional zone, the viscosity of both fluids is important. In the right-hand zone with high viscosity ratio, the displacing fluid has the dominant viscous force (Lenormand et al., 1988). In the current investigation, all simulations were in the middle zone.

Water Flooding in Fault Structure
Cross-beddings in the fluvial environment and faults are the most common structures. These structures are vital in petroleum reservoirs because the heterogeneity is not parallel to the flow direction (Weber, 1982; Ko- where: k 1 and k 2 are the permeability of first and second layers in millidarcy, respectively, θ 1 is the angle of between the streamline in the first layer and normal vector of interface, θ 2 is the angle of between the streamline in the second layer and normal vector of interface. If the angle of the streamlines is zero, the flow line will be normal to the boundary. If the angle of the streamlines is 90, the flow line will be parallel to the boundary. If the angle of the streamlines is between zero and 90, the flow line will be refracted and following the Equation 8. If k 1 > k 2 , the flow line will be refracted to reach normal vector of second media and if k 2 > k 1 , the flow line will be refracted to reach boundary between two media (see Figure  5) (Roti and Dawe, 1993;Bear, 2013).
As shown in Figure 6, the dimensions of the model were 2×1 cm. The interbedded layer has a 45 degree inclination, and its thickness is 15% of the model length.
In the left-hand inlet of the porous media, the Zou-He velocity boundary scheme is used (Zou and He, 1997), and open boundary condition is imposed at the righthand outlet (Mohamad, 2011). For the lower and upper boundaries and grains' surfaces, a half-way bounce-back boundary condition was used. In this section, porous media in all simulations were within a fault zone with different heterogeneities in permeability. They consist of two layers with a porosity of 40% and grain sizes of 673 μm (640-750 μm, grade 6) and 425 (μm) (310-425 μm, grade 9). They correspond to the permeability ratio of 2.5 based on Berg's equation (Berg, 1970). Initial conditions and fluids and solid properties were set based on the Dawe et al. study in 2011(Dawe et al., 2011. For wettability, the contact angle was assigned as 55°. The porous media were saturated with the oil with a viscosity of 1.4 centipoises (mPa s) and a density of 0.792 (g/ cm 3 ). Water was injected with a viscosity of 1 (mPa s) and a density of 1 (g/cm 3 ). The interfacial tension between two fluids was 34.1 (milli N/m). After changing

High Permeability Stripe
The results of water flooding in a high permeable fault zone are shown in Figure 7. During water flooding in a high permeability stripe in the early stage, water penetrates through the low permeability layer by a stable front. When the fluid front reaches the boundary between the low and high permeable layer, if displacing fluid pressure exceeds the capillary pressure between two layers, it enters the high permeable layer. So after filling layer 1, it passes through the shortest path of layer 2 with an angle of around 90 degrees. In the low permeable layer, oil is efficiently swept, while in the high permeable layer, oil is bypassed. W In a high permeable fault zone, water flooding obeys Hubbert's law (see

Low Permeability Stripe
The results of water flooding in a low permeable fault zone are shown in Figure 8. During water flooding in a low permeability stripe, water penetrates the high permeable layer with an unstable front. After reaching the low permeable layer (Layer 2), water completely sweeps the oil in this stripe. After filling this layer, water will bypass the next high permeable layer in the shortest pathway at 90 degrees. In a low permeable fault zone, water flooding obeys Hubbert's law (see Figure 5). Also, fluid front and flow pathway of LBM simulation were comparable with the experimental results of the previous study on a glass bead micromodel (Dawe et al., 2011).
As can be seen in Figure 7 and

Water Flooding in a Discontinuity
A fracture or discontinuity between two porous media consisting of two homogeneous layers with very low and high permeabilities is a hot topic for considering the effect of heterogeneity (Conn et al., 2014;Dong et al., 2017;Xiao et al., 2018). The following paragraphs are experimental and numerical evaluation of immiscible water flooding for sweep efficiency and streamlines tracking purposes in a fracture or discontinuity between two porous media.
Experimentally, a heterogeneous porous medium was fabricated using the method presented in the glass micromodel fabrication section. The experimental setup is shown in Figure 1. The glass micromodel had 1 inlet port and 1 outlet port at the centre of the model with 1.06 mm width. The porous medium had dimensions of 3.987 × 1.8 cm and the pore space height was 50 μm. The first layer has small pillars with a diameter of 450 μm, pore throat size of 120 μm, and 75% porosity. The second layer had large pillars with a diameter of 1350 μm, pore throats size of 360 μm, and 49% porosity. The size of the discontinuity or crack between two layers was 900 μm. The distances between the grains and walls was equal to pore throat sizes. The size of all inlet and outlet channels were 1200 μm. Arvandan crude oil from the south of Iran with a density of 0.875 (gr/cm 3 ) at 25°C, a dynamic viscosity of 5 (mPa s) at 25°C, and total acid number of 0.12 (mg KOH/g) was used for saturation of porous medium, and distilled water was used as the displacing fluid. Before starting the experiment, the glass micromodel was washed with deionized water and toluene. Then, it was dried by air blowing and exposing it to a high temperature. After that, the micromodel was completely saturated with the oil. Finally, distilled water was pumped into the micromodel at a constant flow rate.
The results of the experiment by discharges of 1 and 200 (ml/h) with equivalent capillary numbers of 0.07 and 14 are shown in Figure 9. In a low discharge of 1 (ml/h), the tendency of the fluid to pass through the pores is more than a higher discharge of 200 (ml/h). The crack or discontinuity between the two layers plays the  main role in the flow discharge through the porous medium. Fluids' tendency to pass through the fine and coarse layers decreases as the flow rate increases. An in-crease in flow rate from 1 and 200 (ml/h) causes a reduction in sweep efficiency from 14% to 7.3% by the moment of breakthrough, respectively (see Figure 11). Fluids and boundary conditions in the simulation were similar to the experiment. Zou-He velocity boundary scheme is used at both inlet and outlet of the porous medium (Zou and He, 1997) and, the half-way bounceback boundary condition is used for solid surfaces. After converting the flow units from SI to the lattice unit, the numerical method ran with discharges of 1 and 200 (ml/h) were equivalent to the capillary numbers of 0.07 and 14. The relaxation times of 0.875 and 0.5555 and the densities of 0.875 and 1 were set for the displaced fluid (oil) and displacing fluid (water), respectively. The results (see Figure 10) show fluid movement is in piston form in a very low capillary number. As a result, the increase in flow rate caused the fluid front piston pattern to change to the fingering form. Simultaneously, the crack has the main role in the fluid flow pathway in the heterogeneous porous medium. An increase in flow rate from 1 and 200 (ml/h) caused a reduction in the sweep efficiency from 15.6% to 10% at the moment of breakthrough, respectively (see Figure 11).
The results of the experimental and numerical method were consistent with hydrodynamic concepts. Sweep efficiency in both the micromodel water flooding and the LBM simulation were decreasing. The sweep efficiency curves of them versus flow rate were similar with a shift of about 2 percent (see Figure 11). So, a pseudo-potential model, which is used in this study, can efficiently track fluid front in a pore-scale in different heterogeneous systems. Also, an experimental glass micromodel is a sufficient, precise, cheap, and priceless tool for visualization of the fluid front compared to other experimental methods.

Conclusions
In this current study, the effects of capillary number on viscous fingering pattern in porous media with heterogeneity in permeability, such as a fault zone, and a discontinuity were studied. In the fault, the results showed that when the permeability of the fault zone is less or more than the vicinity layer, sweep efficiency, and fingering pattern of the fault zone and adjacent layer were different. In addition, based on Hubbert's law, the streaming direction depends on heterogeneity and deep angle.
In the discontinuity, the results showed that a crack between the two porous medium significantly contributed to the fluid flow, and the sweep efficiency decreases as the flow rate increases. The LBM and glass micromodel results showed acceptable consistencies with the physics of fluid flow. They showed that a pseudo-potential model has a good ability and accuracy in the simulation of immiscible multi-component fluid displacement. Also, a glass micromodel is an efficient tool to fabricate any structure and pattern of heterogeneous porous media, and visualize fluid flow. It is recommended that more glass micromodels with heterogeneity in both permeability and wettability should be studied in future works, and also, LBM simulation is enhanced for solving real-life problems.