Influence of Wetting on Viscous Fingering Via 2D Lattice Boltzmann Simulations

We present simulations of two-phase flow using the Rothman and Keller colour gradient Lattice Boltzmann method to study viscous fingering when a “red fluid” invades a porous model initially filled with a “blue” fluid with different viscosity. We conducted eleven suites of 81 numerical experiments totalling 891 simulations, where each suite had a different random realization of the porous model and spanned viscosity ratios in the range M∈[0.01,100]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M\in [0.01,100]$$\end{document} and wetting angles in the range θw∈[180∘,0∘]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _w\in [180^\circ ,0^\circ ]$$\end{document} to allow us to study the effect of these parameters on the fluid-displacement morphology and saturation at breakthrough (sweep). Although sweep often increased with wettability, this was not always so and the sweep phase space landscape, defined as the difference in saturation at a given wetting angle relative to saturation for the non-wetting case, had hills, ridges and valleys. At low viscosity ratios, flow at breakthrough is localized through narrow fingers that span the model. After breakthrough, the flow field continues to evolve and the saturation continues to increase albeit at a reduced rate, and eventually exceeds 90% for both non-wetting and wetting cases. The existence of a complicated sweep phase space at breakthrough, and continued post-breakthrough evolution suggests the hydrodynamics and sweep is a complicated function of wetting angle, viscosity ratio and time, which has major potential implications to Enhanced Oil Recovery by water flooding, and hence, on estimates of global oil reserves. Validation of these results via experiments is required to ensure they translate to field studies.


Introduction
In this paper, we study the physics and morphology of immiscible two-phase fluid flow in a 2D model of a porous medium with the aim of increasing understanding of patterns fluid-fluid displacement and their effect on the fluid saturation, and in particular, the effect of wettability, a measure of the degree to which a fluid is attracted to the solid.
We will make use of the Lattice Boltzmann Method (LBM), a well-known method to simulate fluid flow, which is capable of modelling flow of immiscible fluids at the pore scale. Lattice Boltzmann Methods have their origins in Lattice Gas Automata (LGA) in which particles move and collide on a discrete lattice representing a simplified discrete version of molecules moving and colliding in a gas. LGA were first proven by Frisch et al. (1986) to yield the Navier-Stokes equations in the macroscopic limit. The LGA lead to socalled Lattice Boltzmann Methods where one is solving the classical Boltzmann equation on a discrete lattice involving number densities moving and colliding on a discrete lattice. Since an efficient method via relaxation to calculate the collision term due to Bhatnagar, Gross and Krook-the BGK method (Bhatnagar et al. 1954)-was developed (Qian et al. 1992;Chen and Doolen 1998), research and applications of the Lattice Boltzmann Method have undergone an explosion (see Succi 2001Succi , 2018. The Lattice Boltzmann Method has been applied to various studies of viscous fingering. Dong et al. (2010) used the Shan and Chen Lattice Boltzmann Method (Shan and Chen 1993) to simulate viscous fingering in a channel and investigate the effects of wetting and viscous fingering in a Hele-Shaw cell and simple porous media (Dong et al. 2011). The Rothman-Keller (RK) colour gradient multiphase LBM has been applied to study the mode of flow (viscous fingering, capillary fingering and stable displacement) in the phase space of capillary number versus viscosity ratio by Huang et al. (2014). A forcing term Lattice Boltzmann approach has been applied to simulate viscous fingering where the porous medium was mimicked using the grey lattice Boltzmann or the Brinkman force model (Vienne et al. 2019). And the RK LBM was applied to the study of immiscible displacement by a shear-thinning fluid in a porous medium (Wang et al. 2019). In 2019 ), a benchmark study was conducted where 14 teams compared the results of using a range of different pore-scale flow simulation methods including the LBM method in 2D and 3D, and found that no single method excels across all conditions and that corner flow and thin films present significant challenges, although LBM methods seemed to capture such effects at least to some extent. In this paper, we make use of the Rothman-Keller colour gradient Lattice Boltzmann Method (Rothman and Keller 1988;Latva-Kokko and Rothman 2005;Reis and Phillips 2007) to study viscous fingering due to its ability to simulate large viscosity ratios and ease and accuracy of simulating different wettability. The main pitfall of the RK LBM is that the colour gradient required to calculate the effect of fluid cohesion has some numerical errors which leads to anisotropy of the interfacial tension and spurious currents. This can affect details of phenomenology such as droplet and finger formation but the use of more accurate and isotopic colour gradient calculations renders these inaccuracies negligible (Mora et al. 2021a).
It is well established that when a low viscosity fluid such as water is injected into a porous rock matrix saturated with a higher viscosity fluid such as crude oil, patterns of viscous fingering occur (Homsy 1997;Måløy et al. 1985;Chen and Wilkinson 1985;Lenormand et al. 1988). Namely, the formation of patterns at the unstable interface between the two fluids in the porous medium (typically, patterns of narrow tendrils of the injected fluid in the porous medium). This is an interesting physical phenomenon with practical implications to Enhanced Oil Recovery (EOR) where water is injected into one well to help extract the oil from the hydrocarbon reservoir at an adjacent production well. The displacement morphology of immiscible fluid displacement in a porous medium is complex and is affected by many parameters of the fluids including the viscosity ratio of the two fluids and the wettability. Wettability refers to the relative affinity of the solid for one of the fluids with respect to the other. Thus, we say that a porous medium is wetting to water rather than oil if it exhibits a tendency to be coated by water and to repel oil. The angle of the fluid interface relative to the solid surface is called the wetting angle. A wetting angle approaching 180 • means that a droplet is almost circular in 2D (spherical in 3D), and such fluids are termed non-wetting. The opposite case is a highly wetting fluid where the fluid droplet is more attracted to the solid and the angle between the fluid interface and solid approaches 0 • . Our aim here is to apply the Lattice Boltzmann Method to the study of immiscible flow patterns such as viscous fingering in a model porous medium in 2D, and in particular, to study the effect of the wetting angle on the sweep, which is well known in petroleum engineering to play a vital role in determining the "recovery factor"-the total fraction of oil in a given hydrocarbon reservoir that can be produced (Deng et al. 2020). Numerous studies in the field of petroleum engineering have led to a broad consensus that the recovery factor can be increased by increasing the wettability of the invading fluid (Kennedy et al. 1955;Jadhunandan and Morrow 1995;Seethepalli et al. 2004;Morrow and Buckley 2011;Sharma and Mohanty 2013).
There have been many studies, experimental, theoretical and numerical on the influence of wettability on the pattern of flow in porous media. Stokes et al. (1986) found in experimental work that wettability affects the finger width. Cieplak and Robbins (1988);1990 (Cieplak andRobbins 1990) developed a numerical model for quasi-static fluid-fluid displacement and conducted numerical experiments of flow in a 2D array of discs showing that as the wetting angle decreases, a progressive smoothing mechanism occurs and that the width of invading fingers seems to diverge beyond a critical angle which depends on porosity. Trojer et al. (2015) conducted a systematic experimental study of fluid-fluid displacement in a granular pack and found that wettability profoundly affects the invasion morphology. Namely, they observed a compactification of viscous fingering and a regime of compact displacement at low capillary numbers for weak imbibition. Zhao et al. (2016) found in microfluidic experiments involving vertical posts representing a porous medium, that as wettability is increased, there is more efficient displacement and higher saturation up until a critical angle is reached, after which, the system undergoes a wetting transition and the trend is reversed. Other microfluidic research combined with theoretical analysis and pore-scale simulations (Hu et al. 2019;Lan et al. 2020) have studied the phase diagram capturing the viscous to capillary fingering transition to study the impact of medium disorder and wettability on this transition. Research using an invasion-percolation model by Primkulov et al. (2018) extended the Cielpak and Robbins description of quasistatic fluid invasion reproducing the wetting transition in strong imbibition. There have been many core-scale experiments which show that the displacement efficiency can be improved when the wettability is increased towards imbibition such as by addition of surfactants or by using low-salinity water flooding (Kennedy et al. 1955;Seethepalli et al. 2004;Morrow and Buckley 2011;Sharma and Mohanty 2013).
In the following, we apply the Rothman and Keller colour gradient multiphase Lattice Boltzmann Method-RK colour gradient LBM (Latva-Kokko and Rothman 2005) with the aim of studying the effect of wetting on morphology and efficiency of flow of an injected lower viscosity fluid into a porous rock matrix saturated with a higher viscosity fluid. Our reason to choose the RK LBM is its ability to handle large viscosity ratios, a wide rangeover ten orders of magnitude-of surface tensions, and its accuracy as well as convenience of setting the wetting angle.
In previous work, we found that the effect wettability on saturation at breakthrough was complex and demonstrated that the saturation did not necessarily increase with wettability, and that optimal wetting angles w > 0 that maximized saturation could occur at specific viscosity ratios (Mora et al. 2021b). In this work, we extend the previous work and aim to study the effects of wetting angle and pore matrix geometry on the morphology of viscous fingering and the impact of wettability on saturation at and beyond breakthrough, and to improve understanding of the patterns of viscous fingering and its effect on sweep.

Numerical Simulation Methodology
In this paper, we apply the Rothman-Keller (RK) multiphase Lattice Boltzmann model which was originally derived for a Lattice Gas Automaton (Rothman and Keller 1988;Gunstensen et al. 1991), and later extended to the Lattice Boltzmann Method (LBM) by Latva and Kokko (2005). The colour gradient RK Lattice Boltzmann model involves modelling particle distributions denoted f k of two fluids (red and blue for k = 1 and k = 2 ) moving and colliding on a discrete lattice. The total number density of the two-phase fluid is given by where the subscript specifies the direction in the lattice, and the superscript k = 1, 2 denotes fluid 1 and fluid 2 (red or blue fluid).
There are three steps in this method which are (i) streaming, (ii) collision, and (iii) "recolouring". The streaming step is the same as the standard Lattice Boltzmann Method streaming step. Namely, in one time-step, the particle distributions can move by one lattice spacing along the orthogonal axes, or along diagonals. We use the standard LBM notation DnQm for a simulation in D = n dimensions, and with Q = m velocities on the discrete lattice. In the following, we restrict ourselves to 2D and use the D2Q9 Lattice Boltzmann lattice arrangement shown in Fig. 1. In this lattice, we define f k ( , t) as the number density of particles of fluid k moving in the -direction where the Q = 9 velocities are given by This choice means that 0 is the zero velocity vector and therefore represents stationary particles, and for = (1, … , 8) are the velocities in the eight directions shown in Fig. 1 which is defined such that = − +1 for = (1, 3, 5, 7) . The lattice is unitary so the lattice spacing and time step are Δx = Δt = 1 . The streaming step is specified as (1) The collision step for the two-phase LBM involves two terms and can be written as (Latva-Kokko and Rothman 2005) where the superscript denotes the post collision distributions, and (Δf k ) 1 and (Δf k ) 2 are the two collision terms which represent how the particle distributions change during each time step due to collision (Δf k ) 1 while encouraging colour segregation (Δf k ) 2 . The first collision term is nearly the same as the standard collision term of the LBM and is given by where is the relaxation time and f k,eq ( , t) is the equilibrium distribution which is given by where c s = Δx∕( √ 3 is the speed of sound in the lattice. The above equilibrium distribution is the same as the standard equilibrium distribution except for the rest factor C instead of w . The coefficients C are given by Grunau et al. (1993) where k is a parameter that enables the density of the two fluids to be adjusted (Grunau et al. 1993;Reis and Phillips 2007) and is given by The other weights are the same as the standard LBM. Namely, w 0 = 4∕9 , w = 1∕9 for = 1, 2, 3, 4 and w = 1∕36 for = 5, 6, 7, 8 . The macroscopic density of the two fluids are given by the total density of the fluid is given by and the momentum of the fluid is given by
The second collision term is more complex and there are several forms in the literature. Here, we use the term as written in Reis and Phillips (2007) where ( , t) is the colour gradient, is the angle between ( , t) and and A is a parameter that controls the interfacial tension. The colour gradient ( , t) is calculated according to Mora (2021a) which optimizes isotropy of the gradient as where the are the velocities and b are scalar coefficients of the finite difference approximation of the gradient that is accurate to second order, namely W is given by and w = 0.3 is the weight given to diagonal nearest neighbours relative to orthogonal nearest neighbours in the finite difference calculation of the gradient. Choice of the value of w = 0.3 optimizes isotropy of the gradient at small radius of curvature interfaces such as those that occur in flow through a porous medium (Mora et al. 2021a) for the choice of the interfacial thickness parameter = 0.5 which will be described shortly. This choice of w has an order of magnitude better isotropy than the standard choice of w = 1 from the original paper on the RK colour gradient LBM given by Latva and Kokko (2005), and a factor of two superior isotropy for small radius of curvature interfaces such as those occurring in porous flow than the choice of w = 0.25 derived by Leclaire et al. (2011) based on obtaining isotropic errors in the gradient calculation to second order. As such, the choice of w = 0.3 is optimal to obtain the most accurate pore-scale phenomenology which require isotropy of the colour gradient and hence surface tension, to correctly capture behaviour such as viscous fingering, droplet formation and other pore-scale phenomena.
In Eq. (13), the cosine term is given by and B 0 = −4∕27 , B = 2∕27 for = 1, 2, 3, 4 and B = 5∕108 for = 5, 6, 7, 8 . Reiss and Phillips (Reis and Phillips 2007) have shown that the above parameters yield the correct term for interfacial tension in the Navier-Stokes equations. The final step in the Lattice Boltzmann Method for two-phase flow is a so-called recolouring step, which achieves separation of the two fluids. This is achieved as follows (Latva-Kokko and Rothman 2005) and where f * = ∑ k f k * , and A and ∈ (0, 1] are adjustable parameters that affect the interfacial properties. Namely, affects the interfacial thickness, and A and 1 and 2 affect the interfacial tension. That is, parameter A controls the surface tension, but surface tension in the model is also affected by the two viscosities 1 and 2 , and hence 1 and 2 . One must conduct a numerical experiment of a static droplet and apply the Young-Laplace formula for a given set of viscosities 1 and 2 to obtain the exact relationship between A and surface tension at that set of viscosities. This will be further explained later in the paper. In the above recolouring equation, the equilibrium distribution at zero velocity is given by the standard equilibrium distribution, namely The pressure in the flow field is obtained from the equation of state and can be calculated as In the Lattice Boltzmann Method, one achieves no-slip boundary conditions by "bounceback" boundary conditions at the solid interface. Namely, particle number densities bounce back in the direction they came from at fluid-solid interfaces. The RK model for two-phase flow allows any wetting contact angle w to be specified by setting the densities of the two fluids in the solid region through (Latva-Kokko and Rothman 2005) where w1 is the density of fluid 1 in the solid regions, w2 is the density of fluid 2 in the solid regions, and i is the initial density of the majority component = 2 .

Results
In the following numerical experiments, we explore flow of immiscible fluids involving an invading fluid (red fluid = fluid 1) being injected at the left side of a simple 2D model rock matrix saturated with another fluid (blue fluid = fluid 2) over a range of viscosity contrasts M and wetting angles w , where the viscosity contrast M is defined as The goal is to study the flow regimes over the phase space of (M, w ) . In particular, a major goal is to study the effect of the wetting angle on the flow regimes from viscous fingering when the invading fluid has a lower viscosity than the second fluid, to stable displacement when the invading fluid has higher viscosity. We wish to determine whether the wetting angle plays a critical role in determining the flow morphology and the saturation S which is defined as the percentage of red fluid in the pore space at the moment of breakthrough, when the red fluid breaks through the right side of the model. As the percentage of red fluid equals the percentage of blue fluid evacuated from the model, the saturation is equal to the Recovery Factor in Petroleum Engineering. One can consider the case of M < 1 to be analogous to the case of water being injected into a porous rock matrix filled with crude oil which typically has a viscosity several to a hundred times that of water (for light-intermediate crude oil) up to tens to hundreds times that of water for intermediate to heavy crude oil. In the simulations, we ran at viscosity contrasts ranging from M = 0.01 through M = 100 at intervals of Δ log 10 M = 0.5 , so we have M = (0.01, 0.0316, 0.1, 0.316, 1, 3.16, 10, 31.6, 100) , and for wetting angles from nonwetting w = 180 • through to perfectly wetting w = 0 • at intervals of −22.5 • .
The model consists of a square region of size 300 × 300 pixels initialized with nonoverlapping random sized solid circular particles with radii ranging from r = 5Δx through r = 15Δx , with a minimum separation of 4Δx . Figure 2 shows the model rock matrix that is used in the initial suite of simulations.
In our simulations, we set the densities of the two fluids to be identical with r = b = 1 , and we set the numerical parameter controlling interface thickness in Eqs. (18) and (19) to be = 0.5 , which is one of the typically used values in RK colour gradient LBM simulations. The other parameter affecting interface thickness, in the equation for the relaxation time given by Eq. (11), was also set to the value typically used in RK colour gradient LBM simulations of = 0.98.

Wetting Angle Verification
In order to verify the wetting angle implementation of Eq. (20), we ran several experiments in which a semicircular region of red fluid was initialized above a solid wall, and a simulation was run until there was no further change. Hence, the final configuration will be a deformed droplet on the lower wall such that the wetting angle (angle between the red fluid in the droplet and the wall) will be consistent with the specified angle using Eq. (20). In our simulations, we set the density of the red and blue fluids in the solid regions as  20) does not lead to a unique formula for w1 and w2 , and a bad choice that still obeys Eq. (20) can be numerically unstable. Figure 3 shows the result of the wetting angle tests using angles of w = (60 • , 90 • , 135 • , 180 • ) and indicates that the range of wetting angles can be correctly simulated using the RK colour gradient LBM.

Setting the Simulation Parameters
In the following, we wish to investigate the transition from viscous fingering to stable displacement as a function of the viscosity ratio M, and also, to study the effect of the wetting angle w on the morphology of flow and the saturation at and beyond breakthrough. In order to ensure that we are in the viscous fingering regime, the simulations must be performed at a sufficiently high capillary number such that viscous forces dominate over capillary forces. In addition, we wish for the simulations to be at flow rates that are sufficiently low such that we are well below the laminar to turbulent flow transition of Re t ∼ 2300 and also, where inertial effects are negligible where Re denotes the Reynolds number. It is generally known that inertia can be safely ignored for flows with low enough Reynolds numbers such as Re < 1 which is well below the turbulent transition. In addition, the simulations should be tractable (not too computationally expensive) and accurate. These The capillary number is a dimensionless quantity that relates to relative effect of viscous drag forces versus to surface tension or capillary forces and is defined as where r is the dynamic viscosity of the invading fluid, r is the kinematic viscosity of the invading fluid, r is the density of the invading fluid, u in is the injection rate, and is the surface tension. In the RK colour gradient LBM, the value of A in Eq. (13) has been found to relate to the surface tension as where the scalar ∈ [0.55, 0.79] for a wide range of viscosities of the two fluids r and b (Mora et al. 2021a), and where the colour gradient is calculated with Eqs. (14) through (16) using the optimal isotropy weight of w = 0.3 . The exact value of above can be determined through a simulation of a droplet and application of the Young-Laplace formula where ΔP is the difference in pressure inside versus outside a droplet of radius r 0 . When the capillary number is sufficiently high, the viscous forces are dominant. In the following, we will set the interfacial tension such that Ca ≳ 7 which ensures viscous forces are dominant. The other constraint of sufficiently slow flow rates such that inertial effects and turbulence are negligible can be achieved through the Reynolds number which is another dimensionless quantity that relates to the flow pattern. The Reynolds number is defined as where D is the diameter of a tube, is the kinematic viscosity, and u is the flow speed. In the following experiments, we set the injection flow rate u in such that Re = 0.2 for an assumed diameter D = 4 which is the smallest width between grains in the solid matrix. Obviously, the matrix has larger throat widths than D = 4 , but the assumption of D = 4 is sufficient for the order of magnitude calculations done here. Hence, from Eq. (27) 2 , we will have an injection velocity of u in = 0.01 . Note that the reason we do the above calculations is that we wish to minimize computational time, while ensuring accurate results. Use of Eqs. (28) and (29) to set injection velocities and viscosities achieves this goal. Once the viscosities are set from the desired viscosity ratio M using Eq. (29) and the injection velocity is set using Eq. (28) such that the Reynolds number is at a value of Re = 0.2 , we can apply Eq. (24) to determine the surface tension required to meet the criterion Ca ≳ 7 which ensures we are in the regime where viscous forces dominate and viscous fingering will occur. Namely, from Eq. (24), we have where is a parameter that depends on the viscosities that is typically in the range ∈ [0.55, 0.79] Hence, we can calculate the interfacial tension parameter A required by Eq. (13) using Here, we use a value of = 0.55 and Ca = 10 , so for cases where = 0.55 , the interfacial tension parameter A will be set such that the capillary number will be Ca = 10 in the simulation. For cases where approaches the upper limit of = 0.79 , the capillary number of the simulation will be Ca ∼ 7 . Application of Eq. (31) results in A-values ranging from A = 3.6 × 10 −5 through A = 0.0036 respectively for M = 0.01 and M = 100 which are well within the range that can be accurately simulated using the RK colour gradient LBM. Specifically, this viscosity ratio range can be stably simulated and these A-values can be accurately simulated (Mora et al. 2021a). In summary, simulation parameters are set using Eqs. (28), (29) and (31) which ensures simulations are done at sufficiently high capillary number Ca such that viscous forces dominate over surface tension forces ( Ca ∼ 7 → 10 ), and such that inertial effects are negligible ( Re ∼ 0.2 ). We note that these Ca numbers are significantly higher than in microfluidic experiments which are typically done at Ca ≲ 0.1 which means that such experiments may be in the capillary fingering regime at least at some viscosity ratios. Due to our higher Ca numbers than those of experiments, it is entirely possible that our results may not be directly comparable to microfluidic experiments, either due to the possibility that these experiments are in the capillary fingering regime, or that another phenomenon is playing a role such as the possibility of Bretherton films (i.e. when the invading fluid leaves a film of defending fluid behind, adhered to solid walls). Similarly, we recommend caution comparing to field experiments where Ca numbers are significantly lower. In the following numerical experiments using the parameters as explained above, we simulate flow of a red fluid being injected from the left boundary of our model rock matrix shown in Fig. 2 at a rate of u in , and we impose a constant pressure boundary condition at the right of the model. These boundary conditions are achieved using Zou and He velocity and pressure boundary conditions (Zou and He 1997). Periodic boundary conditions are used in the z-direction.

Typical Patterns of Flow and Viscous Fingering at Breakthrough
In an initial suite, we performed a total of 81 simulations using the porous matrix model shown in Fig. 2  . At this viscosity ratio, the saturation at breakthrough is relatively low. At a viscosity ratio of M = 1 , there are no longer the narrow tendrils, leading to a significantly higher saturation, although the front of the red fluid remains irregular so saturation is far from full saturation. Finally, at a viscosity ratio of M = 100 , there is an almost linear front of the red fluid that invades the blue fluid (stable displacement) and the saturation is over 90%. These results are as expected given knowledge of viscous fingering which is known to occur when a low-viscosity fluid invades a higher viscosity fluid at high capillary number.
The cases for the fully wetting fluid ( w = 0 • ) in Fig. 4 show a somewhat different morphology of the viscous fingers at viscosity ratio M = 0.01 than the non-wetting fluid ( w = 180 • ). Specifically, the viscous fingers for the wetting fluid are somewhat broader and more rounded relative to the fingers of the non-wetting fluid which is consistent with the studies of Stokes et al. (1986), Trojer et al. (2015), Zhao et al. (2016) and Primkulov et al. (2019). However, visually, the saturation level seems to be similar for the wetting and non-wetting cases. This will be studied in detail later in this paper. Figure 5 plots the fluid velocity for the same set of viscosity ratios and wetting angles and clearly shows that when M < 1 , the flow is dominantly through narrow fingers of the invading red fluid. As such, it is expected that the flow pattern will tend to continue to flow through the narrow fingers after breakthrough. The evolution of flow and saturation after breakthrough will be explored in a subsequent section. For the lowest viscosity ratio of M = 0.01 in Fig. 4, the non-wetting case has one narrow finger spanning the model, and a second narrow finger spanning about half of the model near the top. In contrast, the fully wetting case has a broad finger spanning the model centred on the same location as the finger for the nonwetting case and is starting to form the second finger at the top of the model. Despite the main finger being broader, the velocity is only high in a narrow channel in the centre of the finger spanning the model (see Fig. 5). It is interesting to note that the morphology of the fingers and velocity flow paths are somewhat different for the different wetting angles. This shows that the relationship between wetting angle and both morphology and the details of the flow pattern is complex.
To gain a more quantitative understanding of the saturation at breakthrough for the different viscosity ratios and wetting angles, Fig. 6 shows the saturation as a function of M and w at breakthrough, where the saturation is computed as S = n r ∕n pore where n r is the number of fluid sites occupied by the red fluid at breakthrough, and n pore is the total number of lattice sites of the fluid filled pore space. This plot clearly shows the effect of the viscosity ratio with low saturation levels of S ∼ 38 % for the lowest viscosity ratio of M = 0.01 , and S ∼ 91 % for the highest viscosity ratio of M = 100 , with an intermediate value of S ∼ 68 % at a viscosity ratio of unity. However, one cannot see on this plot the  To better visualize the variation of the saturation with wetting angle, Fig. 7 shows the difference between the saturation at a given wetting angle relative to the saturation at a wetting angle of w = 180 • , namely, we plot ΔS(M, w ) = S(M, w ) − S(M, 180 • ) . One can see that the saturation has a general tendency to increase as the fluid becomes more wetting. However, the saturation landscape is not simple with the existence of hills, ridges and valleys. For example, at log 10 (M) = −2 , the sweep (saturation at breakthrough) increases with wettability to a maximum at the optimal wetting angle of w = 22.5 • , followed by a decrease in sweep. Furthermore, there is a valley in the saturation difference landscape in the range of log 10 (M) ∈ [−1, −0.5] , and the saturation does not seem to vary significantly with wetting angle. Figure 8 shows profiles of the saturation as a function of wetting angles for the viscosity ratios of log 10 M = −2 and log 10 M = −0.5 , and clearly shows the peak in saturation at w = 22.5 • for the case of M = 0.01 . This increase in saturation with wettability for the low-viscosity case with M = 0.01 up to a certain angle, followed by a decrease in saturation is consistent with the microfluidic experiments by Zhao et al. (2016) involving vertical posts representing a porous medium. In the case of M = 0.01 , the optimal saturation is S = 49.5 % at w = 22.5 • , compared to a saturation of S = 38.7 % for the non-wetting fluid with w = 180 • , a relative increase of 28%. The profile of the saturation as a function of wetting angle for log 10 M = −0.5 which occurs in the valley of the saturation difference phase space is more complex and shows that the saturation is slightly higher for a non-wetting fluid than for a wetting fluid, the opposite to what is expected. Furthermore, there are two peaks in the saturation at w = 22.5 • and w = 90 • , with the highest saturation occurring for the intermediate wetting angle of w = 90 • . However, the saturation range is relatively small between the highest saturation at w = 90 • of 61.6% compared to the lowest saturation of 56.9% at w = 45 • , a relative difference of only 8%.
In summary, for the model shown in Fig. 2, the saturation difference landscape generally has somewhat higher saturation at higher wettability, but in detail, the landscape is complex with a valley in saturation for moderately small viscosity ratios log 10 M ∈ [−1, −0.50] and a hill at w = 22.5 • for the lowest viscosity ratio of log 10 M = −2 . This result suggests a complex relationship between flow morphology, wettability and viscosity ratio. We will explore the effect of variations in the porous model geometry on the saturation difference phase space in a subsequent section.

Evolution of Flow and Saturation After Breakthrough
In this section, we explore how the flow patterns and saturation evolve after breakthrough. Based on the plots of the velocity shown in Fig. 5, for low viscosity ratios of M = 0.01 , most of the fluid flow at breakthrough is occurring only through narrow fingers. As such, one may expect that after breakthrough, flow will continue to be dominant through these narrow fingers and that the flow pattern will cease to evolve significantly. It is for this reason that many studies focus on the saturation at breakthrough. In the following, we investigate whether or not this expectation is fulfilled. Figure 9 shows snapshots of the evolution of a non-wetting invading fluid through time ( w = 180 • ) both before and after breakthrough. Based on this plot, it is clear that the fluid pathways are not static after breakthrough and the red fluid continues to invade more pore space in the model rock matrix after breakthrough. That is, the fingers and flow pathways Fig. 9 Snapshots showing the evolution of a non-wetting invading fluid ( w = 180 • ) through time up to and beyond breakthrough for a viscosity ratio of M = 0.01 . Breakthrough occurred at t = 89333 (Fig. 4 shows the snapshot at breakthrough) continue to evolve and expand. Figure 10 shows corresponding snapshots of the evolution of velocity through time before and after breakthrough and indicates that although the fingers are evolving and expanding, most of the flow initially remains through the initial finger flow paths, although the number of flow paths increases with time. Figure 11 shows snapshots of the evolution of a perfectly wetting invading fluid through time ( w = 0 • ) before and after breakthrough. There are some differences between these snapshots and those for the case of a non-wetting fluid invading (Fig. 9). Mainly, the morphology is somewhat different with wider fingers for the fully wetting case in agreement with experimental observations (Stokes et al. 1986;Trojer et al. 2015). However, the  (Fig. 4 shows the snapshot at breakthrough) general trend is the same as for the non-wetting case with the fingers expanding and the invading fluid eventually occupying most of the pore space. Figure 12 shows corresponding snapshots of the evolution of velocity through time before and after breakthrough for the case of a perfectly wetting fluid ( w = 0 • ). Again, although the fingers are evolving and expanding, most of the flow initially remains through the fingers. As for the case of the non-wetting invading fluid, the finger flow paths increase in number with time and eventually occupy most of the model.
Finally, Fig. 13 shows the saturation as a function of time prior to and after breakthrough for the non-wetting case ( w = 180 • ) and for the perfectly wetting case ( w = 0 • ) with viscosity ratio of M = 0.01 . This plot shows that the saturation continues to increase after breakthrough, although the rate is decreased initially by around a factor of 2 relative to the rate before breakthrough, and the rate continues to decrease with time. The decreased rate of saturation increase post breakthrough indicates that a considerable fraction of the flow continues to go through the narrow fingers that span the model after breakthrough. Nonetheless, after breakthrough, the fingers continue to evolve and widen, and the porous matrix begins to fill up with the invading fluid and eventually occupies most of the pore space ( S ∼ 90%). This continued increase in saturation is a consequence of the continued evolution and growth of fingers that occurs in the numerical model as explained above. Another interesting feature is that the saturation versus time plots for both the wetting and non-wetting case are almost the same for before and after breakthrough for these constant injection rate experiments. Hence, while breakthrough occurs later for the wetting case, and hence, saturation at breakthrough for the wetting case is higher, the saturation and hence oil recovery factor in the simulation, is similar for both the wetting and non-wetting cases at any given time after commencement of injection. This suggests that final EOR production by water flooding may be more related to the time of duration of water flooding rather than the wettability of the injected fluid.
These results have implications to enhanced oil recovery and suggests that an oil field can continue to produce significantly after breakthrough when the water reaches the production well albeit at a reduced rate due to the flow preferentially following the breakthrough finger. Furthermore, based on Fig. 13, the saturation increases for both injection of a non-wetting and wetting fluid, and in both cases, the saturation after 10 6 time steps is similar and exceeds 90% (i.e. ∼ 10 × t breakthrough ). Hence, while the perfectly wetting case has a greater relative saturation at breakthrough than the non-wetting case (45.9% cf. 38.7% ⇒ a 19% relative increase), if one continues to inject beyond breakthrough, the final sweep injecting both a wetting and non-wetting fluid is similar and exceeds 90%. These results suggest that: (1) subject to economic viability due to dropping rates of production, a well could continue to produce up to and beyond 90% of the total oil contained in the reservoir through continued water injection beyond breakthrough, and (2) the time of water flooding may be the primary factor determining total production rather than fluid wettability. These results, if validated with field and experimental tests, could have major significance in terms of extending global oil reserves by around a factor of two or more.

Effect of Variations in the Porous Matrix Geometry
In order to study the effect of the details and statistics of the model on results, we have conducted ten additional suites of 81 runs using models generated with the same algorithm but using a different random number seed. Namely, random sized solid circular particles with radii ranging from r = 5Δx through r = 15Δx , with a minimum separation of 4Δx were dropped into a square region of size 300 × 300 pixels. Figure 14 shows the original model and the ten additional random realizations of the model. Figure 15 shows plots of the saturation difference ΔS(M, w ) for the ten new runs alongside the phase space for the original run where the order of the plots is such that similar phase space landscapes are plotted adjacent to one another. We observe that there are two main groupings with somewhat similar phase spaces Group 1 consists of four cases (the original model and models 1, 2 and 4) with a phase space in which for most viscosity ratios, the saturation has a tendency to increase with wettability in a manner similar to the phase space of the original model. Namely, aside from the valley at log 10 M ∈ [−1, −0.5] , the saturation is higher for more wetting fluids, and at the lowest viscosity ratio of log 10 M = −2 , there is a hill in the range w ∈ [0 • , 22.5 • ] . Group 2 also contains four cases (models 5, 6, 7 and 8), but aside from a ridge at w < 135 • for log 10 M ∼ 0.5 , the saturation does not significantly depend on wettability. In addition to these two major groups, there are three additional cases of which, two cases (models 9 and 10) have a broad plateau for w < 135 • and log 10 M ≥ 0 . Finally, one case (model 3), shows saturation increasing with wettability at all values of M, but with a main plateau for w < 112.5 • and log 10 M ≥ −1.5 . From these plots, one can see that while there is some tendency for saturation to be higher  Fig. 2 and ten different random model realizations that were created using the same algorithm and input parameters as was used to generate the original model for more wetting fluids, there are many cases (models or M ranges) where the saturation is not significantly affected by wettability. We see that in roughly half of the cases (models: original, 1, 2, 3 and 4), there is a strong tendency for saturation to increase with wettability but not for all M's. Specifically, for four of the five cases, there is a valley at log 10 M ∈ [−1, −0.5] . And there is an optimal wetting angle that is low w ∼ 22.5 • but above zero at log 10 M ∈ [−2, −1.5] for three of the five cases (models: original, 2 and 4) . In contrast, the other six of the eleven cases have no tendency for saturation to increase with wettability for small viscosity ratios log 10 M ≲ −0.5 , although at higher viscosity ratios log 10 M ≳ −0.5 , the saturation does increase with wettability for at least some values of viscosity ratio.
These results suggest that the phase-space saturation landscape depends strongly on the porous matrix geometry. For some porous matrix geometries, one will observe what is widely accepted, a tendency for an increase in saturation with wettability, although the optimal wetting angle that maximizes saturation may be greater than perfectly wetting such as w ∼ 22.5 • . But for other porous matrix geometries, this may not occur, or else, there may be a complex phase space saturation landscape with hills, ridges and valleys which means that at certain viscosity ratios, the saturation will be maximized at certain wetting angles which may range from intermediate wettability such as w ∼ 90 • , to high wettability such as w ∈ [0 • , 45 • ]. Figure 16 shows the snapshots at breakthrough for the lowest viscosity ratio of M = 0.01 for the non-wetting case ( w = 180 • ) and the perfectly wetting case ( w = 0 • ). From these plots, one can see that the simple picture of flow morphology variation with wettability for the low viscosity ratio case of M = 0.01 that we saw in our first suite shown in Fig. 4 is not typical. Rather, it is but one picture of many possibilities. Namely, in Fig. 4, we see that for non-wetting fluids, the viscous fingers were narrow, and these became more rounded with wettability which translated into a higher saturation with wettability at M = 0.01 as shown in Fig. 7. In Fig. 16 showing all eleven suites for M = 0.01 and w = (180 • , 0 • ) , we also see that wettability affects flow morphology. Generally, the fingers and flow patterns at w = 0 • are somewhat broader and/or more rounded than at w = 180 • , but the flow patterns are complex and it is not a simple "cartoon" picture of rounding and broadening of fingers with high wettability. Furthermore, the visual appearance of some rounding and/or broadening does not always translate into an increase in saturation. Namely, the change in flow morphology with wettability at M = 0.01 only translated into an increase in saturation in five of the eleven models that we ran.

Discussion on the Influence of Porous Matrix Geometry
In this paper, we were constrained to quite small models sizes for computational reasons. This means that we do not have large enough models for an accurate statistical representation of a porous medium. For this reason, we chose to run suites for eleven different random realizations of a porous matrix to study the variability in the results due to random variations in our small models. In effect, because the models are small, each different case can be considered as a different porous matrix model, with different details and somewhat different statistics. We found that both the flow morphology variation with wettability and the saturation landscape for each different model were different. What is interesting however is that the differences in the saturation landscape were not all entirely different. Rather, at low viscosity ratios, roughly half of the cases had a tendency for saturation to increase with wettability, and the remainder did not. Furthermore, while each saturation landscape was different, there were two main groups with similar features including hills, ridges and valleys. These results clearly demonstrate that the porous matrix structure plays a vital role in controlling saturation. And the optimization of sweep (saturation at breakthrough) is not necessarily a simple matter of trying to maximize the wettability of the injected fluid. Namely, for some porous matrix structures, saturation may be maximized by setting the wettability to a specific value between fully wetting ( w = 0 • ) and partially wetting ( w ∼ 90 • ), where this optimal wettability depends on the viscosity ratio. For other porous media, the wettability may not affect saturation at low viscosity ratios, in which case, changing the wettability of the injected fluid may be futile.
In this research, we were restricted to small model sizes and 2D which means unrealistically high porosity to enable viable permeability. As such, the details of all results can be expected to change if one had the resources to perform large-scale 3D simulations of realistic porous matrix structures. However, some general conclusions will remain such as the possibility of complex saturation phase spaces with hills and valleys, and the sensitivity of the saturation phase space to the porous matrix structure. Ultimately, one can imagine performing large enough 3D simulations of a real medium such that the saturation phase space could be mapped for that hydrocarbon reservoir, which would then enable the optimal wettability for the given reservoir with its viscosity ratio to be predicted.

Effect of Grid Refinement
In the above, the model was quite small due to computational limitations. In order to test whether model resolution plays a role, we have performed simulations for the wetting and non-wetting cases at a viscosity ratio of M = 0.01 and at a resolution of 901 × 901 which is at three times the initial resolution. We then compare the higher resolution results to those obtained at the lower resolution of 301 × 301 . The high-resolution model contained circular grains of the same size at the same location as the low resolution model. Figure 17 shows the fluid flow pattern at breakthrough for the runs using the higher resolution model, and Fig. 18 shows the flow pattern for the runs using the same model at the lower resolution. We see that the finger morphology is somewhat different in the high and low resolution runs but some key features of the flow persist. Namely, there are two main fingers and these are broader for the wetting fluid than for the non-wetting case. This is what we have observed in previous sections and is consistent with microfluidic experiments and leads to the wetting case having higher saturation at breakthrough than the non-wetting case. However, these fingers for the wetting case are broader in the low resolution run than in the high resolution run which means that the saturation increase of the wetting case relative to the non-wetting case for the higher resolution run is not as great as in the low resolution case as shown in Table 1. The main reasons are the decreased finger width coupled with  differences in the fingering pattern. Namely, in the higher resolution run at w = 0 • , the lower finger does not progress as far as in the lower resolution run.

Interpretation and Implications of Grid Refinement
Based on the above, we can see that while the resolution of the run has a significant effect on the details of the finger morphology, the key properties of flow remained the same. Namely, there are two main fingers model, with broader fingers for the wetting fluid and hence, greater saturation at breakthrough for the wetting case. This implies that realistic large-scale 3D simulations will be required to obtain accurate results, but that the general features of flow are likely to remain unchanged as resolution is increased. As such, we expect our main results and conclusions to be qualitatively correct. We suggest that the differences between the two cases mentioned above are primarily related to the relative scale of the fluid interface thickness in the model (a few lattice spacings) relative to the throat size in the porous medium ( 12 + for the high-resolution model and 4 + for the low resolution run). If this is so, then our grid refinement study suggests that care must be taken when applying the RK LBM to study viscous fingering and flow in porous media. For example, is the model large enough to eliminate the finite interface thickness effect in the LBM?

Conclusions
The Lattice Boltzmann Method is a flexible computational tool that allows, among other things, one to simulate multiphase fluid flow in complex heterogeneous media. In this paper, we present 2D simulations of immiscible two-phase fluid flow in a porous medium to study the effect of the wetting angle on the morphology of flow and the evolution of saturation. Eleven suites of 81 simulations totalling 891 simulations were performed of two-phase flow in the regime of viscous fingering (high capillary number) and where inertial effects are negligible (low Reynolds number). Each suite used a different random realization of our 2D porous model. The simulations show, as expected, that when the viscosity ratio is low and a less viscous fluid invades a more viscous fluid, we see viscous fingering behaviour and low saturation levels at breakthrough with narrow fingers of the low viscosity fluid invading the high viscosity fluid. Similarly, as expected, when the viscosity ratio is above unity, the flow tends towards stable displacement of a front of the invading fluid which becomes increasingly linear as the viscosity ratio increases. However, the effect of wetting angle was not as simple. We see a different morphology of the fingers for more wetting fluids, typically with more rounded and somewhat broader fingers. However, there is no simple relationship between wetting angle and saturation. Sometimes, the different finger morphology for wetting fluids translates into higher saturation, but for other cases, this is not so. For example, in nearly half of the suites (5 out of 11) at a low viscosity ratio of log 10 M = −2 , we observe an optimal wetting angle of w ∈ [0 • , 22.5 • ] that maximizes saturation at breakthrough at low viscosity ratios of log 10 M ∈ [−2, −1.5] , with of order 30% higher saturation relative to the case of a non-wetting fluid. However, at moderately low viscosity ratios of log 10 M ∼ −0.5 , four of these five suites exhibit no increase in saturation with wettability. For the remaining half of the suites, there was no increase in saturation with wettability at low viscosity ratios for log 10 M < 0.5 , and for one suite of the eleven, saturation tended to increase with wettability at all viscosity ratios. For low viscosity ratios at breakthrough, most of the flow is through narrow channels centred on the viscous fingers and the saturation is typically of order 35%. After breakthrough, the flow pattern evolves with these fingers expanding, and after ten times longer than the breakthrough time, the saturation exceeds 90%, both for non-wetting and wetting injected fluids. This suggests that while wettability can be optimized to maximize sweep (saturation at breakthrough), at least at certain viscosity ratios and for certain models, wettability matters little to the ultimate productivity of a hydrocarbon reservoir if it was viable to continue water injection beyond breakthrough for a long time.
In summary, our results point to a complex relationship between sweep and wettability dependent on the statistics and details of the porous model and the viscosity ratio. Namely, the phase space of saturation difference as a function of wettability and viscosity ratio may have hills, ridges and valleys. The results suggest that at least for certain porous medium models, one can optimize sweep by a specific choice of wetting angle for the viscosity ratio appropriate to a given reservoir. And furthermore, saturation may continue to increase post-breakthrough albeit at a lower and decreasing rate. The above conclusions require further verification with experiments, especially since our simulations were done at a significantly higher capillary number than microfluidic experiments.
This paper involves small-scale 2D simulations of two-phase flow in a porous medium, and as such, further research using larger scale and 3D models is required to verify how our 2D conclusions translate to realistic 3D examples. These intriguing results require experimental validation and further research to understand. This is particularly so given that our simulations were done at significantly higher capillary numbers than those in microfluidic laboratory experiments and field studies. As such, further research is required to investigate whether the general conclusions in this paper still hold at lower capillary numbers more comparable to those in laboratory and field studies before our results can be applied in practice. Nonetheless, the results demonstrate the potential of the LBM to be used to study multiphase flow and viscous fingering phenomena, and to uncover unexpected behaviour relevant to Enhanced Oil Recovery.