Dynamic stability characteristics of fluid flow in CO₂ miscible displacements in porous media

Abstract

The dynamic characteristics of fluid flow are important in miscible displacement processes in carbon dioxide enhanced oil recovery (CO₂-EOR) projects. And the stability of the in situ mixing zone greatly influences the oil recovery factor, which deserves further research. We investigated CO₂ miscible displacement processes using magnetic resonance imaging (MRI) apparatus. The CO₂ miscible displacement flows were performed at a low injection rate of 0.1 ml min⁻¹ with reservoir conditions of 8.5 to 9.5 MPa and 37.8 °C. The oil saturation evolution, the length of the in situ mixing zone, and the mixing-frontal velocity and CO₂-frontal velocity were quantified. The experimental results showed that the residual oil saturation decreased with pressure and the mixing zone length was independent of pressure. The mixing-frontal velocity and the CO₂-frontal velocity were nearly the same and increased with pressure. The critical velocity of the CO₂/n-decane (CO₂/nC₁₀) system was 1.105 × 10⁻⁵ m s⁻¹. Although the whole mixing zone length had no obvious change with pressure, a higher pressure compressed the mixing zone and led to an unstable mixing front above the critical velocity. The longitudinal dispersion coefficient was calculated by fitting the experimental data with an error function, which had no obvious change with pressure. Additionally, a three-dimensional lattice-Boltzmann method (LBM) was used to simulate pore-scale miscible fluid flows in upward vertical displacements. A front fingering occurred at a low kinematic viscosity ratio (ν_Co2 : ν_o = 1 : 1). At a large kinematic viscosity ratio (ν_Co2 : ν_o = 1 : 15), the high kinematic viscous oil restrained the buoyancy of supercritical CO₂, but also impeded the displacement with a pore-scale backflow which might lead to a low oil recovery factor.

---

1 Introduction

In a light or medium oil reservoir, approximately 50 to 60% of the original oil in place (OOIP) remains after primary and secondary oil recoveries.^1,2 Enhanced oil recovery (EOR) projects are seriously considered by oil producers. CO₂ miscible displacements with high recovery factors have been recognized as one of the most effective EOR technologies during secondary and tertiary oil recoveries.^3,4 For example, during secondary oil recovery, CO₂ displacement significantly reduces residual oil by as much as 60% OOIP, which is more than water displacement (44% OOIP).⁵ During tertiary oil recovery, miscible CO₂ displacement can reduce the residual oil by 8 to 25% OOIP in a field-scale pilot test.^6–8 Additionally, CO₂ EOR is an environmentally friendly technology that can reduce greenhouse gas effects.^9,10

Understanding the dynamic characteristics of fluid flow in CO₂ miscible displacement is important for predicting the in situ mixing position and oil production in oilfields. For example, in Chevron's McElroy Field¹¹ in West Texas, USA, a variety of porosities, permeability and heterogeneity cause the field to have various flooded zones and oil recovery efficiencies.¹² The effects of CO₂ miscible displacements have been previously discussed in various numerical and experimental studies, such as flow rate, heterogeneity, diffusion, viscosity and density differences, grain size distribution and grain shapes.¹³

Experimental studies have been performed to investigate the flow characteristics of miscible displacement processes. Al-Wahaibi et al.¹⁴ studied the effects of the flow rate, the gravity effect, the bead size and the length scale on gas–oil non-equilibrium oil recoveries for multi-contact miscible (MCM) displacement at normal pressures and temperatures. The miscible non-equilibrium increased with flow rate but independent of the permeability and the length of the bead-pack. Torabi and Asghari¹⁵ investigated the effects of connate water saturation, the matrix permeability and the oil viscosity on the performance of gravity drainage from the matrix into the fractures. The ultimate oil recovery was less sensitive to the matrix permeability at pressures near or above the minimum miscibility pressure (MMP). Rashid et al.¹⁶ focused on the factors affecting oil and gas recovery in CO₂-EOR processes. Core-flooding experiments showed that reservoir permeability differences of up to 1 order of magnitude did not affect the CO₂-EOR factor. Wylie and Mohanty^17,18 investigated the effects of water saturation and wettability on the bypassing and mass transfer under near miscible conditions of gas displacements. Torabi and Asghari¹⁹ measured the effects of the operating pressure, matrix permeability and oil in place and the connate water saturation on oil recoveries for CO₂ huff-and-puff processes. They observed a drastic increase in the recovery factor from immiscible to the near-miscible/miscible pressures²⁰ and found that the recovery may decrease far above the miscibility. Trivedi and Babadagli²⁰ examined the effects of flow rate and reservoir back-pressure on the critical rate, oil recoveries and supercritical CO₂ storage during the first-contact miscible (FCM) displacement in artificially fractured cores.

All of these experimental studies used traditional methods identical to black-box experiments.²¹ The oil saturation evolution in miscible displacement process is measured using the fluid qualities at the inlet and outlet in traditional methods. However, mixing zone performances cannot be visualized or measured using traditional methods. Magnetic resonance imaging (MRI)^22–25 and computed tomography (CT)^26–28 can detect and observe CO₂ miscible displacement processes at reservoir conditions with elevated pressures and temperatures. These methods are recognized as effective approaches to solve visual restrictions at reservoir conditions in porous media.^29–31 Using MRI technology, Liu et al.³² performed a displacement of a CO₂-n-decane (CO₂–nC₁₀) system at miscible conditions of 8.5 MPa and 40 °C above the MMP. (The minimum miscible conditions for CO₂/nC₁₀ system are determined to be 7.79 MPa at 37.8 °C.³³) The authors found the stable front in the homogeneous bead-pack core. Zhao et al.^29,30 proposed a method to calculate the average CO₂ front velocity of CO₂ miscible displacements in porous media using MRI. Berg et al.²⁷ investigated the miscible displacement performance of a CS₂–nC₁₀ system at 1.5 MPa and 25 °C using CT scanning. The authors concluded that the front fingering could affect the mixing zone area. Additionally, Song et al.³⁴ investigated the in situ mixing zone performance of CO₂ miscible displacement flows using MRI. The in situ mixing zone parameters, such as mixing-frontal and CO₂-frontal velocities, the mixing zone length, the longitudinal dispersion coefficient and the Peclet number were quantified. They found a volumetric contraction phenomenon during CO₂ miscible displacement processes.

With regard to numerical methods, conventional level set-based approaches have an inherent difficulty in tracking miscible fluids due to its discrete treatment for interface.³⁵ Hence, it is desired to develop a new numerical scheme. Since the early 1990s, the lattice-Boltzmann method (LBM) has been developed as an alternative and promising numerical scheme for simulating flows of viscous, multicomponent, and multiphase fluids.³⁶ The LBM is especially useful for modeling complicated boundary conditions and multiphase and multicomponent interfaces.³⁷ Gunstensen et al.³⁸ introduced a two-color LBM for simulating immiscible flows of binary fluids in microscopic porous media. Holme et al.³⁹ presented another two-color LBM for simulating binary miscible flows and found miscible viscous fingerings in 2-dimension (2D) simulations. Flekkøy⁴⁰ developed a new two-color LBM for two-component miscible fluid flows in both 2D and 3D. Shan and Doolen^41,42 proposed a multiphase and multicomponent LBM model including interparticle interaction and external forces.⁴³ Most of these researches focused on multiphase and multicomponent fluid flows in bulk volume. However, there are less researches on the behaviors of fluid flows in miscible displacements in microscopic porous media, which couple the fluid-solid effects together. In this work, miscible fluid flows in micro-pores were studied by using a multiphase and multicomponent LBM model.

The experiment and simulation methods were both used for investigating the dynamic stability characteristics of fluid flows in CO₂ miscible displacement in porous media in this study. The experiments were performed at a constantly low injection rate of 0.1 ml min⁻¹ and reservoir conditions of 8.5 to 9.5 MPa and 37.8 °C by using an MRI apparatus. Some miscible flow characteristics, such as the oil saturation evolution, the in situ mixing zone length, the mixing-frontal and CO₂-frontal velocities, the longitudinal dispersion coefficient were measured and analyzed. Additionally, a three-dimensional LBM was used for simulating miscible displacement flows in micro-pores.

2 Material and methods

2.1 Experimental approach

2.1.1 Experimental materials. nC₁₀ (99%) was used as the displaced oil and CO₂ (99.99%) was used as the displacing fluid. The core holder was tightly packed with soda glass beads of BZ-02 made by AS ONE, Japan. The bead size distribution was 0.177 to 0.250 mm with an average diameter of 0.20 mm. The effective porosity of the core was 34.8%. The absolute permeability was 13.8 D. The critical point of CO₂ was 31.0 °C and 7.38 MPa.⁴⁴ The MMP of the CO₂–nC₁₀ system was 7.79 MPa at 37.8 °C.³³ To ensure the miscibility of the CO₂–nC₁₀ system, the experimental conditions were chosen to be 8.5 to 9.5 MPa and 37.8 °C. The thermophysical properties of CO₂ and nC₁₀ at 37.8 °C and 8.5 MPa is shown in Table 1.

Table 1 Thermophysical properties of CO₂ and nC₁₀ at 37.8 °C and 8.5 MPa

Substance	Density (kg m^−3)	Volume (10^−3 m^3 kg^−1)	Viscosity (10^−5 Pa s)	Phase
CO2	447.32	2.2355	3.1895	Supercritical
nC10	723.60	1.3820	78.383	Liquid

---

2.1.2 Experimental apparatus. Fig. 1 shows the schematic of the MRI experiments for CO₂ miscible displacement. The apparatus consisted of an MRI system and a displacement system.

Fig. 1 Schematic of the MRI experiments of CO₂ miscible displacement.

The MRI system was used for image acquisition. A 400 MHz, 9.4 T NMR system (Varian, Inc., Alto, CA, USA) with a maximum gradient strength of 50 Gauss cm⁻¹ was used for the displacement experiments. The inner diameter of the NMR wide-bore was 89 mm. A standard spin-echo multislice scan (SEMS) was used to acquire images during the displacement. The echo time (TE) was 1.37 ms, the repetition time (TR) was 2000 ms, and the acquisition time was 3.2 min. The selected field of view (FOV) was 40 × 40 mm, and the image data matrix was 96 × 96.

The displacement system was used for CO₂ miscible displacement. Supercritical CO₂ was injected using a syringe pump with temperature-controlled circulation to maintain the supercritical CO₂ conditions. The oil was injected using another syringe pump. The syringe pumps (Model 260 D, Teledyne Isco. Inc., Louisville, KY, USA) set a flow rate, the back-pressure regulator (Model BP-2080-M, JASCO, Tokyo, Japan) maintained the core pressure at a reservoir high-pressure condition, and a low differential-pressure transmitter was used to measure the pressure drop through the bead-pack core during the displacement process. A special polyimide bead-pack holder with a maximum limit of 15 MPa and 70 °C was used during the displacement process. The bead-pack holder had an inner diameter of 15 mm and was 200 mm long. The holder had a temperature-controlled circulator (Fluorinert FC-40 with no MRI signal⁴⁵ was chosen as the circulating fluid) to maintain a constant reservoir temperature during the displacement process. The temperature control range of the circulator was −45 to 200 °C with a precision of ±0.5 °C. The fluorinert flowed in the outer tube of the holder to control the temperature of the holder. The thermocouples were placed at the inlet and the outlet of the core to log the temperature change.

2.1.3 Experimental procedure. The experimental procedure is as follows: first, a tube was vacuumed for at least 30 minutes using the vacuum pump. The bead-pack core was then saturated with oil at a constantly slow rate of 0.1 ml min⁻¹. The temperature-controlled circulation was maintained at 37.8 °C throughout the displacement process. Supercritical CO₂ was then injected vertically upward into the bead-pack core at a certain injection rate. Sequences of pressures from 8.5 to 9.5 MPa were performed in the BZ-02 core at 0.1 ml min⁻¹ (2.71 × 10⁻⁵ m s⁻¹) and 37.8 °C. The pressure was adjusted using the back-pressure regulator. To minimize the sand packing effects on pore structures, the tests were performed on the same core without dismounting the core holder. The MRI scan started before the displacement and stopped when the MRI signal intensities had no change.

2.1.4 MRI slice selection. A 16 mm and 2 mm slice selection of the MRI scan was used in this study and is shown in Fig. 2. CO₂ and the rock matrix had no signal intensity in the MRI scan and appear as dark regions in the images. The oil that contained ¹H had MRI signal intensities and appears as bright regions in the images as shown in Fig. 2. The MRI scan order was three 16 mm slices alternated with a 2 mm slice. The 16 mm MRI slice was used to reconstruct the pseudo-3D images shown in Fig. 2a. The thick 16 mm slice, which contained 40 × 40 × 16 mm MRI signals, reflected the oil saturation evolution in the FOV. The pseudo-3D images were displayed in a plane, even though it contained 3D information. Therefore, the oil saturation transversely decreased from the center to the edge in the core holder. At the edge of the core holder, the oil saturation was notably lower and was blue, as shown in the 16 mm MRI image in Fig. 2a. The thin 2 mm slice contained 40 × 40 × 2 mm MRI signals, as shown in Fig. 2b. The oil saturation of the 2 mm slice was uniform in transverse, and no blue edge was found at the edge of the core holder. Therefore, the thin 2 mm MRI slice was used to measure the frontal shape of the miscible displacement. The 16 mm MRI slice was used to measure the oil saturation evolution in the FOV.

Fig. 2 Slice selection for MRI image reconstruction.

2.2. Mathematical model of multicomponent lattice-Boltzmann method

Lattice-Boltzmann model, which is a relatively new method for mesoscopic kinetic simulations,⁴⁶ is suitable for modeling processes of fluid flow in micro-channels and micro-pores.⁴⁶ In order to achieve the multiphase and multicomponent fluid behaviors, an attractive force⁴⁷ between nearest neighboring particles was proposed by Shan and Chen.⁴⁸ This LBM is based on the rule of the single relaxation time approximation⁴⁹ proposed by Qian et al.⁵⁰ and Chen et al.,⁵¹ which is called lattice Bhatnagar–Gross–Krook (BGK) model.⁵² In this work, a D3Q19 multicomponent LBM developed by Shan and Chen was adopted for simulating miscible displacement flows in micro-pores.^41,46,47,53

As shown in Fig. 3, there are 19 discrete internal velocities e_i (i = 0, 1, …, 18) in each grid of D3Q19 model. And each velocity has three motor directions x, y and z. The basic rule of the LBM is:

where f_i^k(x,t) is the density distribution function of σ component moving in lattice direction e_i at node x and time t; e_i is the discrete internal velocity, i = 0, 1, …, 18; (t is the time step; t is the time; x is the length; τ_σ is the collision relaxation time of σ component, associating with the kinematic viscosity of fluids; f_i^σ(eq)(x,t) is the equilibrium distribution function of σ component.

Fig. 3 Space of discrete internal velocities e_i in D3Q19 LBM.

The equilibrium distribution function is given by:^46,47,54

where d₀ represents the portion of rest fluid particles at equilibrium when u = 0. For component σ, (c_s,σ)² = (1−d₀)/2 = c_s²/m_σ;^46,47 c_s,σ is the lattice sound velocity in the region of pure component σ; c_s is the lattice sound velocity, c_s² = 1/3 in D3Q19 model; m_σ is the molecular mass of component σ; and n_σ is the number density of component σ.

The other main equations of the multicomponent LBM are listed in Table 2. For component σ, the density, kinematic viscosity and velocity can be calculated by the eqn (3)–(5). In the multicomponent LBM, the interaction forces have effects on the equilibrium velocity (eqn (6)) and change the equilibrium distribution function. The total force formula of component σ (eqn (8)) including three part: the interaction force between two components (eqn (9)), the interaction force between fluids and solid surfaces (eqn (11)), and the gravity force (eqn (12)). In eqn (9), ψ_σ(x) is the effective number density, which the function of the number density n_σ(x), and it can be simplified as ψ_σ(x) = n_σ(x). G is the force parameter between two components of the fluids. is called Green function, and is the function of G which depends on the distance between two lattices, as shown in eqn (10). G_σw is the force parameter between the fluid σ and the solid surface. g is the gravitational acceleration. And n_w is the number density of the solid surface. The macro velocity and pressure can be calculated by eqn (13) and (14).

Table 2 Main equations of the multicomponent LBM in D3Q19 model

Descriptions	Equations
Mass density formula of component σ	(3)
Kinematic viscosity of component σ	(4)
Velocity of component σ	(5)
Equilibrium velocity formula of component σ	u^eqσ = u′ + τσFσ/ρσ (6)
Common velocity formula	(7)
Total force formula of component σ	Fσ = F1σ + F2σ + F3σ (8)
Interaction force between different components	(9)
Green function in interaction force	G0 = 0, G1–6 = 2G, G7–18 = G (10)
Interaction force between fluid k and solid surface	(11)
Gravity force formula	F3σ = mσnσg = ρσg (12)
Macro velocity formula	(13)
Macro pressure formula	(14)

---

3 Results and discussion

The experiment and simulation methods were both used for investigating the dynamic stability characteristics of fluid flows in CO₂ miscible displacement in porous media in this study. The experiments were performed at a constantly low injection rate of 0.1 ml min⁻¹ and reservoir conditions of 8.5 to 9.5 MPa and 37.8 °C by using an MRI apparatus. The oil saturation evolution, the in situ mixing zone length, the mixing-frontal and CO₂-frontal velocities, the longitudinal dispersion coefficient were measured and analyzed. Additionally, a three-dimensional LBM was used for simulating miscible displacement flows in micro-pores.

3.1 The dynamic stability characteristics of the in situ mixing zone

The in situ mixing zone represents dynamic characteristics of CO₂ miscible displacement flows, which is important for CO₂-EOR projects.³⁴ In this section, the dynamic stability characteristics, such as the oil-saturation evolution of the in situ mixing zone and the length evolution of the in situ mixing zone, were calculated and analyzed in the miscible displacement processes.

3.1.1 The oil-saturation evolution of the in situ mixing zone. According to the experimental results of Song et al.,⁵⁵ nC₁₀ density evolution is nearly linear with its saturation. The effects of the mixing density can be ignored in CO₂ miscible processes.²¹ The oil saturation evolution can be calculated in proportion to the pseudo-3D MRI signal intensity as follows:^21,29,30

S_o = I_t/I_ini × S_ini × 100% (15)

where I_ini is the initial pseudo-3D MRI signal intensity in the FOV; I_t is the pseudo-3D MRI signal intensity in the FOV at a certain displacement time t; S_ini is the initial oil saturation in the FOV before CO₂ is injected with a value of 1; and S_o is the oil saturation at a certain displacement time t.

In our experiment, the core contained two fluids (oil and CO₂) in the porous media. The CO₂ saturation S_Co2 can be defined as follows:

S_Co2 = 1 − S_o (16)

Fig. 4 shows the oil saturation evolution along the core in the FOV at 8.5, 9.0 and 9.5 MPa. The total displacement time at 8.5 MPa is longer than at 9.0 and 9.5 MPa, meaning the miscible displacement at 8.5 MPa is slower than at 9.0 and 9.5 MPa. The miscible displacement is quicker at a higher pressures.

Fig. 4 The oil saturation evolution along the core at 8.5, 9.0 and 9.5 MPa.

In oil production engineering, the mean oil-saturation evolution in oil fields is important for monitoring oil recovery processes. To quantify the oil-saturation evolution, the mean oil saturation value in the FOV during displacement is used as the oil saturation measured by eqn (15). Fig. 5 shows the mean oil-saturation evolution in the FOV at 8.5, 9.0 and 9.5 MPa during CO₂ miscible displacement processes. The miscible displacement process contains three stages: the first stage O–A, the second stage A–B and the third stage B–C. During the first stage O–A, the oil saturation slowly decreases because the mixing zone gradually goes into the FOV. During the second stage A–B, the saturation curve sharply and linearly decreases. The mixing zone moves from the bottom to the top of the FOV. The curve slope can be used to calculate the estimated velocity using eqn (18) in Sec. 3.2. During the third stage B–C, the oil saturation slowly decreases because the mixing zone gradually leaves the FOV. The residual oil saturation is considered to be the oil saturation at approximately the 100th minute of displacement.

Fig. 5 The mean oil-saturation evolution in the FOV at 8.5, 9.0 and 9.5 MPa in CO₂ miscible displacement processes.

Fig. 5 shows that the first stage O–A at 9.5 MPa is shorter than at 8.5 and 9.0 MPa. During the second stage A–B, the slope of the oil saturation curve at 9.5 MPa is steeper than at 8.5 and 9.0 MPa. The fitting lines of stage A–B at 8.5, 9.0 and 9.5 MPa are shown in Fig. 5. The residual oil saturations at 8.5, 9.0 and 9.5 MPa are 13.7%, 11.6% and 9.5%, respectively. These results indicate that the residual oil saturation decreases with pressure at miscible conditions. Therefore, higher pressures enhance the miscibility^56,57 and reduce the residual oil saturation in CO₂ miscible displacement processes. Similar pressure effects on miscible residual oil were analyzed by Al-Abri et al.⁵⁷

The recovery factor of the CO₂ miscible displacement is estimated using the following equation:

R = (1 − S_or/S_ini) × 100% (17)

where R is recovery factor and S_or is the residual oil saturation. The residual oil saturation is considered to be the oil saturation at the approximately 100th minute. The residual oil saturation and recovery factors of CO₂ miscible displacements from 8.5 to 9.5 MPa in this study are listed in Table 3.

Table 3 Pressure effects on the residual oil saturation and the recovery factors for CO₂ miscible displacements

Pressure (MPa)	Residual oil saturation (%)	Recovery factor (%)
8.5	13.7	86.3
9.0	11.6	88.4
9.5	9.5	90.5

---

3.1.2 The length evolution of the in situ mixing zone. The mixing zone length represents an important one-dimensional diagnostic measure for capturing the overall dynamics of the miscible displacement process.¹³ Song et al.²¹ calculated the mixing zone length in near miscible displacements using MRI technology. According to the measurement method of Song et al.^21,29,30 the mixing zone length of miscible displacements in this study is considered to be the distance between the oil saturation of (S_or + 0.05) and 0.95 in transversely average profiles of MRI images. S_or is the residual oil saturation during miscible displacement processes.

The measurement results of the mixing zone length in the miscible displacement process are shown in Fig. 6. At 8.5 MPa, the length evolution of the in situ mixing zone is decreased at first, and then constant. The length evolution of the in situ mixing zone is nearly constant at 9.0 MPa, but decreased at 9.5 MPa. The decrease phenomenon implies that the mixing zone is more stable at the end rather than the beginning. The average mixing zone lengths at 8.5, 9.0 and 9.5 MPa are 6.8, 6.7 and 6.6 mm, respectively. There is only a slight decrease of the average mixing zone length with pressure, which can be ignored. The core length is 200 mm, the normalized mixing zone lengths at 8.5, 9.0 and 9.5 MPa are approximately 0.03, meaning the mixing zone length is independent of the pressure at the low injection rate of 0.1 ml min⁻¹ above the MMP. The volumetric contraction³⁴ is not obvious in the entire displacement process, which means that the mixture of a CO₂/nC₁₀ system can be approximately seen as an incompressible fluid at the low injection rate of 0.1 ml min⁻¹ with reservoir conditions of 8.5 to 9.5 MPa and 37.8 °C.

Fig. 6 The in situ mixing zone length and its average value at 8.5, 9.0 and 9.5 MPa.

3.2 The dynamic stability characteristics of the mixing front and the CO₂ front

The stability of the mixing front and the CO₂ front is important in CO₂ miscible displacements, which determines CO₂ breakthrough time and oil recovery factor.²¹ In this section, we measured the mixing-frontal velocity and the CO₂-frontal velocity during CO₂ miscible displacement processes, and the critical velocity for stable miscible displacements.

3.2.1 The mixing-frontal velocity and the CO₂-frontal velocity. As defined in our previous work,³⁴ the mixing front is at the position of the transversely average oil saturation 0.95. The CO₂ front is at the position of the transversely average oil saturation (S_or + 0.05). S_or is the residual oil saturation. Then the mixing-frontal and CO₂-frontal positions were calculated. Fig. 7 shows the experimental results of the mixing-frontal and the CO₂-frontal positions with the displacement time at 8.5, 9.0 and 9.5 MPa during CO₂ miscible displacement processes. The fitting lines of the mixing front and the CO₂ front at 8.5 MPa, 9.0 MPa and 9.5 MPa are given. It is shown that the slopes of the mixing front and the CO₂ front are nearly the same at this low injection rate of 0.1 ml min⁻¹.

Fig. 7 The mixing-frontal and CO₂-frontal positions with displacement time at 8.5, 9.0 and 9.5 MPa during the CO₂ miscible displacement process.

Because the mixing-frontal location is approximately proportion to displacement time,^58–60 an averaged frontal velocity can be obtained from experimental data fitting. Table 4 lists the mixing-frontal and the CO₂-frontal velocities for the in situ mixing zone. It is observed that the mixing-frontal and the CO₂-frontal velocities (v₁ and v₂) are nearly the same at the low injection rate of 0.1 ml min⁻¹. And they both increased with pressure.

Table 4 Migration velocities at 8.5, 9.0 and 9.5 MPa

Pressure (MPa)	v1^a (10^−5 m s^−1)	v2^b (10^−5 m s^−1)	v̄^c (10^−5 m s^−1)	ṽ^d (10^−5 m s^−1)	Δṽ^e (10^−5 m s^−1)
a v1: the mixing-frontal velocity.b v2: the CO2-frontal velocity.c v̄: the mean migration velocity, v̄ = (v1 + v2)/2.d ṽ: the estimated velocity calculated by eqn (18).e Δṽ: the absolute error of the estimated velocity, Δṽ = ṽ − v̄.
8.5	1.09	1.12	1.105	0.895	−0.21
9.0	1.91	1.92	1.915	1.58	−0.335
9.5	2.17	2.30	2.235	1.855	−0.38

---

Zhao et al.³⁰ proposed a simple and convenient way to estimate the migration velocity, here it is called the estimated velocity ṽ. The equation of the estimated velocity is as follows:^21,29,30

ṽ = −dS/dt × h (18)

where ṽ is the estimated velocity, m s⁻¹, dS/dt is the slope of the linear curve of the oil saturation with the displacement time during the second stage A–B (Fig. 5), s⁻¹, and h is the migration length of the mixing zone in the FOV, m. Here, h is chosen to be 30 mm.

In this study, the estimated velocity ṽ given by Zhao et al.³⁰ was compared with the real velocity v̄. The estimated velocity ṽ was calculated by eqn (18). And the comparison results are shown in Table 4. The estimated velocity ṽ is always lower than the real velocity v̄. And the maximum absolute error of the estimated velocity method is −0.38 × 10⁻⁵ m s⁻¹. This indicates that the estimated velocity method given by Zhao et al. is always underestimated.

3.2.2 The frontal shapes and the critical velocity of CO₂/nC₁₀ system. The mixing fronts (the frontier of the red region) and the CO₂ fronts (the frontier of the blue region) at 8.5 to 9.5 MPa are shown in Fig. 8. The shapes of the mixing fronts and the CO₂ fronts are practically linear at 8.5 MPa. It indicates that the vertical CO₂ miscible displacement is stable at 8.5 MPa. At this condition, viscous fingering is discounted because of the gravity stability of the vertically upward displacement.^61,62 Because oil is displaced by a lesser dense fluid, CO₂, with a low injection rate in vertically upward displacements, gravity can suppress viscous fingering.^63,64 However, when the pressure increased to 9.0 MPa and 9.5 MPa, the fronts are likely to become inclined and unstable.

Fig. 8 CO₂ frontal shapes for 16 mm MRI slices in miscible displacements with different pressures: (a) 8.5 MPa; (b) 9.0 MPa; and (c) 9.5 MPa.

According to the definition, the critical velocity is the maximum velocity for the displacement to obtain stable fronts. At 8.5 MPa, the mean migration velocity is 1.105 × 10⁻⁵ m s⁻¹, as shown in Table 4. When the migration velocity is larger than this value (1.915 × 10⁻⁵ at 9.0 MPa or 2.240 × 10⁻⁵ m s⁻¹ at 9.5 MPa), the front becomes unstable. It indicates that the critical velocity of CO₂/nC₁₀ system is actually 1.105 × 10⁻⁵ m s⁻¹. In addition, the gravity-stable flood has close relation with the length of the mixing zone. When the flood is no longer stable, the finger front occurs and the mixing zone will be lengthened.³⁴ Therefore, the mixing zone length must be lengthened at unstable conditions of 9.0 and 9.5 MPa, which are above the critical velocity. In contrast, Fig. 8 shows that the mixing zone region is compressed and shortened to some extent at higher pressures of 9.0 and 9.5 MPa. Subsequently, the thinner area of the mixing zone region accelerates the unstability of the migration, and lengthens the length of the in situ mixing zone. As a result, the unstability and the pressure stress couple together and lead to no obvious change of the whole mixing zone length as shown in Fig. 6 in Sec. 3.1. However, the shape of the mixing front under a high pressure is changed and inclined above the critical velocity.

Dumore⁶⁵ gave a method for judging the critical velocity in vertical displacements. The critical velocity criterion is as follows,⁶⁵

where p_o and p_Co2 are the pressures of the displaced fluid (oil) and the displacement fluid (CO₂), respectively, Pa; μ_o and μ_Co2 are viscosities of the displaced fluid and the displacement fluid, Pa s; v_o and v_Co2 are real velocities of the displaced fluid and the displacement fluid, m s⁻¹; ρ_o and ρ_Co2 are densities of the displaced fluid and the displacement fluid, kg m⁻³; z is the displacement distance, m; ϕ is the porosity, 1; K is permeability, m²; and g is the gravitational acceleration, m s⁻². On the right side of eqn (19), ‘+’ is for downward displacements, ‘−’ is for upward displacements.

If , then the displacement is stable. When , the critical rate v_c is obtained. The equation of the critical velocity v_c is as below,

where v_c is critical velocity, m s⁻¹; Δρ = |ρ_o − ρ_s|; and Δμ = |μ_o − μ_s|. In this study, the parameter values were obtained from Table 1. The calculated critical velocity v_c is 1.55 × 10⁻⁴ m s⁻¹. However, according to our experimental results, the critical velocity of CO₂/nC₁₀ system is actually 1.105 × 10⁻⁵ m s⁻¹. Therefore, the critical velocity criterion given by Dumore is about a magnitude larger than the real value.

3.3 The evolution of the longitudinal dispersion coefficient

In recent years, non-Fickian behaviors on miscible displacements have been identified for some conditions.^66,67 The Fickian law has its limitations. To fully compare the Fickian theory, the miscible displacement experiment of a CO₂/nC₁₀ system was fitted with an error function based on the convection–dispersion equation. Additionally, the longitudinal dispersion coefficients and the Peclet numbers were quantified from the fitting results.

The error function used in the convection–dispersion equation is as follows:^68,69

where C is the CO₂ concentration, mol m⁻³; C_ini is the initial CO₂ concentration, mol m⁻³; K_l is the longitudinal dispersion coefficient, m² s⁻¹; x is the position along the core (in the FOV), m; v̄ is the mean migration velocity (see Table 4) of the mixing zone, m s⁻¹; and t is the displacement time, s.

The miscible displacement of the CO₂/nC₁₀ system could assume that a part of the oil (the residual oil S_or) remains in the core without mixing, suggesting the residual oil S_or is a non-reaction portion throughout the miscible process. The error function is then converted to:

where S_Co2 is the CO₂ saturation, S_Co2 = 1 − S_o, where S_o is the oil saturation; and S_or is the residual oil saturation.

Therefore, the graphs of CO₂ saturation can be fitted using the following error function:

To fully investigate the deviation between the miscible experiment of the CO₂/nC₁₀ system and the Fickian theory, CO₂ saturation vs. core position x at a series of fixed displacement time t were fitted using eqn (23). It was fitted a with two parameters, the longitudinal dispersion coefficient K_lx and displacement time t. Fig. 9 shows the evolutions of the longitudinal dispersion coefficients at 8.5, 9.0 and 9.5 MPa. The evolution of the longitudinal dispersion coefficient slightly decreases during displacement. And with proceeding, the longitudinal dispersion coefficient (K_lx) tends to be stable and constant. The declining rate of K_lx increased with pressure at a constant low injection rate of 0.1 ml min⁻¹, meaning the time for K_lx to balance and stay constant decreases with pressure. This is because the increase of pressure can speed up the migration velocity and shorten the equilibrium time for the longitudinal dispersion coefficient to become stable. The range of the longitudinal dispersion coefficients K_lx at 8.5 to 9.5 MPa with a constant injection rate of 0.1 ml min⁻¹ is 8.6 × 10⁻¹⁰ m² s⁻¹ to 4.0 × 10⁻⁹ m² s⁻¹. In our previous work,³⁴ the range of the longitudinal dispersion coefficients K_lx at injection rates of 0.1 to 0.2 ml min⁻¹ with a constant pressure of 8.5 MPa was 8.6 × 10⁻¹⁰ m² s⁻¹ to 3.5 × 10⁻⁸ m² s⁻¹. The evolutions of the longitudinal dispersion coefficients at 8.5 to 9.5 MPa are small and slow compared to the evolutions of the longitudinal dispersion coefficient at injection rates of 0.1 to 0.2 ml min⁻¹. This implies that the longitudinal dispersion coefficient for miscible displacements relies more on injection rate than pressure.

Fig. 9 The evolutions of the longitudinal dispersion coefficient K_lx at 8.5, 9.0 and 9.5 MPa with a constant injection rate of 0.1 ml min⁻¹.

3.4 Lattice-Boltzmann simulation of pore-scale miscible displacement flows

The LBM can simulate fluid flow through complex boundaries and interfacial dynamics,^70,71 which is an effective way to solve simulation limitations of fluid flows for computational fluid dynamics (CFD) in microscopic porous media. In addition, the LBM model is not only suitable for capturing microscopic phenomena, but also satisfies macroscopic mass conservation, momentum conservation and energy conservation equations,³⁷ which may help us to understand macro fluid flow behaviors. In this section, a multicomponent LBM was adopted for simulating miscible displacement flows in the micro-pores. Fig. 10 shows a pore-scale geometric configuration of the regular homogeneous porous media in the LBM simulation. The porosity is 67.03%. In Fig. 10a, the red particles are represented solid parts of porous media, and the yellow cube is extracted for the LBM simulation. Because the shapes of the yellow cube are the same in three directions (x, y and z), the geometry size for the LBM simulation is shown in Fig. 10b. The large space is called pore, and the narrow space is called throat. Fig. 10c is the geometric result in LBM simulation. The blue part (whose solid density is 0) represents the pore space, and the red part (whose solid density is 1) represents the matrix. We assume that CO₂ is injected in z direction at the bottom surface of the cube, and the outlet is at the upper surface of the cube. The study only focuses on the mass transfer in CO₂ injected (z) direction, and ignores the mass transfer in the x–y planes. To enhance the flow in z direction, please note that other four surfaces of the cube are assumed to have no mass transfer which can be regarded as walls.

Fig. 10 The geometric configuration of the regular homogeneous porous media in the LBM simulation.

The parameters for the LBM simulation in the miscible displacement flow are listed in Table 5. The initial mass density ratio and kinematic viscosity ratio of CO₂ to oil can be calculated by eqn (3) and (4). In order to study the effect of viscosity ratio, two groups of relaxation times for CO₂ and oil are chosen for simulations: the first group is 1.0 and 1.0, and the second group is 0.6 (CO₂) and 2.0 (oil). Therefore, the kinematic viscosity ratios of CO₂ to oil in Group 1 and Group 2 are respectively 1 : 1 and 1 : 15, according to eqn (3) and (4).

Table 5 The parameters for the LBM simulation in the miscible displacement flow

Parameter		Value
Injection velocity at z direction, vinj		0.02
Porosity, ϕ		67.03%
Force between two components of fluids, G		−0.002
Force between fluid k and solid surface, Gσw		0
Number density of solid surface, nw		1
Gravity, g		−0.0001
CO2 molecular mass, mCo2		1.0
Oil molecular mass, mo		1.6
Initial CO2 number density, nCo2		1.0
Initial oil number density, no		1.0
Initial mass density ratio of CO2 to oil, ρCo2 : ρo		1 : 1.6
Group 1	Relaxation time ratio of CO2 to oil, τCo2 : τo	1.0 : 1.0
	Kinematic viscosity ratio of CO2 to oil, νCo2 : νo	1 : 1
Group 2	Relaxation time ratio of CO2 to oil, τCo2 : τo	0.6 : 2.0
	Kinematic viscosity ratio of CO2 to oil, νCo2 : νo	1 : 15

---

Fig. 11 is the LBM simulation results of micro-pore miscible displacement flows with the kinematic viscosity ratio of CO₂ to oil τ_Co2 : τ_o = 1 : 1. The concentration of CO₂ is given in the color bar. The blue color means no CO₂. And the red color means 100% CO₂. The results show that the displaced fluid is completely driven out by the displacing fluid, and ultimately the pores are fully occupied with the displacing fluid (CO₂). At the beginning of the displacement, the displacement front is as a piston in the throat. When the displacement moves into the pore, the front fingering occurs, and the fingering length gradually increases. The CO₂ breakthrough occurs at about 600 step, when CO₂ concentration of the upper surface of the cube is above 90%. The displacement is completely finished at 1120 step. The displaced fluid is completely driven out by the displacing fluid, and the ultimate CO₂ concentration is 1. Therefore, the ultimate oil recovery factor is 100%.

Fig. 11 Simulation results of Group 1 in miscible displacement flows with τ_Co2 : τ_o = 1 : 1 in micro-pores. There are front fingerings and no backflow. The displaced fluid is completely driven out by the displacing fluid. The ultimate concentration of the supercritical CO₂ is 1.

Fig. 12 is simulation results of agravic miscible displacement flows with τ_Co2 : τ_o = 0.6 : 2.0 and g = 0 in micro-pores. The kinematic viscosity ratio of CO₂ to oil is 1 : 15. In this case, the high kinematic viscosity of the oil restrains the buoyancy of supercritical CO₂ and makes the displacement more stable with a piston front. There is no CO₂ breakthrough occurred. And from step 600 to the end, a backflow is found in upward vertical displacements. Therefore, the oil is no longer driven by the supercritical CO₂, which distinctly reduces the oil recovery factor. The ultimate CO₂ concentration is 0.4. And the oil recovery factor is 40%. It indicates that a large kinematic viscosity ratio can greatly impede the displacement and may result in a pore-scale backflow phenomenon.

Fig. 12 Simulation results of Group 2 in miscible displacement flows with τ_Co2 : τ_o = 0.6 : 2.0 in micro-pores. The high kinematic viscosity of the oil restrains the buoyancy of supercritical CO₂ and a backflow occurs. It makes the displacement more stable with a piston front. There is no CO₂ breakthrough occurred. The oil cannot be driven out by the supercritical CO₂ entirely, which distinctly reduces the oil recovery factor. The ultimate concentration of the supercritical CO₂ is 0.4.

Fig. 13 is simulation results of agravic miscible displacement flows with τ_Co2 : τ_o = 0.6 : 2.0 and g = 0 in micro-pores. A backflow phenomenon is still occurred without gravity, which indicates that the backflow does not result from the gravity force. The ultimately displacement time without gravity is 1040 step in Fig. 13 which is shorter than that of 1050 step under gravity in Fig. 12. It indicated that the agravic miscible displacement is slightly quicker than the miscible displacement with gravity in upward miscible displacements.

Fig. 13 Simulation results of Group 2 in agravic miscible displacement flows with τ_Co2 : τ_o = 0.6 : 2.0 and g = 0 in micro-pores. A backflow phenomenon is still occurred without gravity, which indicates that the backflow does not results from the gravity force.

4 Conclusions

The dynamic stability characteristics of fluid flow in CO₂ miscible displacements were investigated at a low injection rate of 0.1 ml min⁻¹ and reservoir conditions of 8.5 to 9.5 MPa and 37.8 °C. The oil saturation evolution, the in situ mixing zone length, and the mixing-frontal and the CO₂-frontal velocities were visually quantified using a magnetic resonance imaging apparatus.

The experimental results showed that the residual oil saturation decreased with pressure and the mixing zone length was independent of pressure. The mixing-frontal velocity and CO₂-frontal velocity were nearly the same and increased with pressure at the low injection rate of 0.1 ml min⁻¹. The estimated velocity method given by Zhao et al. was generally underestimated, and the maximum absolute error was −0.38 × 10⁻⁵ m s⁻¹. The critical velocity of CO₂/nC₁₀ system was 1.105 × 10⁻⁵ m s⁻¹. The critical velocity criterion given by Dumore was about a magnitude larger than this value. Although the whole mixing zone length had no obvious change with pressure, a higher pressure compressed the mixing zone and led to an unstable mixing front above the critical velocity.

The longitudinal dispersion coefficient was calculated by fitting the experimental data with an error function. It slightly decreased during the displacement and tended to be constant at the end of the displacement in the porous media. The declining rate increased with pressure, because that the time for K_lx to balance and stay constant decreased with pressure. However, this pressure effect on the longitudinal dispersion coefficient was small and slow. The decreased evolution of the longitudinal dispersion coefficient during the displacement even occurred below the critical velocity. This may be because of the micro-scale heterogeneity of the porous media.

Additionally, a three-dimensional lattice-Boltzmann method was used to simulate pore-scale miscible fluid flows in upward vertical displacements. At a low kinematic viscosity ratio (ν_Co2 : ν_o = 1 : 1), a front fingering was found and a CO₂ breakthrough occurred quickly. At a large kinematic viscosity ratio (ν_Co2 : ν_o = 1 : 15), the displacement was stable with a piston front. The high kinematic viscous oil restrained the buoyancy of supercritical CO₂, but also impeded the displacement with a pore-scale backflow which might lead to a low oil recovery factor. It was concluded that the displacement could be impeded at a large kinematic viscosity ratio and might resulted in a pore-scale backflow with a low oil recovery factor. This study contributes to an increased understanding of the dynamic stability characteristics of fluid flow in CO₂ miscible displacement processes and the mixing mechanism in porous media under reservoir conditions of elevated pressures and high temperatures.

Nomenclature

cs	The lattice sound velocity, cs^2 = 1/3 in D3Q19 model;
cs,σ	The lattice sound velocity in the region of pure component σ in the LBM;
C	The CO2 concentration;
d0	The portion of rest fluid particles at equilibrium when u = 0 in the LBM;
ei	The discrete internal velocities;
i	0, 1, …, 18;
fi^σ(x,t)	The density distribution function of σ component moving in lattice direction ei at node x and time t;
fi^σ(eq)(x,t)	The equilibrium distribution function of σ component;
F	The total interaction force in the LBM, F = F1 + F2 + F3;
F1	The interaction force between different components;
F2	The interaction force between fluid σ and solid surface;
F3	The gravity force;
g	The gravitational acceleration, m s^−2;
G	The parameter in the interaction force in the LBM;
	Green function, is the function of G as shown in eqn (10);
Gσw	The parameter of the force between the fluid σ and the solid surface;
h	The migration length of the mixing zone, m, here h = 30 mm;
I	The MRI signal intensity;
k	The component;
K	The permeability of the core, m^2;
Kl	The longitudinal dispersion coefficient, m^2 s^−1;
Klx	The longitudinal dispersion coefficient by fitting CO2 saturation vs. core position x at a series of fixed displacement time t, m^2 s^−1;
m	The molecular mass in the LBM;
n	The number density in the LBM;
p	The pressure, Pa;
P	The macro pressure in the LBM;
R	The recovery factor;
S	The saturation;
Sor	The residual oil saturation;
t	The time;
δt	The time step;
u	The velocity in the LBM;
u′	The common velocity in the LBM;
U	The macro velocity in the LBM;
v	The real velocity of the fluid, m s^−1;
v1	The mixing-frontal velocity, m s^−1;
v2	The CO2-frontal velocity, m s^−1;
v̄	The mean migration velocity, m s^−1, v̄ = (v1 + v2)/2;
ṽ	The estimated velocity calculated by eqn (18), m s^−1;
Δṽ	The absolute error of the estimated velocity, m s^−1, Δṽ = ṽ − v̄;
vc	The critical velocity, m s^−1;
x	The length;
z	The displacement distance, m;
μ	The viscosity, Pa s;
ρ	The density of the fluid, kg m^−3;
υ	The kinematic viscosity, , m^2 s^−1;
ψ(x)	The function of the number density, which is simplified as ψ(x) = n(x);
τ	The collision relaxation time;
ϕ	The porosity;

Superscripts

eq	The equilibrium state;
σ	The component of the fluid;

Subscripts

c	The critical value;
Co2	The supercritical CO2 phase;
i	The number of the discrete internal velocities;
i	0, 1, …, 18 in D3Q19;
ini	The initial value;
o	The oil phase;
or	The residual value;
t	The displacement time;
w	The solid or wall;
σ	The component of the fluid.

Acknowledgements

This study was supported by the National Natural Science Foundation of China (51376033, 51206017, 41004031 and 51106019), the Ph.D. Programs Foundation of the Ministry of Education of China (20100041120039), the Fundamental Research Funds for Central Universities (DUT12JN13, DUT11RC (3) 65, and DUT13LAB21), the Specialized Research Fund for the Doctoral Program of Higher Education (20120041120057), the Ph. D. Startup Foundation of Liaoning Province, China (20121022), the Major National Science and Technology Program (2011ZX05026-004), and the National Natural Science Foundation of China (51206018, 51106018, and 51176023).

References

1. G. Moritis, Oil Gas J., 2006, 17, 37–57.
2. S. A. Shedid, J. Pet. Sci. Eng., 2006, 50, 285–292.
3. F. I. Stalkup Jr, J. Pet. Technol., 1983, 35, 815–826.
4. E. T. Staff, Improved Oil Recovery, OPHRYS, 1990.
5. F. Chung, R. Jones and T. Burchfield, Recovery of viscous oil under high pressure by CO₂ displacement: a laboratory study, 1988.
6. K. Pyo, N. Damian-Diaz, M. Powell and J. Van Nieuwkerk, CO₂ Flooding in Joffre Viking Pool, 2003.
7. J. B. Magruder, L. H. Stiles and T. D. Yelverton, J. Pet. Technol., 1990, 42, 638–644.
8. M. Cao and Y. Gu, Fuel, 2013, 109, 157–166.
9. Y. Song, L. Jiang, Y. Liu, M. Yang, Y. Zhao, N. Zhu, B. Dou and A. Abudula, Int. J. Greenhouse Gas Control, 2012, 10, 501–509.
10. A.-C. Aycaguer, M. Lev-On and A. M. Winer, Energy Fuels, 2001, 15, 303–308.
11. Z. Wang, M. E. Cates and R. T. Langan, Geophysics, 1998, 63, 1604–1617.
12. J. M. Harris, R. C. Nolen-Hoeksema, R. T. Langan, M. Van Schaack, S. K. Lazaratos and J. W. Rector III, Geophysics, 1995, 60, 667–681.
13. T. Perkins and O. Johnston, Soc. Pet. Eng. J., 1963, 3, 70–84.
14. Y. M. Al-Wahaibi, Energy Fuels, 2010, 24, 1813–1821.
15. F. Torabi and K. Asghari, J. Can. Pet. Technol., 2010, 49, 61–68.
16. A. Rashid, K. Liu, T. Sayem, A. Honari, M. B. Clennell, A. Saeedi and X. Wei, Laboratory Investigation of Factors Affecting CO Enhanced Oil and Gas Recovery, 2013.
17. P. Wylie and K. K. Mohanty, SPE Reservoir Eng., 1997, 12, 264–268.
18. P. L. Wylie and K. K. Mohanty, SPE Reservoir Eval. Eng., 1999, 2, 558–564.
19. F. Torabi and K. Asghari, Fuel, 2010, 89, 2985–2990.
20. J. Trivedi and T. Babadagli, J. Can. Pet. Technol., 2010, 49, 22–27.
21. Y. Song, N. Zhu, Y. Zhao, Y. Liu, L. Jiang and T. Wang, Phys. Fluids, 2013, 25, 053301.
22. J. J. Tessier and K. J. Packer, Phys. Fluids, 1998, 10, 75–85.
23. Z. Pearl, M. Magaritz and P. Bendel, Transp. Porous Media, 1993, 12, 107–123.
24. J. Mitchell, J. Staniland, R. Chassagne and E. J. Fordham, Transp. Porous Media, 2012, 94, 683–706.
25. E. J. Fernandez, C. A. Grotegut, G. W. Braun, K. J. Kirschner, J. R. Staudaher, M. L. Dickson and V. L. Fernandez, Phys. Fluids, 1995, 7, 468–477.
26. H. J. Vinegar and S. L. Wellington, Rev. Sci. Instrum., 1987, 58, 96–107.
27. S. Berg, S. Oedai, A. Landman, N. Brussee, M. Boele, R. Valdez and K. van Gelder, Phys. Fluids, 2010, 22, 113102.
28. V. Clausnitzer and J. Hopmans, Water Resour. Res., 2000, 36, 2067–2079.
29. Y. Zhao, Y. Song, Y. Liu, L. Jiang and N. Zhu, Pet. Sci., 2011, 8, 183–193.
30. Y. Zhao, Y. Song, Y. Liu, H. Liang and B. Dou, Ind. Eng. Chem. Res., 2011, 50, 4707–4715.
31. L. Jiang, Y. Song, Y. Liu, M. Yang, N. Zhu, X. Wang and B. Dou, Chin. J. Chem. Eng., 2013, 21, 85–93.
32. Y. Liu, Y. Zhao, J. Zhao and Y. Song, Magn. Reson. Imaging, 2011, 29, 1110–1118.
33. S. Yong-Chen, Z. Ning-Jun, L. Yu, Z. Jia-Fei, L. Wei-Guo, Z. Yi, Z. Yue-Chao and J. Lan-Lan, Chin. Phys. Lett., 2011, 28, 096401.
34. Y. Song, W. Yang, D. Wang, M. Yang, L. Jiang, Y. Liu, Y. Zhao, B. Dou and Z. Wang, J. Appl. Phys., 2014, 115, 244904.
35. J. Park, Y. Kim, D. Wi, N. Kang, S. Y. Shin and J. Noh, Computer Animation and Virtual Worlds, 2008, 19, 455–467.
36. M. Yoshino, T. Murayama, A. Matsuzaki and T. Hitomi, J. Fluid Sci. Tech., 2009, 4, 13–24.
37. S. Chen and G. D. Doolen, Annu. Rev. Fluid Mech., 1998, 30, 329–364.
38. A. K. Gunstensen, D. H. Rothman, S. Zaleski and G. Zanetti, Phys. Rev. A, 1991, 43, 4320.
39. R. Holme and D. H. Rothman, J. Stat. Phys., 1992, 68, 409–429.
40. E. Flekkøy, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 47, 4247.
41. X. Shan and G. Doolen, J. Stat. Phys., 1995, 81, 379–393.
42. X. Shan and G. Doolen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 54, 3614.
43. T. Inamuro, M. Yoshino, H. Inoue, R. Mizuno and F. Ogino, J. Comput. Phys., 2002, 179, 201–215.
44. R. Span and W. Wagner, J. Phys. Chem. Ref. Data, 1996, 25, 1509–1596.
45. G. Ersland, J. Husebø, A. Graue, B. Baldwin, J. Howard and J. Stevens, Chem. Eng. J., 2010, 158, 25–31.
46. W. Wang, Z. Liu, Y. Jin and Y. Cheng, Chem. Eng. J., 2011, 173, 828–836.
47. J. Zhang and F. Tian, Europhys. Lett., 2008, 81, 66005.
48. X. Shan and H. Chen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 47, 1815.
49. X. Zhang, A. G. Bengough, J. W. Crawford and I. M. Young, Adv. Water Resour., 2002, 25, 1–8.
50. Y. Qian, D. d'Humières and P. Lallemand, Europhys. Lett., 1992, 17, 479.
51. H. Chen, S. Chen and W. H. Matthaeus, Phys. Rev. A, 1992, 45, R5339.
52. P. L. Bhatnagar, E. P. Gross and M. Krook, Phys. Rev., 1954, 94, 511.
53. N. S. Martys and H. Chen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 53, 743.
54. P. Rinaldi, E. Dari, M. Vénere and A. Clausse, Simulat. Model. Pract. Theor., 2012, 25, 163–171.
55. Y. Song, W. Jian, Y. Zhang, Y. Shen, Y. Zhan, J. Zhao, Y. Liu and D. Wang, J. Chem. Eng. Data, 2012, 57, 3399–3407.
56. D. Yang and Y. Gu, Pet. Sci. Technol., 2005, 23, 1099–1112.
57. A. Al-Abri and R. Amin, Transp. Porous Media, 2010, 85, 743–756.
58. C. Jiao and T. Maxworthy, Exp. Fluids, 2008, 44, 781–794.
59. E. Meiburg and G. Homsy, Phys. Fluids, 1988, 31, 429–439.
60. M. Ruith and E. Meiburg, J. Fluid Mech., 2000, 420, 225–257.
61. W. Schowalter, AIChE J., 1965, 11, 99–105.
62. R. Cardenas, R. Alston, A. Nute and G. Kokolis, J. Pet. Technol., 1984, 36, 111–118.
63. G. M. Homsy, in Disorder and Mixing, Springer, 1988, pp. 185–202.
64. Y. M. Al-Wahaibi, A. H. Muggeridge and C. A. Grattoni, J. Pet. Sci. Eng., 2009, 68, 71–80.
65. J. Dumore, Soc. Pet. Eng. J., 1964, 4, 356–362.
66. B. Bijeljic, P. Mostaghimi and M. J. Blunt, Phys. Rev. Lett., 2011, 107, 204502.
67. B. Bijeljic, A. Raeini, P. Mostaghimi and M. J. Blunt, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2013, 87, 013011.
68. M. Sahimi, Flow and transport in porous media and fractured rock: from classical methods to modern approaches, John Wiley & Sons, 2012.
69. J. G. Seo, PhD thesis, Texas A&M University, 2004.
70. C. U. Hatiboglu and T. Babadagli, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2007, 76, 066309.
71. C. U. Hatiboglu and T. Babadagli, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2008, 77, 066311.

---