Lattice Boltzmann simulation of immiscible fluid displacement in porous media: Homogeneous versus heterogeneous pore network

Carbon capture and storage (CCS) is a method of reducing anthropogenic emission of greenhouse gases into the atmosphere thereby mitigating global climate change. In CCS, carbon dioxide (CO₂) is captured from power plants or other large point-source emitters, purified, compressed, and injected into subsurface formations for long-term sequestration. Deep saline aquifers are considered as the most ideal candidate reservoirs for sequestering CO₂ because they are geographically widespread, have large potential capacities for storage, and are not used for water supply.¹ When CO₂ is injected into deep saline aquifers, it exists in a supercritical state and displaces the formation fluid from the pore space in a variety of possible saturation patterns, depending on the relative strength of capillary and viscous forces, as well as large and small scale geological heterogeneities.^2–4 Fingering and displacement patterns at the pore scale strongly influence the CO₂ storage process within the reservoir in terms of storage capacity, security, and ultimate fate of the injected CO₂. Therefore, it is of paramount importance to study and understand the mechanisms of immiscible fluid displacement in realistic porous media.

Experimental studies of immiscible fluid displacement focus on using two categories of porous media: natural media, such as rock cores and transparent network models collectively known as micromodels. Rock cores are advantageous for characterizing individual formations, but suffer from difficulty in monitoring fluids at the pore scale since sophisticated and unique micro-tomographic facilities are needed to visualize the internal distribution of the fluids within the rock cores. Moreover, it is challenging to independently manipulate porosity, pore size, connectivity, and wetting properties for natural porous media. These limitations can be overcome by micromodels, which are two-dimensional (2D) pore network patterns etched into materials such as silicon, glass, polyester resin, and most recently, polydimethylsiloxane (PDMS).⁵ Micromodels allow for visualization of fluid distribution using cameras with or without fluorescent microscopy, and subsequent quantification of fluid saturation and interfacial area may provide mechanistic insight about physical displacement process at the microscopic level. For example, Lenormand et al.⁶ performed a series of classic displacement experiments for several fluid pairs in an oil-wet micromodel constructed of a polymer resin, and established a phase diagram delineating parameter domains for stable displacement, capillary, and viscous fingering. They also observed “crossover” behavior in intermediate regions of their phase diagram corresponding to flow morphologies with characteristics of more than one regime. Phase diagram behavior was later demonstrated by Zhang et al.⁷ in a homogeneous water-wet micromodel, and their experimental results showed that the exact locations of the various region boundaries and crossover zones are dependent on the pore network. Micromodels, however, are criticized for the lack of the complex geometry of real media, which often has multiscale and random characteristics that will dictate fluid and solute transport. Numerical simulations can complement experimental studies, providing an economic and efficient pathway to explore the influence of flow and physical parameters in various complicated porous media. However, the numerical methods based on the continuum description are insufficient to consider the influence of pore-scale parameters on the macroscopic bulk properties and, hence, the details of fingering pattern in the porous media cannot be resolved.⁸ Statistical models, i.e., Invasion Percolation (IP), Diffusion-Limited Aggregation (DLA), and anti-DLA, are able to describe certain “specialized” displacement regimes, but they cannot capture transitions from one regime to another.^6,9

The lattice Boltzmann method (LBM) has been developed as an attractive and promising numerical tool for pore-scale simulation of multiphase flows in porous media.^10–14 Unlike pore-network models,^6,15,16 which use a simplified representation of pore geometry and approximate transient flow with a steady-state Poiseuille law, LBM models complex multiphase flows in domain with realistic pore geometries. The fundamental idea of the LBM is to construct simplified kinetic models that incorporate the essential physics of microscopic or mesoscopic processes to ensure that the macroscopic averaged properties obey the desired macroscopic equations. LBM has several advantages over the conventional grid-based computational fluid dynamics (CFD) methods, such as volume-of-fluid (VOF)^17,18 and level-set (LS) methods,^19,20 especially in dealing with complex boundaries, incorporation of microscopic interactions, flexible reproduction of interface between multiple phases, and parallelisation of the algorithm. In the LBM community, a number of multiphase models have been proposed. These models can be classified into four types, i.e., the color-fluid model,^21–23 the interparticle-potential model,^24–26 the phase-field-based model,^27,28 and the mean-field theory model.²⁹ A detailed review of these models for pore-scale flows can be found in Ref. 30.

In this paper, a recently improved color-fluid LB model²³ is used to simulate immiscible displacement of a wetting fluid by injecting a non-wetting fluid at a constant flow rate under a wide range of flow conditions in two microfluidic flow cells, one with a homogeneous pore network and the other with a randomly heterogeneous pore network. We investigate the effect of capillary number, viscosity ratio, and media heterogeneity on displacement patterns. We also quantify the fluid saturations and interfacial areas, and compare the simulation results with the experimental results of Zhang et al.,⁷ who conducted a series of displacement experiments in a homogeneous pore network micromodel with precisely microfabricated pore structures.

The model we are using in the LB simulations is a D2Q9 two-phase lattice Bhatnagar-Gross-Krook (LBGK) algorithm, where we introduce the interfacial tension and phase segregation in accordance with the method of Liu et al.²³ In this model, red and blue distribution functions $f i R$ and $f i B$ are used to represent two different fluids. The total distribution function is defined by $f i = f i R + f i B$⁠, which undergoes a collision step as³¹ 

where $f i ( x → , t )$ is the total distribution function in the i-th velocity direction at the position $x →$ and time $t, f i eq$ is the equilibrium distribution function of $f i , f i †$ is the post-collision distribution function, τ is the dimensionless relaxation time, and Φ_i is the perturbation term. Conservation of total mass and momentum requires

where ρ = ρ_R + ρ_B is the total density with the subscripts “R” and “B” referring to the red and blue fluids, respectively, u is the local fluid velocity, and e_i is the lattice velocity in the i-th direction,³² 

The equilibrium distribution function is calculated by

where w_i is the weight factor with w₀ = 4/9, w₁₋₄ = 1/9, and w₅₋₈ = 1/36. In the above equations, c = δ_x/δ_t is the lattice speed, where δ_x is the lattice spacing and δ_t is the time step.

The perturbation term contributes to the mixed interfacial regions and creates an interfacial tension. The perturbation term reads as^23,31

where A is a parameter directly related to the interface tension, $σ= 4 9 Aτ c 4 δ t$⁠, and ρ^N is the phase field function defined as

The generalized expression for B_i is given by Liu et al.³¹ and it is taken as B₀ = 2c²/9, B₁₋₄ = c²/9, and B₅₋₈ = c²/36 in this work.

The evaluation of ∇ρ^N is required to calculate the perturbation term, i.e., Eq. (6). To minimize the discretization error, ∇ρ^N is evaluated by Ref. 23,

Then, we use the recoloring algorithm proposed by Latva-Kokko and Rothman³³ to produce the phase segregation and guarantee the immiscibility of both fluids. Following their algorithm, the recolored distribution functions of the red and blue fluids are

After the recoloring step, the distribution functions propagate to the neighboring lattice nodes, known as propagation or streaming step,

with the post-propagation distribution functions being used to calculate the densities of both fluids by $ρ k = ∑ i f i k$⁠.

Through the Chapman-Enskog multiscale analysis, Eqs. (1), (9), and (10) can be reduced to the Navier-Stokes equations in the low frequency, long wavelength limit with Eqs. (5) and (6). The resulting equations are²³ 

where $p= 1 3 ρ c 2$ is the pressure, $ν= 1 3 c 2 ( τ − 1 2 ) δ t$ is the kinematic viscosity of the fluid mixture, and F_S is a volume-distributed interfacial force and is given by

Here, $δ Σ = 1 2 |∇ ρ N |$ is the Dirac delta function and I the second-order identity tensor.

To allow for unequal viscosities of the two fluids, we determine the viscosity of the fluid mixture by a harmonic mean

where ν_k (k = R or B) is the kinematic viscosity of fluid k. It has been shown that Eq. (14) can ensure the continuity of viscosity flux across the interface.³¹ 

No-slip boundary conditions at solid walls are implemented by the halfway bounce-back scheme.³⁵ The wettability of the solid walls can be imposed by assuming that the solid wall is a mixture of two fluids, thus having a certain value of the phase field. The perturbation term in Eq. (6) becomes dependent on the properties of the neighboring solid lattice sites, resulting in a special case of the wetting boundary condition. The assigned value of the phase field at sites neighboring the wall sites can be used to modify the static contact angle of the interface.³¹ Similar approaches have been widely adopted by researchers in various multiphase LB models.^36–40 

This two-phase LBGK model (i.e., color-fluid model) has been programmed into a computational software package in our research group and extensively validated through various test cases including verification of Laplace law, droplet deformation and breakup in simple shear flow, single bubble rising under buoyancy force, and the dynamic capillary intrusion.²³ Recently, this model was used to simulate liquid CO₂ displacement of water in a dual-permeability pore network, which reproduces three different displacement patterns observed in the previous micromodel experiments.³¹ More details concerning this two-phase LBGK model can be found in Refs. 23 and 31, which clearly indicate its capability for accurately handling immiscible fluids with variable density and viscosity ratios, as well as fluid flows in complex porous media. Therefore, the two-phase LBGK model is employed for pore-scale simulation of the drainage processes in microfluidic flow cells initially saturated with a wetting fluid, focusing mainly on the interfacial phenomena.

According to the phase diagram proposed by Lenormand et al.,⁶ two-phase displacement in a micromodel, in the absence of gravitational force, can be characterized by two dimensional numbers: the capillary number (Ca) and the viscosity ratio (M). The capillary number describes the relative magnitude of viscous to capillary forces and is defined by Ca = u_nη_n/σcos(θ), where u_n and η_n are the mean velocity and dynamic viscosity of the advancing non-wetting fluid, respectively, and θ is the contact angle between the fluid-fluid interface and the pore wall. The viscosity ratio is defined as the ratio of the advancing non-wetting fluid viscosity to the defending wetting fluid viscosity (η_w), i.e., M = η_n/η_w. In addition, recent laboratory studies^41–44 have shown the strong influence of subcore scale heterogeneities on steady-state migration patterns, spatial distributions, and fluid saturations. To gain a better understanding of pore-scale two-phase displacement mechanisms, a series of numerical simulations are conducted to study the effect of Ca and M on displacement stability and fluid saturation in a homogeneous and a heterogeneous pore networks, and the obtained results are compared to indicate the effect of media heterogeneity.

In our LB simulations, the densities of both fluids are assumed to be unity, the interfacial tension σ = 0.03 and the contact angle θ = 15°. Unless otherwise stated, the dynamic viscosity of the non-wetting fluid is fixed at 0.024. The inlet mean velocity and the dynamic viscosity of wetting fluid are determined by the capillary number and the viscosity ratio, respectively. To match these LB simulation parameters to their physical values, one needs to choose three reference quantities: a length scale L₀, a time scale T₀, and a mass scale M₀. In this study, L₀ = 5 × 10⁻⁶ m, T₀ = 3.6 × 10⁻⁷ s, and M₀ = 1.25 × 10⁻¹³ kg. A simulation parameter with dimensions [m]ⁿ¹[s]ⁿ²[kg]ⁿ³ is multiplied by [L₀]ⁿ¹[T₀]ⁿ²[M₀]ⁿ³ to obtain the physical value. Following this criterion, for example, we can obtain the physical value of density ρ^phy by $ρ phy =ρ M 0 L 0 3 = 1 . 25 × 1 0 − 13 ( 5 × 1 0 − 6 ) 3 =1 0 3 kg/ m 3$⁠, the physical value of interfacial tension σ^phy by $σ phy =σ M 0 T 0 2 =0.03 1 . 25 × 1 0 − 13 ( 3 . 6 × 1 0 − 7 ) 2 =0.0289N/m$⁠, and the physical value of non-wetting phase viscosity $η n phy$ by $η n phy = η n M 0 L 0 T 0 =0.024 1 . 25 × 1 0 − 13 ( 5 × 1 0 − 6 ) ⋅ ( 3 . 6 × 1 0 − 7 ) =1.67×1 0 − 3 Pas$⁠.

First, we investigate the invasion process of a non-wetting fluid in a two-dimensional homogeneous porous media geometry, as shown in Fig. 1(a). The porous media geometry consists of one inlet section and one outlet section, connected by a pore network, which contains a staggered periodic array of uniform circular grains, 80 μm (16 lattices) in diameter, with 97 μm pore bodies, 45 μm pore throats, and a porosity of 0.68. In order to produce fingering, small perturbations are introduced in the simulations. Specifically, the spatial location of each grain center is generated by the perfectly arranged coordinate plus a random perturbation (which obeys a uniform distribution) between −5 and 5 μm in both horizontal and vertical directions. Initially, the pore network is fully saturated with the wetting fluid. The non-wetting fluid is injected continuously at a constant flow rate at the left inlet, while a constant pressure is imposed at the right outlet. The lateral boundaries are assumed to be solid walls. After careful grid independence studies, the computational domain is chosen as 1920 × 1200 lu² (lu: lattice unit) with the pore network 1800 lattices in the horizontal direction. It is worth noting that the present pore network can be comparable to the one used by Zhang et al. in their micromodel experiments⁷ except that we have neglected the effect of depth and simulated the pore network with smaller grains and higher porosity in order to minimize the computing cost. In addition, it is expected that the choice of equal densities for both fluids has negligible influence on the results, since the displacement experiments are typically operated in a horizontally placed micromodel with small depth, where the gravity effect is largely suppressed.

A series of simulations are performed for capillary numbers ranging from logCa = − 5 to logCa = − 3 at two different viscosity ratios, i.e., M = 1/12.5 and M = 1. For each of the cases under consideration, the simulation is run until saturation of the non-wetting fluid reaches the steady state. Fig. 2 shows the final fluid distributions in the entire homogeneous flow cell at various capillary numbers for (a) M = 1/12.5 and (b) M = 1 when the non-wetting fluid is injected from the left inlet. Fingering is clearly observed during the displacement of non-wetting fluid, which is attributed to the small perturbation in position of the grains, whereas in an experiment, it may be attributed to non-uniformity in grain size, grain position, flow channel depth (3D effects), or wettability variability, all associated with the microfabrication technology. Three representative displacement patterns, namely, capillary fingering, viscous fingering, and stable displacement, are identified strongly depending on the values of capillary number and viscosity ratio. At low capillary numbers, i.e., logCa = − 5, the pore body is almost occupied completely by the advancing fluid before it can reach a neighboring pore, consistent with the assumption introduced by Lenormand et al.⁶ in pore-network simulations. The advancing fluid flow perpendicular to the main flow direction is evident and fingers sometimes progress into new pore bodies in the backward direction, indicating capillary fingering. In capillary fingering, the invasion of non-wetting fluid is dominated by the capillary force, which causes the non-wetting fluid to progress preferentially from a pore through the largest pore throat in any direction due to a lower entry pressure. For M = 1/12.5 at high capillary numbers (logCa > − 4), multiple loosely connected or disconnected fingers mainly progress forward toward the outlet boundary with limited or no lateral flows, and are referred to as viscous fingers. In viscous fingering, the capillary force and pressure drop over the invading fluid is negligible, and the dominant force is due to the viscosity of the wetting fluid. The dominant viscous force causes a significant fraction of fingers to occupy only partial pore bodies they pass through. At moderate capillary numbers, the unstable displacement by the less viscous non-wetting fluid shows the features of both capillary and viscous fingers. The zone where at least two types of fingers are observed is typically named as the crossover zone. For the two-phase displacement with M = 1, an increase in capillary number leads to the development of new flow paths that are connected to the initial finger and, hence, the interface fronts become increasingly flat, revealing some features associated with stable displacement. However, the displacement cannot fall into the regime of stable displacement even at the highest capillary number, i.e., logCa = − 3, since the viscosity induced pressure drops in both invading and defending fluids are equally important. As anticipated, increasing the viscosity ratio from M = 1/12.5 to M = 1 can enhance stability during displacement, so for M = 1, the fingers keep occupying completely the pore bodies that they flow through and the trapped blobs of wetting fluid are smaller in size. Finally, we can notice that the trapped blobs of wetting fluid occupy considerably more than five pore bodies at low capillary numbers; whereas, their sizes decrease with an increase in capillary number for a fixed M.

Based on the displacement patterns shown in Fig. 3, we position these simulations (represented by hollow circles) in the logM − logCa stability phase diagram, where we also give the boundaries of stable displacement, and capillary and viscous fingering regions established by Lenormand et al.⁶ (bounded by dashed lines) and by Zhang et al.⁷ (bounded by solid lines). As Ca increases, the immiscible displacement changes from capillary to viscous fingering for M = 1/12.5 (logM = − 1.1), but it changes from capillary fingering to near stable displacement for M = 1 (logM = 0). This trend is consistent with the phase diagrams by Lenormand et al. and Zhang et al., which were developed on the basis of a large number of drainage experiments for several fluid pairs in an oil-wet micromodel and in a water-wet micromodel, respectively. However, region boundaries obtained in the present pore network simulations differ from previous experimental observations. For example, our simulation with the lowest capillary number at logM = − 1.1 reveals considerable capillary fingering and, hence, the upper boundary of capillary fingering region should be located between logCa = − 5 and logCa = − 4.3, suggesting that the capillary fingering covers a region significantly different from those obtained by Lenormand et al.⁶ and Zhang et al.⁷ The differences are mainly attributed to the variations in the geometry of porous media, including porosity, pore-throat size distribution, and the size of the pore network, which will be further discussed in Sec. III B. Since we use an almost homogeneous, isotropic pore network, similar to the one used by Zhang et al.,⁷ it is not surprising that our obtained region boundaries are closer to Zhang et al.⁷ than to Lenormand et al.⁶ 

In order to quantify the displacement efficiency in the homogeneous pore network, Fig. 4 plots the steady state saturation of the non-wetting fluid (S_nw) as a function of logCa for M = 1/12.5 and M = 1, which are represented by the solid and hollow circles, respectively. For each fixed M, the non-wetting fluid saturation increases with an increase in logCa, and the rate of increase is roughly a constant. The same approximately linear increase in non-wetting fluid saturation with logCa was also reported by Cottin et al.⁴⁵ and Zhang et al.⁷ in drainage experiments, where the immiscible two-phase displacement occurs in the form of capillary and viscous fingering. In addition, it can be seen that an increase in M results in higher non-wetting fluid saturation for a fixed Ca. The increased S_nw is attributed to the enhanced stability of two-phase displacement, which is in accord with our numerical observations (see Fig. 2) and the phase diagram.

In addition to the fluid saturations, we also quantify the interfacial areas between the immiscible fluids. The interfacial area is an important factor influencing mass and energy transfer among phases, e.g., it can strongly influence CO₂ dissolution in the formation brine and subsequent geochemical reactions with brine and host rock in CO₂ sequestration. Following the work of Zhang et al.,^7,43 we use the interfacial length instead of the interfacial area in the present 2D simulations. This length includes not only the length of interfaces between wetting and non-wetting fluids in the pore bodies and throats but also the length of the wetting films between the non-wetting fluid and the grain surfaces. Fig. 5(a) shows the variation of specific interfacial length with non-wetting fluid saturation during the displacement process, where the specific interfacial length is defined as the ratio of the interfacial length to the pore area. For each case, the specific interfacial length varies roughly linearly with S_nw in the entire process of displacement. All data are collapsed onto a single line except two highest Ca for M = 1/12.5, which exhibit a higher slope in the relationship between the specific interfacial length and S_nw. The higher slope is a result of viscous fingering in which some fingers are stretched very thin and even broken into small blobs, thus having a higher interfacial length than in capillary fingering and stable displacement. Similar to Refs. 7 and 43, we also plot the steady state specific interfacial length (l_s) as a function of S_nw, which is shown in Fig. 5(b). It is observed that l_s is linearly proportional to S_nw for each M, and the slope for low M (M < 1) is slightly higher than that for high M (M ≥ 1), which is in good agreement with previous experimental findings of Zhang et al.⁷ 

Having verified that our LB simulations can reproduce the experimental observations in a homogeneous pore network micromodel, we investigate the immiscible two-phase displacement in a heterogeneous pore network. A schematic illustration of the computational domain and boundary conditions is shown in Fig. 1(b). The heterogeneous pore network is composed of randomly distributed circular grains with the grain diameter (d) determined by

where x is a uniformly distributed random number within the interval [0, 1], d_max and d_min are the maximum and minimum limited values of grain diameter, and δ is a parameter between 0 and 1, which is selected to obtain the prescribed standard deviation. In our simulations, the standard deviation is 7.5 μm, and d_min and d_max are taken as 60 μm and 100 μm, respectively, which gives a mean grain diameter of 80 μm. To ensure sufficient computing accuracy and grid independent solutions, we assign at least 5 lattices in the gap between two neighboring grains.³¹ The pore network has the same size and porosity as the above homogeneous pore network, which allows for a direct comparison between the two-phase displacements.

Figure 6 shows the final fluid distributions in the heterogeneous microfluidic system at various Ca for M = 1/6 and M = 1. As Ca increases, for M = 1/6, we can clearly observe the transition from capillary to viscous fingering, which occurs at around logCa = − 3.3; whereas for M = 1, the displacement pattern changes from capillary fingering to near stable displacement with fingers revealing some typical features associated with stable displacement when logCa ≥ − 3.3, such as the invaded pore bodies fully occupied by the non-wetting fluid and relatively flat advancement of the interface fronts. These observations are also reflected in logM − logCa phase diagram (see hollow triangles in Fig. 3). In viscous fingering, i.e., logCa ≥ − 3 and M = 1/6, the dominant viscous force from the wetting fluid causes some fingers to break up, so we can see some isolated blobs of non-wetting fluid trapped in the pores. Also, these blobs decrease in size as Ca increases. In contrast, we have not seen the breakup of fingers at logCa = − 3 and M = 1/12.5 in the homogeneous pore network (see Fig. 2) even though a lower viscosity ratio can lead to an increasingly instability. This difference is attributed to the media heterogeneity, which can promote instability of two-phase displacement. As shown in Fig. 6(b), one can observe considerable capillary fingering (evident lateral flows and back loops) for the capillary number as high as logCa = − 4. Under the same flow conditions (logCa = − 4 and M = 1) in the homogeneous pore network, however, the displacement is located in the crossover zone where both capillary fingering and stable displacement are observed. The different regimes also suggest that the unstable displacement is enhanced by the media heterogeneity. Finally, the residual blobs and pools of wetting fluid decrease in size with an increase in Ca for each constant M in the heterogeneous pore network, consistent with our findings in the homogeneous pore network. But it is found that the heterogeneous pore network usually leads to larger blobs and pools of wetting fluid left behind than its homogeneous counterpart.

Figure 4 also plots the steady state saturation of non-wetting fluid as a function of logCa in the heterogeneous pore network, which is compared with the results in the previous homogeneous pore network. As observed in the homogeneous pore network, S_nw exhibits an approximately linear dependence on logCa in the heterogeneous pore network for each M and also, S_nw is higher for a larger value of M at a fixed Ca. However, it is noticed that S_nw is lower in the heterogeneous pore network than in the homogeneous pore network, which results from the increased instability caused by the media heterogeneity. This can reasonably explain the conventional wisdom in reservoir engineering that heterogeneity would reduce sweep efficiency.

Fig. 7(a) presents the specific interface length as a function of S_nw during the invasion of a non-wetting fluid in the heterogeneous pore network for all cases investigated. Consistent with our previous observations, all data points are collapsed onto a single straight line through the origin except those indicating considerable viscous fingering, i.e., logCa ≥ − 3 and M = 1/6. It is believed that the results will deviate from the collapsed line more significantly with a further increase in Ca or a decrease in M, since both can help the growth of viscous fingering as per the logM − logCa phase diagram. In addition, the steady state specific interfacial length l_s exhibits a linear dependence on S_nw for either unfavorable (M < 1) or favorable (M ≥ 1) displacement, i.e., l_s = kS_nw, where the slope k is a fitting parameter. Based on least-square fitting, the resulting equations are l_s = 331.25S_nw and l_s = 320.83S_nw for M = 1/6 and M = 1, respectively. Obviously, the slope for the unfavorable displacement is slightly higher than that for the favorable displacement, consistent with our numerical observation in the homogeneous pore network. From the fitting equations, it is obtained that the maximum values of l_s are 331.25 cm⁻¹ and 320.8 cm⁻¹ for M = 1/6 and M = 1, respectively. The maximum values of l_s are both close to the specific solid surface length 316.33 cm⁻¹ (calculated as the solid surface length divided by the pore area), which occurs when the non-wetting fluid occupies all of the pore space, i.e., S_nw = 1. The linear relationship between specific interfacial area and non-wetting fluid saturation has also been shown experimentally in one-dimensional column studies during drainage.^46–48 

As we know, it is still challenging to simulate a wide range of viscosity ratios for many existing multiphase LB models.³⁰ To examine the range of viscosity ratios that the present model can access, we simulate the displacement of a more viscous fluid by a less viscous one in the heterogeneous pore network at logCa = − 5. The dynamic viscosities of both fluids are chosen as η_n = {0.005, 0.01, 0.02, 0.04} and η_w = {0.1, 0.3, 0.5} in which η_w is limited to not more than 0.5 because as η (or τ) increases, so does the Knudsen number (defined as the ratio of the mean free path to the characteristic length scale). It is found that regardless the value of η_w, the model remains stable for ν_n ≥ 0.02, but becomes unstable for ν_n < 0.02. Thus, the lowest viscosity ratio that the present two-phase LBGK model can access is $M= 0 . 02 0 . 5 =1/25$⁠. Fig. 8 shows the final fluid distributions in the heterogeneous microfluidic system at logCa = − 5 and M = 1/25. It is clearly seen that the two-phase displacement exhibits evident capillary fingering. As expected, the steady state saturation of non-wetting fluid S_nw = 0.235 is much lower than that for logCa = − 5 and M = 1. In addition, the specific interfacial length l_s is linearly dependent on S_nw during the displacement (see Fig. 9), consistent with the previous observation in Fig. 7(a). To simulate lower viscosity ratios, one of feasible approaches is to replace the present LBGK model with a multiple-relaxation-time (MRT) model.⁴⁹ It has been recently demonstrated that the color-fluid MRT model can simulate the immiscible two-phase displacement with the viscosity ratios ranging from 1/200 to 500.⁵⁰ 

In this work, our recently developed color-fluid LB model²³ is used for pore-scale simulation of the drainage process in two microfluidic flow cells, one with a homogeneous pore network and the other with a randomly heterogeneous pore network, which have the same size, porosity, and mean grain diameter. The extent and behavior of preferential flow (i.e., fingering) is found to depend on the capillary number (Ca), viscosity ratio (M) and media heterogeneity, and three different displacement patterns observed in the previous micromodel experiments are reproduced: for M ≥ 1 at low Ca, displacement occurs in the form of capillary fingering due to the dominant capillary force, and the displacement gradually becomes stable with increasing Ca; for M < 1 at low Ca, capillary fingering is again exhibited, but viscous fingering dominates at higher Ca because of the lower viscosity of the displacing non-wetting fluid. In the logM − logCa stability phase diagram, our simulation results show that the boundaries of displacement regimes (i.e., predominantly stable displacement, and capillary and viscous fingering) are different for the homogeneous and heterogeneous pore networks, and their boundaries of regimes also differ from those obtained by Zhang et al.⁷ and Lenormand et al.⁶ due to difference in the configuration of the micromodel pore network. This suggests that the location of each regime boundary needs to be estimated by performing experimental or numerical studies on each specific pore network. For a fixed M, the non-wetting fluid saturation (S_nw) increases roughly linearly with logCa, consistent with the experimental observations by Cottin et al.⁴⁵ and Zhang et al.⁷ S_nw is higher for a larger value of M since an increase in M can enhance the stability of the displacement process. On the contrary, media heterogeneity weakens the displacement stability, resulting in lower S_nw in the heterogeneous pore network than in the homogeneous one. During the invasion process of the non-wetting fluid in either pore networks, the specific interfacial length is linearly proportional to S_nw with an identical proportionality constant for all cases except those primarily showing viscous fingering. This is because some viscous fingers can only occupy partial pore bodies that they pass through and even break up to form small blobs, leading to a larger interfacial length. As observed in the micromodel experiments of Zhang et al.,⁷ the steady state specific interfacial length and S_nw can be well correlated by a linear relationship for either favorable (M ≥ 1) or unfavorable (M < 1) displacement, and the slope for the unfavorable displacement is slightly higher than that for the favorable displacement.

A recent pore-scale experimental study⁵¹ quantified the strong velocity variations in single- and multiphase flows within a three-dimensional porous media. Although the pore space is highly disordered and complex, it was found that the velocity magnitudes and the velocity components both along and transverse to the imposed flow direction are exponentially distributed. In future, we will develop a three-dimensional high-performance computing code to examine whether our LB simulations can reasonably capture the experimental findings.

The authors gratefully acknowledge the support of the International Institute for Carbon Neutral Energy Research (WPI-I2CNER), sponsored by the Japanese Ministry of Education, Culture, Sports, Science and Technology. This work is also financially supported by the Engineering and Physical Sciences Research Council (EPSRC) of the UK under Grant Nos. EP/I036117/1 and EP/I011927/1. Y.H.Z. would also like to thank the UK’s Royal Academy of Engineering (RAE) and the Leverhulme Trust for the award of a RAE/Leverhulme Senior Research Fellowship.