Numerical Investigation of Fluid Flow Instabilities in Pore-scale with Heterogeneities in Permeability and Wettability

Quadrant geometry with permeability and wettability contrast occurs in different events, such as faults, wellbore dam age, and perforation zones. In these events, understanding the dynamics of immiscible fluid displacement is vital for enhanced oil recovery. Fluid flow studies showed that viscous fingering occurs due to viscous instabilities that depend on the mobility of fluids and capillary forces. Besides, the porous domain heterogeneity is also effective on the formation of fingering. So, the purpose of the current research is to numerically investigate the effect of heterogeneity in wettability and permeability, and flow properties in Saffmann-Taylor instabilities. Numerical simulations with different flow rates in the permeability contrast model illustrated the nodal crossflow, growth of viscous fingering in the nodal part, and by-pass flow in the second zone. In the wettability contrast model, a capillary fingering pattern is observed and fluid patch es are isolated because of capillary force and the end effects are trapped within the quadrant. Moreover, the consequenc es of wettability on apparent wettability that alters the fluid-front pattern and displacement efficiency are shown. effects; Lattice Boltzmann method.


Introduction
Understanding the true fluid movement path and displacement pattern is necessary for efficient hydrocarbon reservoir production management (Dawe, 1998). Reservoir water/oil displacement in heterogeneous porous media has been investigated in large and pore-scales, but fewer studies are done in mesoscale, i.e. centimetre-scale (Kryukov et  A fault is a narrow discontinuity with shear displacement in rock. Based on the magnitude of three principal stresses, three tectonic regimes of normal, thrust and strike-slip faults were classified (Anderson, 1951; Fossen, 2020). Normal faults occur in extensional tectonic systems and usually are formed with two shear breaks in rock with angles of around 60 o . Cross-cutting (conjugate), "X" shape, or "hourglass" fault structures are common normal fault systems that occur in any scales ranging from centimetres to kilometres (Bretan et al., 1996;Ferrill et al., 2000Ferrill et al., , 2009. The quadrant geometry with permeability and wettability contrast occurs in different events such as faults, wellbore damages, and perforation operations (McCarthy, 1991). Conjugate faults are those events that brings different geological strata towards each other and make heterogeneous porous media. Heterogeneity in wettability and permeability depends on rock properties and changes the distribution of fluid saturations. Heterogeneity in permeability is the spatial variation of permeability from one point to another and heterogeneity in wettability is the spatial variation of wettability from one point to another. Faults, especially conjugate faults, are the events that create these contrasts: (1) permeability heterogeneity by adjoining a high permeable sandstone and a low permeable shale and (2) wettability heterogeneity by adjoining oilwet limestone and water-wet silica clastic rocks (Woods, 2015). Taiyuan Formations in the north eastern Ordos Basin and Longtan Formation in south-western Guizhou are two cases of conjugate faults with quadrant structure in gas fields having heterogeneities in both wettability and permeability. This field is developed in swamps/lagoons in the delta plain and coastal environments. Lithology and sedimentary facies, pore spaces, structural geometries, and long-term tectonic movements are the main factors of shale gas enrichment, hydrocarbon migration, and entrapment (Wang and Guo, 2020). Paradox Basin is another example of silica-clastic rocks contiguous to chemical rocks (e.g. limestone) and evaporates. In this basin, a rapid subsidence, a periodic sea level change, and a high evaporation rate led to the formation of 3-km thick sediments consisting of dolomite, evaporate, and black organic shale. Salt-core anticline and Paradox fault-fold belt are the two main tectonic active systems in this region (Torabi et al., 2019). Asmari Formation is another case of adjoining silica-clastic rocks and chemical rocks. Enormous amounts of folding and faulting of nearly 11000 meters of sediments due to the subduction of the Arabian Plate underneath central Iran's plate, along with the effects of salt diapirism, made this region the most hydrocarbon-rich province of Iran. Heterogeneity in permeability and heterogeneity in wettability are the most common events in this highly active tectonic system. In this region, Dezful Embayment with an area of only 60,000 km 2 is one of the world's richest oil provinces containing 8% of global oil reserves. The lithology of the Asmari Formation consists of highly fractured carbonate rocks, loose terrigenous sandstone, interbedded shale and marl, and some evaporate sediments of anhydrite and gypsum (Bordenave and Hegre, 2005;Haidari et al., 2020).
One ideal geometry of nodal flow is the quadrant structure, in which the end effect boundaries cause fluid movement to flow from a low capillary pressure region to a high capillary pressure region. This process leads to the effect that is correlated with heterogeneity in wettability and permeability. Nodal flow in quadrant geometry is one of the structures of fluid flow mechanisms in this geological event that has been a popular research topic (Caruana and Dawe, 1996 (Caruana and Dawe, 1996). Silva et al. and Dawe et al. experimentally studied the effects of wettability and permeability heterogeneity in both miscible and immiscible displacement in a quadrant structure (Silva and Dawe, 2003;Dawe and Grattoni, 2008). The results were in agreement with Yeo and Zimmermann's study and showed that different fluid-fronts will form in both miscible and immiscible displacement (Yeo and Zimmerman, 2001). Also, Al-Hadhrami et al. experimentally investigated miscible fluid flow through a quadrant model and a porous medium with two discontinuous shale barriers. The results showed the solvent moved preferentially through high permeable layers regardless of the amount of permeability. So, the flow path is controlled by the relative permeability rather than the absolute permeability (Al-Hadhrami et al., 2014). Moreover, Gomes et al. modelled multi-phase flow in a quadrant model by a force-balanced control volume finite element method. Their results were in good agreement with experimental results and approved their numerical approach (Gomes et al., 2017).
Fluid flow in a porous medium is a function of both local flow conditions and wettability of solid surfaces. Wettability is the interaction of two fluids to maintain contact with a solid surface. This special characteristic of solid surfaces plays an important role in porous media, such as oil production, soil mechanics, coating, etc. Wettability has a vital role in immiscible displacement efficiency (Lenormand et  Rabbani et al. studied the effects of pore geometry on apparent wettability. They defined a new apparent wettability number based on the chemical properties of a solid surface and pore geometry. In addition, they numerically showed that the interface curvature of fluid is strongly affected by pore geometry (Rabbani et al., 2018). Avendaño et al. studied water flooding in waterwet and oil-wet porous networks constructed by dolomite microfluidics. The results showed uniform fluidfront displacement and cooperation of capillary force for pore filling phenomenon in a water-wet network (Avendaño et al., 2019).
When it is not possible to analytically solve fluid flow, experimental and numerical simulations are the solutions. The most common tasks for numerical simulation are discretizing the flow equation by finite element methods and solving them in macroscopic scales. In these methods, the interfaces of fluid can be obtained, but no data in the porescale is available. Another solution is the Lattice Boltzmann Method (LBM), a relatively new useful mesoscopic numerical modelling in fluid dynamics (Mohamad, 2011). Using various boundary conditions and capability of simulating any complex porous media make this method favourable for numerical simulations of heat transfer and fluid flow (Yilbas et al., 2017;Filatov and Yakimov, 2018;Zakirov et al., 2018).
Therefore, the main purpose of the current work is to show viscous instability numerically in a quadrant structure with wettability and permeability heterogeneity in a mesoscale by LBM, and then to compare the results with previous experimental methods. The two main goals of the study are: the effect of flow rate on flow instability in quadrant model with permeability contrast (section 3.1), and a close interaction of viscous force and capillary

Fluid Flow through Porous Media and Viscous Instabilities
Estimating the initial saturation of oil and gas and the production efficiency are challenging for petroleum reservoir engineers. Hydrocarbon production efficiency is controlled by combinational forces of viscosity, capillary, gravity, etc. The effects of these forces depend on wettability, interfacial tension, viscosity, density, saturation, pressure changes, porosity, and permeability (Amyx et al., 1960). In hydrocarbon reservoirs, there are several immiscible fluids within pore spaces and therefore, capillary pressure is one of the main forces in hydrocarbon recovery (Dawe, 1998;Nguyen et al., 1999). In immiscible water/oil displacement, imbibition is the result of an increase in wetting fluid saturation, and drainage is the result of an increase in the non-wetting fluid saturation. Spontaneous imbibition happens when a wetting fluid displaces a non-wetting fluid only with capillary forces (Towler et al., 2017).
There are three types of displacement processes: (1) miscible with the first contact, in which two fluids are mixed in all their properties; (2) miscible with multiple contacts, in which two fluids are not mixed simultaneously and miscibility happens by mass transfer during the process; and (3) immiscible displacement, in which two fluids are separated by an interface and there is no chemical interaction between them (Kantzas et al., 2015). The mobility of the fluids is usually defined as a permeability/dynamic viscosity ratio, so the units will be m 2 /(Pa×s). The mobility ratio (M), according to the following formula, is the ratio of the mobility of a displacing fluid to the mobility of the displaced fluid. If M>1, then the flow is unfavourable as the displacing fluid is more mobile than the displaced fluid. It is identified by an early breakthrough, low recovery factor, and viscous fingering pattern. If M<1, then the flow is favourable and in this condition, if the viscous force is dominant, then piston-like displacement occurs. On the other hand, if the capillary force is dominant, then the flow regime is capillary fingering. For miscible displacement, in which fluids are mixed, the mobility ratio is the ratio of the viscosity of the displaced fluid over the viscosity of the displacing fluid (Lenormand et al., 1988). In the current study, the simulations of section 3.1 are in the viscous fingering region and the simulations of section 3.2 are in the stable front and capillary fingering regions.

Quadrant Model
The quadrant or checkerboard model is an ideal 2×2 block geometry, in which any section or region has its own properties to illustrate the effects of heterogeneities caused by faulting (see Figure 1). As stated in the introduction, heterogeneities occur in both permeability and wettability contrasts, e.g. the Asmari Formation in the Dezful Embayment in the south-west of Iran (Haidari et al., 2020). We investigated two kinds of heterogeneities in permeability and wettability models in sections 3.1 and 3.2, respectively. In section 3.1, we assumed sandstone in cell 1 and shale in cell 2 to investigate heterogeneity in permeability. In section 3.2, we assumed waterwet silicaclastic rock in cell 1 and oil-wet carbonate rock in cell 2 to investigate heterogeneity in wettability. Other assumptions of the quadrant models are as follows: a solid boundary condition is set at the top, and bottom boundaries and flow boundary conditions are set at both the left-inlet and the right-outlet of the model (see Figure 1). Initially all porous media are saturated with oil and then the water imbibes through the porous media. To consider the effect of spontaneous imbibition in all simulations, both the inlet and outlet of the models are subjected to water contact with an initial confined pressure. The inlet and outlet pressure gradient of the models are different in each of the simulations.
The checkerboard geometry usually exists in petroleum reservoirs as a combination of layering and crossbedding. The challenging issue for petroleum engineers is where the injected fluid tends to flow. This model showed heterogeneity in permeability and wettability in many cases in petroleum engineering: layers are subjected to faulting; the mud filtrate is invaded through the well-bore (in an over-balanced drilling system); the remedial fluids are injected into the formation by acidification; and permeability changes by any fracturing in the damage zone (perforation operation). wettability. In this simulation, cell 1 and cell 2 were sandstone (35% porosity) and shale (45% porosity) with grain sizes of 1 mm and 0.1 mm and a permeability of 100 and 4 mD (9.87e-11 m 2 and 3.95e-12 m 2 ), respectively. The second simulation was run for an investigation of water flooding in a quadrant model with wettability contrasts. In this simulation, cell 1 and cell 2 had an equal grain size of 640 μm (fine grains limestone and clastic rocks with 35% porosity) with contact angles of 60 and 120 degrees, respectively. In LBM, the contact angle is control by the fluid-solid interaction coefficients and set by the simulation of a droplet on a solid surface. These sizes and conditions were chosen to compare the results of the simulation with the previous experimental results.

Numerical Method
LBM originates from Lattice Gas Automata and can be derived from discretizing the Boltzmann equation. Several LBMs exist and in this study, we use a two-component single-phase pseudo-potential method with single relaxation time or Bhatnagar Gross and Krook (BGK) approximation. In this method at any lattice node, each fluid acts as a layer of the distribution function and the energy is transferred by the collisions and streaming operations. A detailed description of the LBM model is explained by Shan and Chen. The distribution functions for each component obey the Boltzmann equation as follows (Shan and Chen, 1993): (1) where: i -index of velocity components, k -two fluid components of a and a', -the identity vector that indicates the velocity direction of different components in a lattice node, τ k =1 -the relation time for each component to achieve an equilibrium state, f k(eq) i (x,t) -the equilibrium state of each fluid distribution. The equilibrium state of each fluid distribution for recovering Navior-Stokes equation as follows (Shan and Chen, 1993): (2) In this study, we use a D2Q9 model with discrete velocities (e i ) and the weights of the velocity distribution functions in the equilibrium state (w i ) as follows (Shan and Chen, 1993): where: a, b, c, d -lattice constant, c = Δx / Δt (lu / ts) -the ratio of lattice spacing (in this study Δx=1 lu and time step of Δt=1 ts), F k -total force on k fluid (g.lu/ts 2 ), i -node number, t -time (ts), ρ k -the density of each phase (gr/lu 3 ), u ' -the overall velocity (lu/ts), u k(eq) -the equilibrium velocity of each phase (lu/ts). In all simulations, LBM models are D2Q9 and lattice dimensions are 500×1000. A grid size independent test showed these lattice dimensions were suitable for the simulations. The porous media were initially saturated with oil with a viscosity of 0.0015 Pa×s (1.5 cP) and a density of 0.8 g.cm -3 . Moreover, water was injected with a viscosity of 0.001 Pa×s (1 cP) and a density of 1 g.cm -3 . G coh is 3 for immiscible displacement and the interfacial tension between the two fluids is 0.033 (N m -1 ). All of these parameters are in SI units that are then converted to lattice units. Time in lattice units is time step (ts). Porous media were saturated with oil and the left and right boundary conditions were subjected to water contact with an initial confined pressure of 1 mu. lu -3 . Constant pressure gradients of 0.2, 1, 0.1, and 0 mu.lu -3 proposed by Huang (Huang, 2016) were applied at the left-inlet and the right-outlet of

Water Flooding in Quadrant Model
with Permeability Contrast Figure 2 and Figure 3 show water flooding stages with different flow rates in an oil-saturated quadrant model with heterogeneity in permeability. In this model, shale and sandstone rocks are placed in the quadrant geometry. In these two models, pressure gradients are different and both of them are in the transition zone between viscous fingering and stable front. The outcome of the current simulation is as follows: • for each case, fluid initially entered both low and high permeable sections of the left side but after a while, the displacing fluid moved faster in the high permeable section. It is due to boundary pressures and diffusivities that are caused by large permeability gradients. The low adjacent permeability region has low saturation or even a stagnant gradient. So, the preferential flow pathway is a high permeable region and some part of other low permeable regions remain intact; • after fully saturating the first high permeable region, the displacing fluid moved toward the node at the centre of the quadrant model and entered the second high permeable section. In this stage, depending on the flow rates and the cross-sectional area of the nodal part, oil was replaced by the displacing fluid with different fingering patterns and different degrees of saturation around the node; • after water passed from the node and entered the second high permeable section, it bypassed towards the outlet and only pushed a portion of saturated oil based on the flow rate and the permeability; • after the breakthrough, the injected fluid and its streamline came to the steady-state condition and oil and water saturations in the model remained constant. For more information, the finger's formation and its growth are sensitive to heterogeneity that is induced by permeability contrast and also, the pressure gradient that increases the flow rate. The pressure gradient value has three effects on the fluid flow that can be compared in Figure 2 and Figure 3: (1) nodal cross-sectional saturation (the solid boundaries at the centre of the model) is dependent on the flow rate; (2) The viscous fingering pattern is more obvious when the pressure gradient is high; (3) Recovery efficiency is reduced when the flow rate increases. Figure 4(a) shows the conceptual quadrant model with heterogeneity in permeability.  Table 1.
In the study of Silva et al., intermittent water and oil injections were performed experimentally on a sand pack (Silva and Dawe, 2003). As shown in Figure 4(c), oil almost equally moves the water in two media with different permeability. Heterogeneity in permeability has little effect on the direction of the injecting fluid. It does not show the tendency of the fluid flow in a preferential path way in a higher permeability layer.
In the study of Dawe et al., oil was displaced by water experimentally on a sand pack (Dawe and Grattoni, 2008). As shown in Figure 4(d), the tongue phenomenon is obvious in fluid displacement. Nodal flow occurs at the centre of the model. The channelling of fluid flow is clear when passing through the central node and then in the more permeable cell. Also, there is fluid displacement in a low permeable cell, but it is slower.
In the studies of Gomes et al. and Christou et al., the displacement of oil by water has been done numerically using the CV-FEM method (Gomes et  The results of the present study are in good agreement with the results of the studies of Gomes et al., Dawe et al., and Christou et al. Considering that the finger phenomenon is affected by both fluid properties, and heterogeneity of porous medium, the slight differences in the results of these 4 models are due to the difference in the amount of these factors, which are expressed in Table 1.

Water Flooding in Quadrant Model with Wettability Contrast
In the current simulation by instantaneous imbibition, we can evaluate capillary force and wettability heterogeneity effects on production performance. Like the con-ceptual model of Figure 1, in a quadrant with wettability contrast, the lower-left and upper-right sections of the quadrant model are considered water-wet regions with a contact angle of 60 o . The upper-left and lower-right sections of the model are considered oil-wet with a contact angle of 120 o . In this model, carbonate and fine grains clastic rocks are placed in the quadrant geometry. All sections have the same permeability and porosity. For a better understanding of the effects of wettability and instantaneous imbibition, in the first model (see Figure 5), the pressure gradient is too small and in the second model (see Figure 7), it is zero.
In Figure 5, the constant low pressure gradient across the porous medium and wettability are the two forces that affect fluid flow. We discuss five stages of flow: the initial stage (see Figure 5(b)), the first cell saturation stage (Figure 5(c)), the nodal flow stage (Figure 5(d)), the breakthrough stage (Figure 5(e)), and the final saturation distribution (Figure 5(f)).
• In the initial stage, water starts to push oil faster in the lower-left water-wet region by the imbibition process. In the upper-left region, water enters and then moves downward toward the lower-left waterwet region. These conditions are fixed as water is being pumped into the left part of the model. • After that, in the first cell saturation stage, water in the lower-left region pushes oil until it reaches the middle boundaries. Then, the fluid-front is inclined towards the node at the centre of the model, and water is trapped in the lower-left section. • In the nodal flow stage, by filling the lower-left water-wet region, water overwhelms the capillary pressure contrast in the nodal part. The difference of capillary pressure in the nodal part is smaller than the difference of capillary pressure between oil-wet and water-wet sections' boundaries. Therefore, water enters the upper-right water-wet section from the nodal part and flows toward the right-outlet boundary for the breakthrough. As it is obvious in this case, the preferential flow pathway is a channel between two cells because of the simultaneous action of wettability force and pressure gradient. • At the breakthrough, the second water-wet region is not fully saturated and the waterfront moves in the channel.
• Final saturation distribution shows entrapment of three oil patches in porous medium and also, the tendency of oil for leaving out from the right-centre of porous medium. Figure 6(a) shows the conceptual quadrant model with heterogeneity in wettability. Figure 6(b) shows the results of the present study, which are presented here for comparison with other studies, and its explanations are provided in the previous paragraphs. The results of this section is comparable with the numerical and experimental results of Silva et al. and Dawe et al. (Silva and Dawe, 2003;Dawe and Grattoni, 2008). The fluids and solid properties are listed in Table 2.
In the study of Silva et al., water was displaced by oil experimentally on a sand pack. In this experiment, cell 1   was oil-wet and the oil flow entered the model from the left at a very low speed (Silva and Dawe, 2003). As shown in Figure 6 (c), the preferred oil path is cell 1 (oil-wet region). The oil entered the water-wet region at the entrance due to the viscous force, but its movement in this region stopped and then deviated towards the oilwet region. In the nodal part, viscous force is greater than the nodal resistance and capillary pressure and the oil in the oil-wet region moves towards the output of the model. In the study of Dawe et al., oil was displaced by water experimentally on a sand pack (Dawe and Grattoni, 2008). As shown in Figure 6(d), the predominant force in oil movement is the capillary force, and the viscous force has no effect on the displacement of oil from the water-wet region.
The results of the present study are in good agreement with the results of Silva et al. study. The pattern and flu-id displacement path showed simultaneous roles of both capillary and viscous forces in fluid flow in porous medium with heterogeneity in wettability.
In Figure 7, only wettability contrast exists in the model, and oil-saturated sandstone is surrounded by water with an equal amount of pressure in both the left-inlet and the right-outlet boundaries. As a result, we can see three stages of flow: the initial stage (see Figure 7(c)), first water-wet cell saturation stage and nodal stage (see Figure  7(d)), and the final saturation stage (see Figure 7(e)).
• In the initial stage (see Figure 7(c)), we observed the effects of wettability and capillary force on infiltration and saturation of oil-wet and water-wet cells. In a water-wet region, water starts to push oil, and in the oil-wet regions, water enters the large pore spaces among the grains and a thin film of oil remains on the grains' surfaces. This film of oil plays the main role for oil outflow. These conditions of 50,000, c) saturation at time steps of 2,700,000, d) saturation at time steps of 3,700,000, e) saturation at time steps of 6,700,000, e) final saturation distribution, the blue colour is water, the yellow colour is oil, and the black colour is solid.
are fixed until the waterfront in the oil-wet regions are connected to the waterfront in the water-wet sections. Then, water leaves from the large pore spaces of oil-wet sections and moves backward towards the water-wet sections. • In the first water-wet cell saturation stage (see Figure 7(d)), the instantaneous imbibition process of the lower-left region is stopped due to weaker capillary force in comparison of the upper-right waterwet region. The wettability force of the upper-right water-wet region pushes oil until the displacing fluid (water) fully saturates the region. After that, fluid-front in the lower-left region starts to push oil until this region becomes fully saturated. This is due to the capillary pressure contrast in each step.
In the nodal part, the node acts as a barrier and does Water-wet 60 ----Oil-wet 120 ---- Figure 6: Comparisons of fluid displacement in porous media with heterogeneity in wettability: (a) conceptual model, (b) numerical (LBM) results of the current study, the blue colour is water, the yellow colour is oil, and the black colour is solid (c) adjusted experimental results of Silva et al. study (2003) not allow the fluid-front from the upper-right waterwet section to pass through this area. • Finally, in the final saturation stage (see Figure  7(e)), water wet sections become fully saturated with water and there is a little invasion of water through the oil-wet sections. The instantaneous imbibition simulation time in LBM was very long and had 40,000,000 time steps. Due to the finite velocity speed of oil exit from the oil-wet region, water preferentially enters from the right or left of cell 1. As can be seen, the two factors of the node at the centre of the model, and very small differences in particle sizes and pore throats, cause a change in the saturation order of the two seemingly similar cells. Performing this simulation, which shows the importance of the central node and the effect of minor heterogeneities

Effects of grain geometry on fluid front shape
Grain geometry affects the fluid front shape and its apparent wettability. Dimensionless apparent wettability number is introduced as (Rabbani et al., 2018): where: θ -static contact angle, β -capillary angle.
In this equation, capillary angle (β) depends on the pore throat geometry. If W>1, then flow is in strict drainage and the fluid-front is convex both before and after passing through the throats. If W<-1, then flow is in strict imbibition and the fluid-front of both inlet and outlet of the throats are concave. If -1<W<1, then the pore geometry controls the apparent wettability and therefore, concave and convex interfaces exist simultaneous-ly. This is a usual phenomenon in contact angle between 60 o to 120 o (Rabbani et al., 2018). In our study, static contact angles are 60 o for water-wet regions and 120 o for oil-wet regions. β is 30 o for the porous medium of the current study that it is derived from: (8) where: r b -pore body radius (lu), r t -pore throat radius (lu).
So, in the current study, the apparent wettability numbers are between -1 and 1 (-1≤W≤1). As seen in Figure  5(a) and Figure 7(b), apparent wettability in the oil/ water-wet regions shows a concave fluid-front before passing through the pore throats and convex fluid-front in the front of pore throats. The flow rate does not have a significant influence on the apparent wettability and these fluid-front curvatures are formed as a result of pore geometry and contact angle. Apparent wettability is related to the capillary angle (β) which depends on the pore throat geometry. The capillary angle varied from the centre of the throats so that the apparent wettability changed (see Figure 5(a) and Figure 7(b)).

Conclusions
Density and viscosity ratio, and in addition, heterogeneity are three factors that are important in immiscible fluid flow and fingering formation. In the current study, the effects of wettability and permeability heterogeneities in fluid front and preferential flow path way were investigated numerically by LBM and compared with previous experimental results. The quadrant model is an ideal model for fluid flow instability investigation in heterogeneous porous media. The main new outcome of the current study can be summarized as follows: • the first evaluation studied the effects of permeability heterogeneity on fluid front instability. The results showed that fingering formation is dependent on fluid properties and pore throat geometry of a Figure 7: Water saturation (spontaneous imbibition process) of the quadrant model with heterogeneity in wettability (4 cells with different wettability a) conceptual model, b) effects of grain geometry and wettability on fluid front shape, c) saturation at time of 10,0000, d) saturation at time steps of 22,000,000, e) saturation at time steps of 40,000,000, blue colour is water, yellow colour is oil, and black colour is solid  porous medium simultaneously. A little difference in each property may weaken or intensify the fingering formation; • the second evaluation evaluated the effects of wettability heterogeneity on fluid front instability. Wettability contrast between all sections of the quadrant model causes oil patches to become trapped in oilwet regions. This phenomenon reduces the displacement efficiency by the capillary entrapment mechanism. The displacing fluid (water) in the wettability contrast model moves preferentially forward, backward, and bilateral. It depends on the amount of wettability force to absorb water molecules through the regions with more capillary pressure, and it also depends on the strength of viscous force. In instantaneous imbibition, pore throat size, grain geometry and wettability are vital for the preferential flow pathway of the displacing fluid; • grain and pore throat geometries affect the apparent wettability, and fluid front pattern. Here, models with a contact angle between 60 o and 120 o showed a co-existence of a concave and a convex fluid-front in the inlet and the outlet of the pore throats.