In the oil development process, an immiscible third-phase slug can be injected to the formation temporarily to assist the water flooding, resulting in a three-phase flow underground. In this work, we study slug-assisted water flooding at the pore scale using the three-phase pseudopotential lattice Boltzmann model. We first briefly describe the three-phase pseudopotential model and propose a concise scheme to set the contact angles of the Janus droplet on the solid wall. Then, we simulate the slug-assisted water flooding process in different porous media structures, i.e., a single pore-throat channel, parallel throats, and a heterogeneous porous medium. The simulation results show that oil recovery can be improved effectively with the addition of the third-phase slug. The addition of the third phase results in much more interfacial interaction between different phases, which helps recover trapped oil in pore corners, narrow throats, and the high permeability zone in the porous medium. Moreover, the injection volume, injection timing, contact angle, and viscosity of the third phase influence the oil recovery in different ways. The injected slug can also be trapped in the porous medium, which may result in formation damage. The study explains the enhanced oil recovery mechanisms of slug-assisted water flooding at the pore scale and provides an effective way to design the injection scheme during industrial production.
I. INTRODUCTION
Multiphase flows with ternary components are ubiquitous in our daily life1 and are of great significance in industrial and environmental processes.2,3 The most common ternary flows may be two-phase flow with miscible or immiscible additives (e.g., the interfacial flows involving surfactants) and immiscible three-phase flow (e.g., the oil–water–gas flow underground). These types of flows often involve more complex interfacial phenomena and dynamics due to the presence of the third component and have aroused significant interest to establish different models that can simulate the relevant processes at pore scale.
Pore scale numerical models mainly include pore network modeling (PNM), conventional computational fluid dynamics (CFD) methods, and the lattice Boltzmann method (LBM).4 Although some ternary flow dynamics can be characterized by the PNM and CFD methods,5–7 the LBM has advantages over others in dealing with complex boundaries, ease of implementation and parallelization, and incorporating microscopic interactions.8 Moreover, it has been widely used in the complex flow modeling, such as thermal flows,9,10 particle flows,11,12 ferrofluid flows,13,14 electro-hydrodynamic flows,15,16 as well as multiphase flows. There are several popular multiphase LBM models, such as the pseudopotential (Shan–Chen) model,17 the color-gradient model,18 the free-energy model,19 and the interface tracking model,20 which are systematically introduced in the book by Huang et al.21 However, in early studies, these models generally were presented in two-phase flow cases, and the rapid increase in ternary LBM work has occurred over just the last five years, especially for three-phase immiscible models.
As we focus on immiscible multiphase flow in this work, we only introduce the two-phase flow model containing three components briefly, taking the surfactant LBM as an example. The earliest surfactant LBM may be the free-energy multiphase model or the phase-field model, in which three distribution functions are introduced to simulate the flow field, interface evolution, and convection–diffusion of the surfactant.22–24 Then, Chen et al. proposed a ternary LBM model based on the pseudopotential model, in which the amphiphilic surfactants are treated as dipoles.25 The model was further developed to study the phenomenology of amphiphilic fluids,26,27 surfactant enhanced aquifer remediation,28 and emulsion flow through porous media.29 As for the color-gradient surfactant model, the LBM scheme is used to evolve the multiphase flow, while the surfactant convection–diffusion equation is solved by the finite difference method. The coupling between the LBM and the finite difference scheme is achieved through the LBM variables and the surfactant equation of state.30–32 Very recently, Liu et al. also showed that the hybrid model can handle contact line dynamics and ultra-low interfacial tension situations very well.33
To the best of our knowledge, the first study on the three-phase LB model was the color-gradient model presented in the Ph.D. thesis of Gunstensen in 1992.34 The numerical model is validated against the analytical solutions of the three-phase pipe flow and was further used to simulate the three-phase separation and flow through porous media. Then, van Kats and Egberts studied the oil–water–gas flow at the pore scale using a similar model,35 and Halliday and co-workers generalized the model to N-phase (N > 2) flow.36–38 Leclaire and co-workers further improved density and viscosity ratios of this model,39,40 which are up to O(1000) and O(100), respectively. Based on Leclaire’s model, Fu et al. studied various problems involving three-phase flow, such as double emulsion formation in a microchannel, modulated liquid mixing, and mass transfer and reaction of dilute solutes in a ternary phase system.41–43 Based on the improved color-gradient model proposed by Liu et al.,44 Yu et al. proposed a versatile lattice Boltzmann model for immiscible ternary fluid flows, which is applicable to the fluids with a full range of interfacial tensions. They further incorporated contact line dynamics in the next paper to achieve the desired contact angles.45 In addition, Xie studied the three-phase viscoelastic fluid flow in complex geometries using this kind of model.46
The other immiscible ternary LB models were proposed slightly after the color-gradient model. Travasso et al. proposed a free-energy model to study the morphology and mechanical properties of sheared ternary mixtures in 2005. Semprebon et al. then developed a notable three-phase free-energy model in 2016, in which the surface tension and contact angle are both tunable.47 Using this model, Wang et al. investigated double emulsion formation in planar flow-focusing microchannels and constructed a three-dimensional formation regime diagram based on the capillary numbers of different phase fluids.48 Combined with the entropic lattice Boltzmann scheme,49 the density ratio of this kind model was improved to the order of O(1000), and the model was successfully used to characterize drop dynamics on liquid-infused surfaces.50 Moreover, Bala et al. incorporated wetting boundaries at solid walls for this entropic ternary model based on forcing and geometric schemes.51
Bao and Schaefer may have given the first three-phase simulation result from the pseudopotential model when they tried to improve the density ratio of the multi-component multi-phase (MCMP) model in 2013;52 however, the third phase was not pure and was a mixture of all three components. To study the meniscus-induced motion of bubble and droplets, Wei et al. extended the MCMP pseudopotential model to an immiscible three-phase system, in which the density ratio was adjusted equivalently by a gravity forcing term.53 Using a similar model, Zhang et al. investigated the three-phase relative permeability of 3D porous media54 and Tang et al. studied the CO2 displacement process at the pore scale.55
Liang et al. proposed an alternative LB model for three-phase flows based on the multicomponent phase-field theory in 2016, in which two LB equations are used to capture the interfaces among three different fluids, and another LB equation is adopted to solve the flow field.56 Then, the wetting boundary conditions were well designed,57,58 similar models were generalized to N-phase flows,59–61 and the density ratio and viscosity ratio were highly improved.8,62,63 Moreover, some phase field models chose to capture the fluid–fluid interface by conservative Allen–Cahn equations and can also handle high density and viscosity contrasts.64 In addition, it is also possible to couple the heat transfer equation in the model to study phase-change phenomena in a ternary system.65 Again, the study of immiscible ternary LB models has been active only over the past 5 years, and various models are well established. However, these models were only applied to study a few fundamental problems to verify the models themselves. From our perspective, another period of rapid growth of studies focusing on successful applications of these models to various industrial scenarios is anticipated.
During the oil development process using the chemical flooding method, a temporary immiscible slug may be injected after water flooding to enhance the oil recovery; gel slugs, gas slugs, and polymer slugs66,67 are possible. As the slug chemicals are generally expensive, the slug cannot be injected continuously. Therefore, it is necessary to investigate the underlying enhanced oil recovery (EOR) mechanisms and provide an optimal injection scheme. In this paper, we use the three-phase pseudopotential model we built previously53 to investigate this problem at the pore scale. Although there are some disadvantages of the pseudopotential model, being a “bottom up” model, it is the most computationally efficient multiphase model to simulate flow through porous media. We first briefly describe the ternary pseudopotential model and present various wetting states of the static Janus droplet on a solid wall. Then, we simulate the slug-assisted water flooding process in a single pore-throat channel, a parallel throat, and a heterogeneous porous medium. Finally, we summarize the flow mechanisms of three-phase flow at the pore scale and find that the injection timing and contact angle of the third phase play important roles in the EOR process.
II. NUMERICAL METHOD
A. Pseudopotential model for immiscible three-phase flow
The immiscible three-phase flows are modeled by a pseudopotential LBM, which was developed in our previous work.53 Three distribution functions are introduced to represent the three immiscible fluids, and each distribution function follows the Lattice Bhatnagar-Gross-Krook (LBGK) evolution equation,68
where fi,σ(x, t) is the density distribution function for the component σ, ei is the discrete velocity along the ith lattice direction, τσ is the relaxation time, and is the equilibrium distribution function. The equilibrium distribution function for the D2Q9 model is
where x is the spatial position; ρσ and u are the macroscopic density and macroscopic velocity vector, respectively; the weight factors ωi are given by ω0 = 0, ω1–4 = 1/9, ω5–8 = 1/36; and cs is the sound speed.
Different phases are separated by the long-range repulsive force between different components in the pseudopotential model. Therefore, there are N-1 interaction forces from other fluids in an immiscible N-phase system for each component, and the interaction force on the σk component is
where σk and σj represent different components and is a parameter that controls the interaction force between these two components.
The solid surface force on each component is calculated using the scheme proposed by Martys and Chen,69
where s(x + eiΔt, t) is an indicator function that is equal to 1 or 0 for a solid or fluid domain node and Gads,σ is a parameter that controls the strength of the interfacial tension between the solid and fluid.
Note that the macroscopic velocity ueq is used instead of u when calculating the equilibrium distribution function for each component with Eq. (2), which is given by
where u′ is a composite macroscopic velocity defined as
The kinematic viscosity can be calculated by , and the pressure in this three-phase MCMP model is
B. Wetting boundary condition
Taking the oil contact angle in water as an example, Young’s equation for computing the contact angle satisfies
where γow, γsw, and γso are the interfacial tensions of the fluid interface and the fluid–solid interfaces, and the subscripts o, w, and s represent the oil phase, water phase, and solid wall, respectively.
As the interfacial tension in the MCMP model is proportional to the strength of the interaction force, Huang et al. proposed a straightforward equation to predict the contact angle by70
where ρo,eq and ρw,dis are the equilibrium main density and the associated dissolved density, respectively.
With a third phase added, we assume that there is a Janus droplet surrounded by the water phase sitting on a smooth wall (Fig. 1). Then, the relationship between the static contact angles and interfacial tensions satisfy the Girifalco–Good relation,
The interfacial tension is proportional to in the pseudopotential model, where ρi is the sum of main and dissolved densities of two different components. For simplicity, assuming that the densities of each phase are identical, we can replace the interfacial tension with the parameters of interaction force in Eq. (9),
Therefore, we can set the contact angles by three steps. First, choose reasonable parameters of the interaction force. Then, specify any two contact angles based on Eq. (8). Finally, the remaining contact angle can be set according to Eq. (10). Here, a concise scheme to set the contact angles is presented as follows: We first set all the densities of each phase and all the relaxation times as 1, and let Gint,tw = Gint,ow = 2. By using this set of parameters, the dissolved density is O(0.01) and can be ignored when we calculate the contact angles in Eq. (8). Then, assign a reasonable value to Gads,w randomly and make Gads,o = Gads,w − cos θow, Gads,t = Gads,w − cos θtw. Finally, calculate the Gint,ot using Eq. (10). However, there are strict limitations for the values of the interaction force parameter in the pseudopotential model to ensure numerical convergence. Another alternative way is to assign all the interaction force parameters in the first step and calculate the remaining contact angle in the last step. The disadvantage of this scheme is that we can only set two contact angles in advance and the last contact angle cannot be set arbitrarily. As we mainly focus on the wettability of the discontinuous phase (oil, the third phase) in water, the model is still suitable to study the slug-assisted water flooding process. Some equilibrium shapes of the Janus droplet for different groups of static contact angles are presented in Fig. 2, which shows that the pseudopotential model can simulate different wetting states in immiscible three-phase flow.
III. RESULTS AND DISCUSSION
In this section, we simulate the slug-assisted water flooding process in different porous media structures. Unless stated otherwise, the densities for all three components are set to 1.0, as well as all the relaxation times. Moreover, the interaction force parameters are all set to 3.0. Zou–He pressure boundaries are applied in the inlet and outlet by specifying the densities of different phases. The sum of the densities of all components in the inlet is maintained at 1.2 per lattice, while that in the outlet is kept at 1.0. The porous medium is saturated with oil initially, and then, the water phase is injected into the porous medium by increasing the water density at the inlet. A slug of the third phase is subsequently injected after water flooding by increasing the density of the third phase, and the water phase is re-injected following the third-phase slug.
A. Slug-assisted water flooding in a single pore-throat channel
We first simulate the water flooding process in a channel with constricted pore-throat structures. The domain size is 229 × 89, and the saturation evolution curve with different oil contact angles is shown in Fig. 3. It can be seen that the oil is displaced quickly by the water, and the residual oil mainly distributes in the pore corners or dead ends. Moreover, the smaller the oil contact angle, the higher the residual oil saturation.
Then, we conduct the simulations of slug-assisted water flooding. The third-phase slug is injected from 3000 to 6000 time steps, and the contact angle of the third phase is the same as the oil contact angle. The saturation evolution is shown in Fig. 4 when the oil contact angle is 60°, and the comparison of the final phase distributions between the water flooding and slug-assisted flooding is also presented. It is notable that the time for the system to achieve the stable state becomes longer with the slug injection. Only 10 000 time steps are needed to reach the stable state in water flooding, while it requires more than 40 000 time steps in slug-assisted flow. This is because the number of phase interfaces increases with each extra phase added, and the flow resistance, for example, the pressure difference across the liquid surface or the capillary forces, also increases consequently. Moreover, the third-phase slug usually cannot be fully recovered and may also be trapped in the porous medium. Finally, the oil recovery in all different wetting conditions is improved after the third-phase slug injection. The main EOR mechanisms in the single channel for the slug-assisted water flooding are that part of the slug is trapped in the corner instead of the oil [Fig. 4(a)] and the multiple interface interactions may pull the oil out of the dead ends [Fig. 4(c)]. Although the oil recovery is improved by only 1% in some cases, this can represent an appreciable amount of oil in industrial applications where the overall geological reserve is huge.
The properties of the third-phase slug may influence the EOR process. We further investigate the effects of injection volume, injection timing, viscosity, and contact angle on the saturation evolution. Keeping other parameters the same, snapshots of the slug-assisted flooding with different injection volumes of the third phase are presented in Fig. 5, in which the injection volume is measured in pore volumes (PV). The corresponding saturation evolution curves are shown in Fig. 6. It is interesting that the larger the injected volume, the less the residual of the third phase in Fig. 6(a). This is because the third phase is more continuous when it has a larger volume and tends to flow as a whole slug that is less likely to be trapped in the corners. However, as shown in Fig. 6(b), the oil recovery is the worst for the flooding with the largest injection volume, for the injected fluid hardly replaces the oil in the corners. Meanwhile, the optimal injection volume for the best oil recovery is 0.24PV, neither the minimum injection volume nor the maximum. Considering the other factors such as the slug cost and formation damage caused by the slug detention, achieving the optimal goal in industrial applications is a complex optimization problem.
Figure 7 presents the saturation evolution curves and final phase distributions for slug-assisted water flooding at different injection timings. It shows that the injection timing has little influence on the slug residual saturation, but the earlier injection can lead to better oil recovery.
Figure 8 presents the saturation evolution curves and final phase distributions for slug-assisted water flooding with different contact angles of the third phase. It shows that the contact angle of the third phase has a great influence on both the slug residual saturation and the oil residual saturation. The smaller the contact angle of the third phase, the larger the slug residual saturation, and the smaller the oil residual saturation. This means that the third phase replaces a lot of trapped oil in the corner when the contact angle is small, as shown in the phase distribution results when the contact angle of the third phase is 30° in Fig. 8. In contrast, when the porous medium is not strongly wetted by the slug phase, the injected slug tends to be fully recovered, as shown in Fig. 8 when the contact angle of the third phase is 120°.
By changing the relaxation time of the third phase, we can obtain the saturation evolution curves and final phase distributions for slug-assisted water flooding with different viscosity ratios (the ratio of slug viscosity to water viscosity), as shown in Fig. 9. To ensure the numerical stability of the model, we only present some results with limited viscosity ratios. However, it can clearly be seen that the slug with larger viscosity not only achieves better oil recovery but also has less third-phase residual left in the corners. This is because it requires more energy to trap or pinch off the more viscous fluid in the corners. The third phase will drag out more oil from the corners before it is trapped.
B. Slug-assisted water flooding in parallel throats
A channel with parallel throats is constructed to study the flow behaviors of slug-assisted water flooding. The domain size is 258 × 62. The third phase is injected from 3000 to 8000 time steps, and the contact angle is the same as that of the oil phase. The final phase distributions of water flooding and slug-assisted flooding are compared in Fig. 10. It shows that the residual oil is mainly trapped in the narrow channels during the water flooding, and the slug improves the oil recovery under different wetting conditions. Moreover, with the addition of an extra phase, the three-phase system involves more interfacial interactions and the time for the three-phase flow to become stable is still longer in the simulation, although we do not show the saturation evolution curves here. In addition, as a price to improve the oil recovery, some slug phase has to stay in the throat.
Similarly, we further investigate the effects of slug properties on the EOR process. The displacement processes of slug-assisted flooding with different injected volumes of the third phase are shown in Fig. 11, and the corresponding saturation evolution curves are presented in Fig. 12. The third phase changes the flow events at the channel crossings; for example, the slug may split at a channel intersection and enter into different parallel throats, which promotes oil recovery. Different from the results in the single channel, the more the injected volume of the slug, the better the oil recovery in the parallel throats. This is because the slug phase cannot enter into all the narrow throats if the slug volume is very small. For example, in Fig. 11(a), the slug splits into two parts in the first channel intersection, while it only enters into the upper larger throat in the next intersection. Then, the oil is trapped in the narrow throat and cannot be displaced out. Even worse, the slug flows backward into the narrow throat after it flows out of the upper throat. Moreover, the slug may also be fully trapped in the porous medium if the injected volume is very small, as shown in Fig. 12(a).
As the local flow events in the intersections have great effects on EOR, the injection timing is another important factor that needs to be studied. The displacement processes of slug-assisted flooding at different injection timings are shown in Fig. 13, and the corresponding saturation evolution curves are presented in Fig. 14. The figures show that the earlier the injection timing, the larger the residual saturation of the third phase. This usually also means better oil recovery. It is because the early appearance of the third phase makes the slug more likely to change local flow events at the channel intersections. For example, in Fig. 13(b), the slug enters into the first narrow channel and pushes out part of the oil. Then, it splits at the second intersection and enters into both the narrow and large throats. In the upper large throat, the slug pinches off the oil band and displaces the oil out. In the narrow throat below, the slug pushes all the residual oil out. However, in Fig. 13(d), the slug only enters into the upper large throat, for there is already a dominant pathway formed during the earlier water flooding.
Figure 15 presents the saturation evolution curves and final phase distributions for slug-assisted water flooding with different contact angles of the third phase. It is clear that the slug with a smaller contact angle can get better oil recovery but larger slug retention. We also compare the saturation evolution for slug-assisted water flooding with different viscosity ratios (M = 1, M = 1.5, M = 2), but we find that the saturation curves are similar in this limited range.
The third phase changes the flow events at the channel crossings, and these events may promote both oil recovery and slug retention. We conduct a large number of simulations under different conditions and summarize four basic flow events, i.e., piston-like displacement, pinch off, trailing drag, and wrapping migration, as the pore-scale EOR mechanisms for slug-assisted flooding in the parallel throats, as shown in Fig. 16. The first type is piston-like displacement, i.e., the third phase will enter the narrow channel containing the trapped oil and displace the oil out, as shown in Fig. 16(a). If the oil band enters earlier than the slug when they enter the same channel, the oil between the solid wall and the third phase may be pinched off, and the front part will be pushed forward, which is shown in Fig. 16(b). Moreover, it is also possible that the slug will drag out a part of the resting oil and move it forward, as shown in Fig. 16(c). In addition, if the oil band and the slug enter the channel at the same time, the third phase may wrap the oil and then migrate, carrying the oil along with it, as shown in Fig. 16(d).
C. Slug-assisted water flooding in a heterogeneous porous medium
The domain size of the heterogeneous porous medium used in this section is 335 × 133, in which the upper part is the high permeability zone and the lower part is the low permeability zone. The permeability ratio of two different zones is 12.40 according to the single phase LB simulation. We will explain the analysis of influencing factors briefly, for some conclusions are similar to those in Sec. III B. The comparison of the water flooding and slug-assisted water flooding in a heterogeneous porous medium is shown in Fig. 17. The third phase is injected from 2000 to 6000 time steps, and the contact angles for the oil phase and the third phase are the same. Although we only present the simulation results under one wetting condition, the conclusions for other wetting conditions are similar. First, the oil recovery is improved by the slug and the stabilization time for the slug-assisted water flooding is still longer. In water flooding, the oil is mainly trapped in the boundary corners, vertical space between solid grains, and low permeability zone, among which the low permeability zone contains the most amount of residual oil. With the slug-assisted flooding, the slug recovers all the oil between the solid grains in the high permeability zone and part of the oil in the low permeability zone. Although the oil in the corners near the inlet is displaced out by the third phase, it accumulates in the corners near the outlet again, which may result from insufficient injection volume or wrong injection timing.
Then, we inject the third phase in slug-assisted water flooding at different timings with different injection volumes; the displacement process and corresponding saturation evolution curves are shown in Figs. 18 and 19. It shows that earlier injection timing helps obtain a better oil recovery. Moreover, a larger injection volume does not necessarily mean that more oil can be displaced out.
Figure 20 presents the slug-assisted water flooding in a heterogeneous porous medium with different contact angles. The oil in the boundary corners and high permeability zone is fully recovered when the contact angle of the third phase is small, as shown in Fig. 20(a). With increasing contact angle, the oil begins to be trapped in the boundary corners and the space between solid grains, along with the third phase fluid. We can see that there is a large amount of residual oil and the slug fluid in the porous media when the contact angle of the third phase is 120°, as shown in Fig. 20(c). Note that the retention saturation of the slug is larger with larger contact angles, which is contrary to previous conclusions in the single channel and parallel throats. This is because there are more droplets of the third phase captured in the porous medium by the Jamin effect when the slug phase is non-wetting.
Figure 21 presents the slug-assisted water flooding in a heterogeneous porous medium with different viscosity ratios. Again, although it shows that the larger viscosity of the third phase leads to better oil recovery and larger amount of slug retention, the viscosity ratio range is very small in this study. It is expected that the low permeability zone could be fully swept if assisted by a high viscosity slug.
In summary, the slug-assisted water flooding method can improve the oil recovery by 3%–8% compared to the conventional water flooding method. Although the same influencing factor may influence the final oil recovery in a different way in different porous medium structures, the wetting conditions and injection timing play dominant roles in the EOR process. Moreover, the viscosity may also be significant in a heterogeneous porous medium. In a real porous medium, it may consist of all three types of structures, and the injection scheme should be optimized according to a specific goal, for example, achieving a high oil recovery by injecting a small volume slug with little slug retention.
IV. CONCLUSION
We simulate slug-assisted water flooding at the pore scale using a three-phase lattice Boltzmann method and summarize the EOR mechanisms in different porous media structures.
First, injecting a third-phase slug during water flooding is an effective way to improve oil recovery. In the single pore-throat channel, the slug mainly helps recover the trapped oil in pore corners. In parallel throats, the slug can change the local flow events at channel intersections and recover the trapped oil in narrow throats by piston-like displacement, pinch off, trailing drag, and wrapping migration. In the heterogeneous porous medium, the oil trapped in the boundary corners, space between solid grains, and low permeability zone could be released under the multiple interfacial interactions of the slug.
Furthermore, the EOR process of slug-assisted flooding is influenced by injection volumes, injection timing, contact angles, and viscosity ratios of the third phase. Earlier injection timing, smaller contact angles, and larger viscosity of the third phase usually lead to better oil recovery. A large amount of the injected phase has to stay in the porous medium, which may result in formation damage. Hence, the injection volume needs to be optimized according to the specific situation.
Finally, the viscosity ratio in the current pseudopotential model is limited to below 5 to ensure numerical stability. Simulation of the slug-assisted water flooding with a high viscosity ratio should be investigated further, especially in heterogeneous porous media.
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
ACKNOWLEDGMENTS
The authors greatly appreciate the financial support of the National Key Research and Development Program of China (Grant No. 2018YFA0702400), the National Science Foundation for Distinguished Young Scholars of China (Grant No. 51625403), and the Important National Science and Technology Specific Projects of China (Grant No. 2016ZX05011-003). B.W. was financially supported by the China Postdoctoral Science Foundation (Grant No. 2020M672176), the Qingdao Postdoctoral Applied Research Project (Grant No. qdyy20190033), and the Fundamental Research Funds for the Central Universities (Grant No. 20CX06029A). Funding to M.C.S. for this study was provided by the US National Science Foundation’s Sustainability Research Network Cooperative Agreement (No. 1444758).