Impact of Pore-Scale Characteristics on Immiscible Fluid Displacement

Department of Civil Engineering, University of Akron, Akron, Ohio, USA School of Sustainable Engineering and the Built Environment, Arizona State University, Tempe, Arizona, USA Department of Energy Resources Engineering, Stanford University, Stanford, California, USA Department of Civil and Environmental Engineering, Yonsei University, Seoul, Republic of Korea Department of Civil Engineering, Kyung Hee University, Yongin, Republic of Korea Department of Civil and Environmental Engineering, Hanyang University, Seoul, Republic of Korea


Introduction
Immiscible fluid flow is frequently observed in various engineering applications, such as enhanced oil and gas recovery, geological CO 2 sequestration, and soil remediation. While the invading fluid is injected into porous media through an injection well, the defending fluid in the pore space is displaced until the invading fluid reaches a drainage well. Consequently, multiple preferential flow channels are formed. When a nonwetting invading fluid displaces a wetting defending fluid, called drainage displacement, the pattern of fluid invasion and displacement depends on several parameters, such as the viscosity ratio, interfacial tension, injection rate, and wettability, as well as the characteristics of the porous media [1].
The viscosity ratio M is defined as the ratio of the viscosity of the invading fluid (μ inv ) to that of the defending fluid (μ def ). The capillary number C is defined as the ratio of the viscous force to the capillary force: Capillary number : C = qμ inv AT s cos θ : Here, q (m 3 /s) represents the flow rate of the invading fluid, A (m 2 ) represents the cross-sectional area normal to the fluid direction, T s (N/m) represents the interfacial tension between the invading and defending fluids, and θ represents the contact angle at the interface between the invading and defending fluids on the pore surface.
Previous studies have indicated that three displacement patterns can be distinguished, which are illustrated in the log Mlog C plot shown in Figure 1 [1][2][3]. The boundaries for each displacement zone proposed by Lenormand et al. [1], Zhang et al. [2], and Liu et al. [3] are superimposed on Figure 1. Viscous fingering occurs when a less viscous invading fluid displaces a more viscous defending fluid in porous media (i.e., negative log M value). The displacement pattern is characterized by the formation of tree-like fingers with no loops. The principal force is the viscosity of the defending fluid, because capillary effects in the invading fluid are negligible. Capillary fingering prevails at low capillary numbers where viscous forces in both fluids are negligible compared with capillary force. During displacement, an invading fluid forms fingers, which grow in all directions and eventually form loops that trap the clusters of a defending fluid. Stable displacement occurs at relatively high capillary numbers and viscosity ratios, where the dynamics are primarily controlled by the viscous forces of the fluid invading at a relatively high flow rate. A straight invading front perpendicular to the invading direction is formed, and few trapped clusters of the defending fluid are observed.
After the invading fluid reaches the drainage well and displacement no longer occurs, the final saturation of the invading fluid is known as the displacement efficiency, which is the highest in stable displacement and the lowest in viscous fingering [4].
Numerous methods have been used to study immiscible fluid flows, including numerical simulations, e.g., tube network models, or lattice Boltzmann simulations [1,3,5,6] and laboratory displacement experiments using microfluidic chips [2,4,7,8]. The invasion pattern is affected by several parameters, such as the pore morphology (e.g., pore size distribution and connectivity), fluid property (e.g., viscosity and density), fluid-wall interaction (e.g., wettability), and flow rate. This complexity makes both experiments and numerical simulations challenging [1]. In experiments, it is difficult to control various parameters and observe three-dimensional (3D) pore-scale phenomena. In numerical simulations, solving the Navier-Stokes equation under multiphase flow conditions is computationally demanding. Alternative approaches, such as lattice Boltzmann simulations, have a limited capability to properly model all the fluid properties and implement boundary conditions. Owing to these difficulties, few experimental and numerical studies considering various pore-scale characteristics and fluid properties have been performed in the 3D domain [9]. Tube network models are simple yet effective for simulating fluid invasion and displacement by applying the various parameters. Therefore, the objectives of this study were to provide the displacement efficiency values in the 3D domain and to investigate the effects of the pore-scale characteristics and fluid properties on the displacement pattern.
In this study, we used 3D tube networks to simulate the drainage displacement for wide ranges of capillary numbers and viscosity ratios. We also examined the effects of the tube size distribution and connectivity on the phase boundary, displacement pattern, and displacement efficiency.

Numerical Method and Procedure
2.1. Tube Network Extraction. A sandy sediment recovered from the Mallik 5L-38 site (gas hydrate reservoir) at a depth of 1091 m was selected for X-ray computed tomography (CT) scanning and tube network extraction. Detailed information regarding X-ray CT scanning was presented by Mahabadi et al. [10]. In the present study, the volume of the CT image was 8 mm 3 (2 mm × 2 mm × 2 mm) (Figure 2(a), left). A 3D tube network was extracted from the CT images using the maximal ball algorithm developed by Silin and Patzek [11] ( Figure 2(a), right). The maximal ball algorithm searches for all possible spheres that can be inscribed within a pore. After the search, the big spheres located at the centers of the pore spaces become the pores of the pore-network model, 2 Geofluids and the small spheres inscribed in pore throats become tubes connecting adjacent pores. A detailed description of the method was presented by Dong and Blunt [12]. The 3D tube network model extracted from the X-ray image of in situ sediment represents the real pore structure well, because the pore-scale characteristics, such as the connectivity, pore size, and pore length, are used to generate the tube networks. In this study, a network model consisting of only tubes, called tube network or bond-network model, was used. Thus, no pore volume was considered, and the tube length was the distance between the centers of two neighboring pores. The obtained 3D tube network model consisted of 5775 tubes, with a mean tube-coordination number (connectivity) of cn = 8:0, a mean tube size of μ½R T = 12 μm, and a standard deviation of the logarithmic tube size of σðlog R T Þ = 0:31.

Tube Size Distributions.
Four tube size distributions were used in this study. The average value μ½R T in the tube size distribution and the rank of the individual tubes were kept constant, and the standard deviations μ½R T were σ = 0:1σ 0 , 0:5σ 0 , σ 0 , and 2σ 0 , where σ 0 represents the standard deviation of the original tube network model extracted from the CT images. Lower σ values correspond to more uniform tube size distributions. Figure 2(b) shows the four tube size distributions used in this study. The other parameters of the tube size distribution, such as the coordination number, tube location and direction, and tube length, were kept constant.

Coordination Numbers.
Several studies have used twodimensional (2D) tube networks with a homogeneously distributed pore space and a fixed coordination number, e.g., cn = 4 [1][2][3]. The degree of heterogeneity and the coordination numbers typically tend to be higher for in situ sediments, compared with the regular 2D tube network model. The coordination number of the original tube network extracted from the CT image is cn = 8:0. In this study, the effect of the connectivity was investigated using tube networks with four different coordination numbers (cn = 8:0, 6.2, 5.1, and 3 Geofluids 4.6), which are similar to the cn values employed in other studies [13][14][15][16].
The coordination number is reduced by randomly selecting tubes and removing them from the original tube network to generate tube network models with a lower coordination number while keeping the tube size distribution constant.

Simulation
Algorithm for Fluid Invasion. The algorithm developed by Aker et al. [17] was used in this study. This algorithm considers the dynamics of the temporal evolution of the fluid flow and the time dependency of the interface between the invading and defending fluids. Initially, the tube network model is fully saturated with the defending fluid with viscosity μ def . Then, the invading fluid with viscosity μ inv is injected through an inlet boundary at a constant injection rate. Here, we study the drainage displacement; therefore, the invading fluid is nonwetting and the defending fluid is in the wetting phase. It is assumed that the two fluids are immiscible; thus, an explicit interface exists between the defending and invading fluids.
Consider a cylindrical tube containing both invading and defending fluids. The capillary pressure P c at the interface between the two fluids is given by the Young-Laplace equation: where T s (N/m) represents the interfacial tension, R T (m) represents the tube radius, and θ (°) represents the contact angle. Aker et al. [17] modified the Young-Laplace equation to consider the capillary pressure in hour glass-shaped tubes where the capillary pressure depends on the position of the meniscus inside the tubes. Accordingly, the capillary pressure can be expressed as a function of the meniscus position inside the tube: where p represents the relative position of the meniscus, which ranges from 0 to 1 (0 at the beginning of the tube and 1 at the end of the tube). In this study, the contact angle was set as zero (θ = 0) for a perfect wetting condition. The proposed modification in Young-Laplace equation allows the model to simulate the burst drainage displacement dynamics of a multiphase flow system: nonwetting fluid advances to the narrower part of a tube with increasing capillary pressure. After the meniscus passes the pore throat, the invading fluid abruptly occupies some parts of connected tubes. The volume flux q through a cylindrical tube is calculated using the Washburn equation: where μ ′ represents the effective viscosity and ΔP tube represents the pressure difference between the two sides of a capillary tube. For a cylindrical tube without a meniscus, Equation (6) reduces to a Hagen-Poiseuille equation with μ′ = μ def or μ inv . For an incompressible fluid, the mass-balance equation is i.e., the sum of all fluxes into and out of a given node should be zero. Application of Equation (7) to every internal node in the tube network enables the calculation of the pressures at every internal node via the conjugate gradient method. Thus, the pressure difference across the network can be determined. Because the pressure gradient at a constant flow rate must be obtained, one must determine the nodal pressures at a constant injection rate. The flow rate Q in the tube network model is given as where A and B are parameters that depend on the properties of the fluid and the geometry of the porous medium and ΔP network represents the pressure difference between the inlet and the outlet of the network model. These parameters A and B can be calculated by solving the nodal pressures for two differential pressures ΔP 1 network and Δ P 2 network and corresponding flow rates Q 1 and Q 2 . After determining the coefficients A and B, ΔP network is calculated for the desired constant flow rate. Then, all the nodal pressures are computed by solving Equation (7) for a given pressure difference ΔP network . A detailed explanation was presented by Aker et al. [17].
The timestep Δt is selected such that only one meniscus among the menisci in all tubes in the network model can reach the other end of the tube in each timestep during invasion: Here, ΔL ij represents the distance between the meniscus and the end of the tube, and V ij represents the flow velocity in the tube, which can be calculated using the flow rate and the tube radius: V ij = Q/πR T 2 . At each timestep, the minimum Δt is calculated and is used for the next timestep. The new position of the meniscus is calculated and is updated using a second-order Runge-Kutta scheme.

Typical Characteristics of Invasion Simulation
Typical displacement patterns for selected capillary numbers and viscosity ratios are shown in Figure 3. The simulation results clearly indicate that the displacement patterns were affected by the values of log C and log M. At a low viscosity 4 Geofluids ratio (log M = -5, Figure 3(a)), tree-like displacement patterns and a clear single preferential flow path developed. In this case, the saturation of the invading fluid at the end of the invasion increased slightly, from 0.066 to 0.074, as the capillary number decreased from log C = 5 to log C = -5. At a low capillary number (log C = -5, bottom row of Figure 3), the displacement pattern changed from viscous fingering to capillary fingering as the viscosity ratio increased from log M = -5 to log M = 5. Patterns of stable displacement were observed for a high viscosity ratio (log M = 5) and high capillary number (log C = 5). The final saturation for stable displacement was approximately 0.7 for log M = 5 and log C = 5.
After the invasion process was completed, the tube network model with dimensions of 2 mm × 2 mm × 2 mm was divided equally into five sections along the flow direction. The "local" saturation of the invading fluid in each internal section was calculated for nine values of log M and log C (Figure 4(a)). Depending on the specific phase boundary proposed in the literature, the points (log M, log C) in the phase diagram were categorized into either one of the three displacement zones (viscous fingering, capillary fingering, and stable displacement zone) or the transition zone (between the three displacement-pattern zones). Here, the noticeable observation is that in the intermediate zone at log M = 0, the displacement efficiency decreases as the capillary number increases, not like at the condition of log M = 5.
As shown in Figure 4(b), the local saturations for the nine investigated cases were spatially inhomogeneous and decreased from the inlet to the outlet. For the three cases with log M = 5 (which correspond to capillary fingering when log C = -5, the transition zone for log C = 0, and stable displacement when log C = 5), the saturations of the 1 st and 2 nd internal sections were almost S = 100% and then decreased to approximately 0.5% in the last section near the outlet boundary. As the log capillary number increased from -5 to 0 and 5, the local saturations in the 3 rd and 4 th internal sections increased, resulting in an overall increase in the "global" saturation, because the displacement pattern changed from capillary fingering to stable displacement, through the transition zone. However, for the three cases of log M = 0 (in the transition zone), the "local" saturation for a given internal section decreased as the capillary number increased from -5 to 5. This effect was most pronounced for the 1 st , 2 nd , and 3 rd internal sections, because the displacement pattern changed from capillary fingering to the transition zone in the proximity of the viscous fingering zone. As shown in Figure 5, the saturation was higher in the capillary fingering zone than in the transition zone, particularly near the viscous fingering zone. For the three cases of log M = -5, the "local" saturation started decreasing in the first section and remained almost constant between the second and last sections for a given log C value. Note that the saturation averaged for

Effects of Viscosity Ratio and Capillary Number on Displacement Efficiency
The viscosity ratios and capillary numbers used in this study ranged from log M = -10 to 10 and from log C = -10 to 10, respectively. A total of 441 simulation runs were performed using 21 viscosity ratios and 21 capillary numbers, with the original tube network model. The obtained displacement efficiency with respect to the viscosity ratio and capillary number is shown in a contour plot ( Figure 5). Three displacement zones, i.e., stable displacement (denoted by S), viscous fingering (V), and capillary fingering (C), are identified with regard to the displacement efficiency. A fourth zone (denoted by "V-C," lower left in the phase diagram) is separated from the conventional capillary fingering zone, because the displacement efficiency in the "V-C" zone is significantly lower than that in the capillary fingering zone. The "V-C" zone has an intermediate displacement efficiency between those of the viscous fingering and capillary fingering zones. The displacement efficiency values in the zones divided by the dashed lines in Figure 5 were used to calculate the average value of saturation (displacement efficiency). The displacement efficiencies in the transition zones and near these boundaries were excluded from the calculation. The average value S ave of the displacement efficiency was 0.68 for the stable displacement zone, 0.53 for the capillary fingering zone, and 0.07 for the viscous fingering zone. The average displacement efficiency in the "V-C" zone was S ave = 0:17. In the three displacement-pattern zones, i.e., stable displacement, capillary fingering, and viscous fingering, the coefficient of variation for the displacement efficiency values for nine neighboring log Clog M points (obtained for three consecutive values on each log C and log M axis) was <0.1. However, the "V-C" zone exhibited an increasing trend in the displacement efficiency from the viscous fingering zone to the capillary fingering zone. Furthermore, the boundary between the stable displacement and capillary fingering zones could not be clearly identified according to the displacement efficiency only.

Impact of Pore-Scale Characteristics on Displacement
Additional simulations were performed using three different tube size distributions and three different coordination numbers (cn) to investigate the effects of the pore-scale characteristics on the displacement patterns. Figures 6(a) and 7(a) show contour maps of the displacement efficiency for four different tube size distributions and coordination numbers, respectively. Results for the original tube network model are also presented.

Effect of Tube Size Distribution on Displacement
Efficiency. In Figure 6(a), each contour map contains the results of 441 simulation runs for 21 viscous ratios and 21 capillary numbers. For stable displacement and viscous fingering, the displacement efficiency decreased as the standard deviation of the tube size distribution decreased. However, this trend was not observed in the capillary zone; the tube network with the original standard deviation (σ = σ 0 ) exhibited a saturation of S = 0:53, whereas the values of S were 0.57, 0.58, and 0.59 for σ = 2σ 0 , 0:5σ 0 , and 0:1σ 0 , respectively. Thus, the difference in the displacement efficiency between the "C" and the "S" zones decreased. For these simulations using different tube size distributions (Figure 6), the rank of tube size in the tube network model was maintained constant even though the tube size distribution varies. As the tube size distribution increases, there occur more numbers of larger-and smaller-sized tubes. Due to the same rank of tube size, the number of tubes invaded at the end of the invasion might be similar or the same for the tube network models with different tube size distributions, which results in higher saturation (or displacement efficiency) for the condition of higher tube size distribution. In order words, the same number of tubes is invaded, but they might be larger tubes. However, this phenomenon is not pronounced in the capillary fingering zone. The possible reason might be that multiple tubes with the similar capillary pressures could be invaded simultaneously, which results in higher displacement efficiency for lower tube size distribution condition.
Moreover, as the standard deviation of the tube size distribution decreased, the vertical boundaries separating the viscous fingering and stable displacement zones moved to higher values of log M (from -5 to -2) in the log Mlog C plot. Furthermore, the displacement efficiency in the "V-C" zone decreased and approached saturation value in the "V" zone. Thus, when the tube network model had a uniform tube size distribution, the displacement pattern was determined predominantly by the viscosity of both fluids, rather than by capillarity effects in the tube network. The boundary between the stable displacement and capillary fingering zones became unclear as the difference in the average displacement efficiency between stable displacement and capillary fingering decreased, for a uniform tube size distribution. This trend is highlighted in Figure 6(b), which shows the gradient of the displacement efficiency. High gradients (i.e., large changes in the displacement efficiency) are observed in the transition zone, which is shown in red. For σ = 2σ 0 , high gradients between four displacement-pattern zones (S, V, C, and V-C) are observed, and the boundary between stable displacement and capillary fingering is clearly identifiable. As the standard deviation of the tube size distribution decreases, the boundary between the "S" and "C" zones becomes unclear, and viscosity effects increase.

Effect of Coordination Number on Displacement
Efficiency. Contour maps of the displacement efficiency are shown in Figure 7(a). The simulation results indicate that the displacement efficiency for each displacement pattern (V, S, and C) decreased as the coordination number decreased. However, the displacement efficiency in the "V-C" zone exhibited the opposite trend. As the coordination number decreased, the difference in the displacement efficiencies between the "V-C" and "C" zones decreased. Additionally, the boundary between the viscous fingering and capillary fingering zones became clearer. The shape and location of the boundaries for the lowest coordination number, i.e., cn = 4:6 ( Figure 7(a), right), were very similar to those of the boundaries suggested by Lenormand et al. [1], Liu et al. [3], and Zhang et al. [2], where cn = 4:0. This trend is also observed in Figure 7(b). For a given standard deviation σ = σ 0 , the gradient of the displacement efficiency clearly exhibited a change in the phase boundaries as the coordination number decreased.

Limitations.
The results of this study are applicable to several immiscible fluid flow conditions, such as those in CO 2 sequestration, soil contamination and remediation, fluid control using Pickering emulsions, enhanced oil and gas recovery, and gas hydrate production [4,[18][19][20][21]. However, a few assumptions and limitations must be considered for applications. (1) The tube network is millimeter-scale, in contrast to kilometer-scale in situ applications. This limitation may be resolved by obtaining constant local saturation values using a longer tube network. For field application, as long as the log C and log M values are accurately identified depending on the distance from the input well, the corresponding displacement efficiency can be predicted using the constant local saturation values obtained from the network model simulation. (2) Even in the small tube network, the local saturation varies significantly owing to the capillary end effect. This phenomenon can be minimized by using a longer network or controlled outflow boundary [22]. (3) The shape of real pore spaces is not perfectly cylindrical. The capillary pressure for irregularly shaped pores found in in situ conditions differs from that for perfectly circular pores [23]. (4) The contact angle and interfacial tension vary depending on the in situ conditions [24]. (5) The invading fluid with a high flow rate may affect the configuration of sediment particles (e.g., the redistribution of sediment particles) [25,26]; however, in this study, it was assumed that the pore shape did not change during the invasion. (6) A constant flow rate condition was selected in this study to evaluate the displacement efficiency for a given log C condition, and the results might be different if other boundary conditions are used. Moreover, in in situ applications, the invasion rate may decrease as the  Geofluids invading fluid spreads in the radial direction from an injection well.

Conclusions
The displacement efficiency values of an immiscible fluid flow were obtained for a wide range of log C and log M values in the 3D domain, using tube networks. The effects of the pore-scale characteristics, such as the tube size distribution and tube connectivity, on the displacement efficiency were investigated. The main findings of this study are presented below: (i) The dynamic invasion process resulted in saturation gradients inside the tube network. After invasion is complete, the saturation of the invading fluid is higher near the inlet and decreases toward the outlet (ii) Three displacement patterns were observed, similar to previous studies. Additionally, a fourth "V-C" zone that was previously part of the capillary fingering zone was defined by the displacement efficiency gradually increasing from the viscous fingering zone to the capillary fingering zone (iii) The displacement efficiency obtained via a 3D tube network simulation was S ave = 0:68 for the stable displacement zone, S ave = 0:07 for the viscous fingering zone, S ave = 0:53 for the capillary fingering zone, and S ave = 0:17 for the "V-C" zone (iv) As the tube size distribution became narrower distribution, the boundaries between the stable displacement and capillary fingering zones and between the viscous fingering and "V-C" zones became unclear. Additionally, the boundary between the viscous fingering and stable displacement zones moved to a higher log M value. The displacement efficiency is affected dominantly by the viscosity ratio in the tube network especially when the tube size distribution becomes uniform. The displacement efficiency decreased with the increasing uniformity of the tube size distribution for stable displacement and viscous fingering (v) As the coordination number cn (indicating the tube connectivity) decreased from 8.0 to 4.6, the boundary between the viscous fingering and "V-C" zones became clear. This is similar to the phase boundary reported in the literature. The displacement efficiency decreased with the decreasing coordination number for the stable displacement, viscous fingering, and capillary fingering zones