Abstract
We perform pore-scale resolved direct numerical simulations of immiscible two-phase flow in porous media to study the evolution of fluid interfaces. Using a Smoothed-Particle Hydrodynamics approach, we simulate saturation-controlled primary drainage in heterogeneous, partially wettable 2D porous microstructures. While imaging the evolution of fluid interfaces near capillary equilibrium becomes more feasible as fast X-ray tomography techniques mature, imaging methods with suitable temporal resolution for viscous-dominated flow have only recently emerged. In this work, we study viscous fingering and stable displacement processes. During viscous fingering, pore-scale flow fields are reminiscent of Bretherton annular flow, that is, the less viscous phase percolates through the core of a pore-throat forming a hydrodynamic wetting film. Even in simple microstructures wetting films have major impact on the evolution of fluid interfacial area and are observed to give rise to nonnegligible interfacial viscous coupling. Although macroscopically appearing flat, saturation fronts during stable displacement extend over the length of the capillary dispersion zone. While far from the dispersion zone fluid permeation obeys Darcy’s law, the interplay of viscous and capillary forces is found to render fluid flow within complex. Here we show that the characteristic length scale of the capillary dispersion zone increases with the heterogeneity of the microstructure.
1. Introduction
Assessing the stability and evolution of saturation fronts, or, from a pore-scale point of view, interfaces between immiscible bulk fluid phases, is key with respect to understanding a multitude of subsurface processes. Examples include sequestration of carbon dioxide in geological media, groundwater contamination remediation, or enhanced oil recovery. Depending on the governing capillary number (Ca), viscosity ratio (), properties of the porous microstructure, and boundary conditions, primary drainage results in flow regimes as diverse as viscous fingering, stable displacement, or capillary finger branching [1]. In traditional macroscopic continuum models, a phenomenological extension of Darcy’s law is assumed to be applicable with relative permeability and capillary pressure functions representing constitutive model inputs [2]. Calibration of these constitutive relations in light of a particular flow regime typically renders them nonlinear and hysteretic [3–5]. In an attempt to face the latter, contemporary models acknowledge the role of interfacial areas in hysteresis [6–10] or explicitly account for mass-exchange between hydraulically reservoir-connected and reservoir-disconnected subphases [11, 12].
Considerable effort has been devoted to studying two-phase flow patterns at the length scale of pore-networks, both experimentally [13–18] and numerically [19–25], providing a reliable set of data for the - phase diagram of drainage displacement patterns as introduced by [1]. However, the pore-scale dynamics of fluid interfaces and the mechanisms by means of which they evolve remain poorly understood.
To this end, recent advances in pore-scale imaging-based characterization methods (see review [26]) that enable the fast visualization of two-phase flow at pore-scale resolution, most notably microscopy imaging of thin micromodels [18, 27–29], X-ray computed tomography [30–34], and confocal microscopy [35, 36], have provided valuable insights into the interplay of viscous, capillary, gravitational, and inertial forces constituting the complexity of interface dynamics at the pore-scale. For instance, free-energy driven Haines jumps have been confirmed as dominant displacement mechanism for flow at small capillary numbers [29, 32]. Clearly, these observations deviate from the assumptions underlying generalized Darcy flow. Besides experimental approaches, we consider direct numerical simulations (DNS) to be an important complementary tool for quantitative characterization of multiphase flow in porous media [26, 37].
Our method of choice is a quasi-incompressible Smoothed-Particle Hydrodynamics (SPH) model [38–41] which incorporates the Navier-Stokes equations together with the Continuum Surface Stress method [42, 43] to account for the interfacial balance equations. Advantages of SPH in the context of two-phase flow in porous media include its mesh-free nature. The reproducing kernel interpolation of SPH does not require nodal integration points (particles) to be distributed on grids or meshes. The latter renders spatial discretization of complex pore spaces less computationally expensive as compared to traditional grid or mesh-based methods. Typical SPH methods use an updated Lagrangian approach; that is, particles are advected in space according to their local advection velocity. Hence, nonlinear convection terms are not required to be modeled. The latter further implies that phase indicator fields are advected through particle motion. As a result, implementation of the Continuum Surface Stress method does not require interface-tracking such that SPH constitutes an attractive approach to model formation, coalescence, and fragmentation of fluid interfaces in complex pore spaces. Moreover, the updated Lagrangian formulation simplifies modeling of locally large Reynolds numbers [44, 45]. Disadvantages of SPH include its high computational costs associated with repetitive use of neighbor searching algorithms. In the context of using explicit time integration schemes, time stepping stability criteria such as the CFL-condition may be rather restrictive. Moreover, since uncorrected SPH interpolation stencils lack the Kronecker delta property, the application of essential boundary conditions remains nontrivial [46, 47].
In this work, we present pore-scale resolved DNS of two-phase flow in 2D partially wettable porous media of particulate microstructure. We perform numerical experiments of saturation-controlled primary drainage for various magnitudes of and . While much effort has been devoted to studying the pore-scale distribution of fluid interfaces near capillary equilibrium, pore-scale numerical studies for viscous-dominated flow remain comparatively sparse. Rather than studying the emerging macroscopic displacement patterns we discuss pore-scale flow fields. For viscous fingering we discuss the formation of hydrodynamic wetting films and their impact on the evolution of specific interfacial area and interfacial viscous coupling effects. For stable displacement we discuss the microscopic dispersion of stable saturation fronts (capillary dispersion zone) and the role of microstructure heterogeneity. In an attempt to elucidate the trade-off between modeling assumptions, topological constraints of 2D simulations, and physical meaningfulness, we critically discuss our numerical approach.
2. Methods
2.1. Direct Numerical Simulation Approach
We model isothermal flow of a Newtonian, quasi-incompressible wetting () and nonwetting () fluid phase through the pore-space () of a rigid solid matrix (). Bulk fluid balance equations for mass density and linear momentum in phase domains readrespectively. The Newtonian viscous extra stress tensor reads and denotes local velocity. Dynamic viscosity is considered constant. A superimposed dot denotes the material time derivative. Fluid pressure is related to bulk mass density by a stiff barotropic equation of state such that density fluctuations are small implying quasi-incompressibility. A constant, positive backpressure is included to avoid negative pressures that otherwise lead to a numerical tensile instability which SPH methods are prone to [48].
Immiscibility and solid wettability are accounted for using interfacial balance equations for all Gibbs material interfaces . Interfacial balance equations are expressed in terms of jump conditions using the Rankine-Hugoniot jump operator . The interfacial mass balanceswhere denotes interfacial velocity, imply kinematic coupling along the interfacial normal vector . Interfacial balances of linear momentum readwhere denotes the surface divergence operator and the interfacial Cauchy stress tensor . Interfacial tensions are considered constant which is considered reasonable in the absence of surface-active agents. Solid wettability properties are prescribed by setting parameters , , and appropriately.
An alternative expression of the interfacial balance of linear momentum (4) for the fluid-fluid interface readswhere is twice the mean curvature. The first term in (5) implies interfacial viscous coupling whereas the second term introduces a pressure jump condition according to the Young-Laplace equation [49].
Our DNS approach [41, 45, 50] is based on the Continuum Surface Stress method [42, 43] whereby bulk balance equations are formulated for the whole-fluid domain . The latter is achieved by immersing jump conditions (4) using interface Dirac delta distributions which are compactly supported on points of only. The resulting whole-fluid balance of linear momentum is expressed aswhere the solid surface identity tensor is introduced to remove the otherwise unbalanced contact line stress normal to the solid surface [45, 50], that is, to reduce contact line forces to planes tangent to the solid surface. Discretization of (6) using SPH leads to the intuitive motion equationfor a fluid particle with lumped mass . Upon discretization, internal bulk forces due to viscous diffusion and pressure gradient result in discrete interaction forces ( and ) that act between and its neighboring particles .
2.2. Simulation Setup and Microstructures
The total computation domain (Figure 1) is comprised of wetting phase (WP) and nonwetting phase (NWP) reservoirs denoted by and , respectively, as well as the porous computational unit cell . The unit cell has side lengths mm and mm in direction of the unit vectors and , respectively. The microstructure exhibits -periodicity in direction of such that periodic boundary conditions are applied with respect to . The latter avoids boundary artifacts that otherwise arise due to wall confinement. Initially, the sample is completely saturated by the WP. Saturation levels inside the porous sample are controlled by imposing the velocity of pistons that displace the reservoir fluids relative to . As a result of choosing constant, saturation rates remain constant during drainage. Simulations are halted as soon as the NWP penetrates into the WP reservoir, that is, at breakthrough.
Three different microstructures (Appendix) comprised of disordered, nonoverlapping circles (hereinafter referred to as grains) have been generated using an event-driven particle dynamics algorithm [51, 52]. Pore-throat sizes and grain diameters are normally distributed. Mean values and standard deviations are denoted by and for pore-throat size distribution and and for grain diameter distribution, respectively. Microstructures are constructed such that they only differ by the degree of heterogeneity as measured by and , while mean values and , intrinsic permeabilities, and porosities are in close agreement (Table 1). Mean pore-throat size serves as characteristic length scale for the reference capillary pressure . Given its small size, the computational unit cell does not constitute a representative unit cell. However, clipping the sample width by was observed to yield changes in and of less than . In other words, we consider the domain size large enough to study pore-scale effects related to the length scale of microstructure heterogeneity. Nonetheless, a suitable representative size of the computation domain needs be studied using an in-depth converge analysis.
Varying the piston velocity , which also represents the prescribed specific discharge, we control the capillary number
Primary drainage processes are studied for . Fluid viscosities are varied to prescribe the viscosity ratio:
The simulated set of viscosity ratios is by considering the viscosity 2-tuples Pa s, respectively. The fluid-fluid interfacial tension is chosen N/mm whereas Young’s equilibrium contact angle . The latter is achieved by setting fluid-solid interfacial tensions and following Young’s equation. Since the contact angle is fixed for all simulations, effects related to variations in solid surface wettability are not studied herein. Volumetric forces such as gravity are absent and macroscopic inertial effects are considered negligible, that is, small macroscopic Reynolds numbers. The parameterization of fluid densities is thus expected to play a minor role such that we set . Neither the mean pore-throat size mm nor the interfacial tension is chosen with regard to a particular reservoir system but rather in an attempt to optimize computational costs. Both, for larger values of and smaller values of , time stepping constraints become increasingly restrictive.
SPH particles are initially placed on the vertices of a Cartesian grid with grid spacing . Our choice of numerical resolution depends on mean pore-throat size as well as and . For viscous fingering, the pore-scale flow field is found to resemble an annular flow with wetting films which require a fine spatial resolution to ensure numerical accuracy. Hence, for viscous fingering we choose μm which leads to wetting films being represented by approximately 4 particles in direction of film thickness. For all remaining cases, we choose μm. It was previously shown in [45, 50] that the incorporated choices of can reproduce local pressure, menisci, and velocity profiles with reasonable accuracy. The total number of SPH particles is approximately . Average computation time for a single time step was s using a moderately optimized code running 6 threads in parallel on an Intel® Xeon® CPU E5-1650. The total number of time steps varied between for small and for large capillary numbers.
Animated simulation results showing phase, pressure, and velocity magnitude distributions during primary drainage are available as supplementary information (Movies S01–S09). While these animations demonstrate that the numerical model is capable of representing Haines jumps during capillary-dominated flow, in the following we focus on viscous-dominated flow. Phase distributions at breakthrough (Figure 2(a)) are consistent with the expected displacement patterns (capillary finger branching for sufficiently low , stable displacement at sufficiently large and , and viscous fingering at sufficiently large and ); however, their statistical characterization is not insightful herein due to the limited size of the unit cells. A prevailing feature of the present microstructures is that capillary-dominated primary drainage is found to occur at near-constant macroscopic capillary pressure and linear growth of specific fluid interfacial area (Figure 2(b)) until breakthrough occurs. Macroscopic capillary pressure is computed as the difference in reservoir-averaged fluid pressures whereas is the ratio of total fluid-fluid interfacial area to total unit cell volume. While constant relates to a narrow entry pressure distribution along the percolating capillary finger, linear growth of interfacial area is due to morphological simplicity of the microstructures as previously observed in both modeling and experimental studies of primary drainage in particulate microstructures [18, 53–55].
(a)
(b)
In an attempt to highlight methodological limitations, we discuss the most crucial restrictions of the present simulations. (a) The present model is only applicable to ideal solid surfaces with absent chemical imperfections, surface roughness, or dust particles. In the presence of inhomogeneous solid surfaces a phenomenon referred to as contact line hysteresis has to be taken into account. (b) Given the computational costs associated with pore-scale resolved SPH simulations, we restrict ourselves to two-dimensional problems. The latter implies topological constraints related to connectivity in 2D microstructures as well as a lack of representing depth-related effects, for example, due to curvature of menisci in the out-of-plane direction. (c) When a gas bubble invades a capillary tube that is initially saturated with a WP of nonnegligible viscosity, a residual wetting film is hydrodynamically entrained between gas bubble and channel wall regardless of solid wettability characteristics [56]. The latter constitutes a special case of the Landau-Levich-Derjaguin problem [57]. When the thickness of the wetting film (0.1 m) [58], mutual interaction of meniscus and solid surface separated by the molecularly thin wetting film gives rise to an additional contribution to pressure inside the film referred to as disjoining pressure [57, 59] which has a crucial effect on the stability of molecularly thin wetting films. While the present model reproduces hydrodynamic entrainment of wetting films well within the viscous fingering regime, it does not take into account the effect of disjoining pressure.
3. Results and Discussion
3.1. Viscous Fingering Regime: Hydrodynamic Wetting Films
Viscous fingering, i.e., the Saffman-Taylor instability [60], occurs when a less viscous fluid displaces a more viscous fluid () at sufficiently large capillary numbers. While invasion percolation is considered a suitable stochastic model for capillary fingering, viscous fingering is stochastically reminiscent of diffusion-limited aggregation (DLA) [61–63]. In this section, we empirically study the annular pore-scale flow profiles during viscous fingering and the impact of wetting films on the evolution of specific interfacial area and related dynamic effects. Hydrodynamic entrainment of wetting films at large was previously shown in micromodel experiments of [18].
As the present model does not take into account disjoining pressure, wetting film effects are discussed well within the viscous fingering regime () where films are sufficiently thick. For core-annular gas-liquid flow in straight capillary tubes, i.e., Bretherton flow [56], lubrication theory [57] predicts the thickness of liquid films , where the microscopic capillary number incorporates a characteristic microscopic velocity rather than . While lubrication theory has provided closed-form expressions for wetting film thicknesses in curved capillary tubes as well [64], the complexity of pore-networks is expected to necessitate either laboratory experiments or DNS to study hydrodynamic film entrainment in porous media.
If sufficiently fine spatial resolutions are used the present simulations reveal the underlying annular flow fields. Menisci are of spherical shape near fingertip regions whereas the curvature of wetting films is determined by the solid surface (Figure 3). Transition from spherical tip to wetting film is observed less obvious when the finger percolates through a pore-body (Figure 3, ). The mean curvature flow assumption clearly does not apply for viscous fingering. In 2D microstructures, dynamic wetting films become hydraulically reservoir-disconnected (trapped) when viscous fingers coalescence (Figure 3, ). The latter does not hold for 3D microstructures where flow through wetting films that extend over the surface of the entire pore-network might constitute a relevant transport mechanism [65]. The interface of an entrapped wetting film exhibits a ridge in the wake region of the wetted grain as previously observed in micromodel experiments [18]. Due to viscous coupling, recirculating fluid flow takes place within the ridge (Figure 3, ).
Following [66], we consider our choice of viscosity parameters for a rough representation of a carbon dioxide- (-) water system in deep sedimentary formations. In particular, the modeled NWP has small but nonnegligible viscosity ( mPa s). Annular flow may thus constitute a relevant transport mechanism during supercritical sequestration if local flow rates are sufficiently high. The presence of wetting films implies not only a decrease in the effective channel width, but also a change in apparent boundary conditions. Before drainage the kinematic no-slip boundary condition applies with respect to the percolating WP whereas after drainage the dynamic interfacial viscous coupling condition applies with respect to the percolating NWP. Despite the fact that streamline patterns before and after drainage are rather comparable to each other due to laminar flow, the dynamic boundary condition is expected to affect energy dissipation.
It is intuitively understood that wetting films have considerable effect on the evolution of specific interfacial areas. For viscous fingering at and , the rate of change of is observed to be significantly larger as compared to stable displacement and capillary fingering (Figure 4(a)). The magnitude of at breakthrough is found nearly six times the corresponding value for capillary fingering. The latter difference is expected to be even more pronounced for domain sizes large enough such that fractal branching can be statistically reproduced.
(a)
(b)
In an attempt to quantify the dynamic role of wetting films, we measure the total forces the reservoir pistons have to do work against to displace the NWP into the pore-space. For viscous-dominated flow we expect the total resistance force against drainage to be additively composed of the total skin friction force due to viscous shear between fluid and solid phase, that is, Darcian drag, and the total viscous coupling force acting between WP and NWP, that is, the Yuster effect [67, 68]. Total skin friction and viscous coupling force are evaluated asrespectively, where denotes a vector tangent to the interface . The ratio is a measure for the relative dominance of viscous coupling forces in total resistance. Both for capillary fingering () and for stable displacement (, ) viscous coupling forces are negligible (Figure 4(b)). The latter does not apply for viscous fingering where the relative dominance of increases monotonically with and in a manner qualitatively similar to . At breakthrough, nearly one-sixth of total resistance is due to viscous coupling and this trend is expected to be even more pronounced for finer microstructures, larger domain sizes, and higher saturations. These results suggest that viscous coupling terms should be accounted for in macroscopic models for viscous fingering in porous media since otherwise lubrication effects are lumped into relative permeability functions, which, by definition, are related to momentum exchange between solid and fluid phases only.
3.2. Stable Displacement Regime: Capillary Dispersion Zone
Viscosity stabilizes interfacial perturbations when a more viscous fluid displaces a less viscous fluid () at sufficiently large capillary numbers. The latter gives rise to compact displacement patterns, stochastically referred to as anti-DLA patterns [62, 63], and interfaces that macroscopically appear as being flat. Stable displacement is hence desirable during enhanced oil recovery (EOS) due to optimal sweeping efficiency. In our simulations, compact patterns are observed for and (Figure 2(a)). While large implied that capillarity effects are macroscopically negligible, the latter does not apply at smaller length scales. On the contrary, capillarity greatly affects pore-scale interfacial dynamics within the capillary dispersion zone (CDZ) [69–71].
Disregarding boundary effects, saturation profiles that arise during stable displacement are sigmoidal in shape and, in a rough approximation, upper and lower asymptotes to can be considered the saturation limits and , respectively (Figure 5). Considerable saturation gradients are observed only in the localized CDZ. Since we associate macroscopic saturation gradients with the presence of microscopic interfaces, sigmoidal profiles evidence compact patterns. However, rather than being sharp as one would anticipate on the basis of the Buckley-Leverett equation [72] and its admissible shock wave solutions, capillarity is observed to regularize the shock wave within the CDZ.
(a)
(b)
We define the width of the CDZ as spatial difference between both points where saturation profiles approach the asymptotes. While these asymptotes are observed to be barely sensitive to properties of the microstructures, the width , on the other hand, appears to increase with the degree of microstructure heterogeneity. In an attempt to quantify the latter, we compute the temporal average of during drainage, hereafter denoted by . In order to avoid boundary artifacts, we consider the averaging window for equal to the time frame during which the CDZ is entirely contained within the porous unit cell. The average width of the CDZ is found to increase linearly with standard deviation of grain diameter distribution (Figure 6(a)). In particular, yields a reasonable fit to the present data. We emphasize that drainage rates, interfacial tensions, and viscosities, which are all expected to affect , are kept constant meanwhile. Moreover, microstructures A–C only differ by the degree of heterogeneity expressed in terms of while porosities, intrinsic permeabilities, and mean values and are comparable (Table 1). This suggests that, in addition to permeability, flow velocity, viscosity, and characteristic capillary pressure [73], the capillary dispersion length scale also depends on the second central moment of the grain size distribution (~).
(a)
(b)
Within the CDZ, fluid interfaces and hydraulically reservoir-disconnected WP clusters are subject to the interplay of viscous and capillary forces. This is discussed in terms of the temporal evolution of mean pressure profiles (Figure 6(b)). Mean pressure profiles are found to be continuous, piecewise linear functions composed of two line segments. Each line segment can be intuitively attributed to Darcy flow of WP and NWP, respectively. As the kink point moves downstream with increasing NWP saturation, pressure gradients appear invariant to the degree of saturation. The latter is due to the fact that intrinsic permeability, bulk viscosities, and specific discharge remain constant during drainage. While all of the above is consistent with the macroscopic assumption that fluid permeation obeys Darcy’s law a broad fluid pressure distribution is found within the CDZ which we attribute to transient pore-scale events related to disconnected WP cluster dynamics. The latter is most pronounced when the entire CDZ is contained within the porous unit cell (Figure 6(b), blue area).
We empirically discuss the types of pore-scale events that occur near the saturation front and within the CDZ (Figure 7). At the saturation front, adjacent menisci are observed to overlap, a mechanism referred to as Melrose event [74, 75]. The point of overlap is generally located at some distance downstream of the solid surface whereby the enclosed WP becomes hydraulically disconnected and forms a spherical cap (Figure 7, ). In close vicinity to the saturation front, the local viscous pressure drop is negligible (Figure 6(b)) as the local saturation of the less viscous WP is comparatively high. Due to the resulting local dominance of capillary forces, entry pressure thresholds govern flow near the saturation front; that is, invasion percolation takes place at the saturation front. Besides wetting cap formation, capillary trapping mechanisms, for example, pendular bridge formation, thus contribute to disconnection of the WP near the saturation front. As the saturation front further advances, the local viscous pressure drop across a wetting cluster increases due to higher viscosity of the surrounding NWP. The latter leads to fragmentation of large wetting clusters such as the transition of pendular bridges into multiple wetting caps (Figure 7, ). Wetting caps that are of spatial extent large enough for the pressure distribution to be nonuniform are found to be mobilized due to viscous entrainment and eventual coalescence with other disconnected clusters (Figure 7, ). If the size of a wetting cap falls below a critical threshold size, mobilization due to viscous entrainment does not occur. While these results are consistent with previous reports [36], we emphasize that meaningful predictions for the critical threshold size necessitate that contact line pinning and dynamic contact angle effects are accounted for.
(a)
(b)
(c)
3.3. Critical Discussion on Connectivity in 2D Simulations
Recent studies emphasize the importance of integral geometry [76] and, in particular, use of the Euler characteristic to relate topological properties of microstructures [30, 77] or displacement patterns [34] to macroscopic hydraulic properties. For instance, [34] recently demonstrated the fundamental role of fluid topology in permanent hysteresis effects during drainage and imbibition. The lack of representing topological properties of 3D processes is the main drawback of 2D simulations. Indeed, while the solid matrix of a 3D consolidated porous material is generally composed of a single connected component, the solid matrix in 2D simulations has to be composed of disjoint sets for the fluid phase to percolate. We concisely study the evolution of the Euler characteristic of the total WP to emphasize topological implications of 2D simulations. In 2D, the Euler characteristic can be expressed as , where denotes the number of isolated components and denotes the number of circular holes (Betti Numbers). For further introduction to integral geometry, the reader is referred to the above references. The notation applies in the following.
During stable displacement, reservoir-connected WP as well as NWP remain compact clusters and changes in connectivity only occur within the narrow CDZ. On the one hand, the number of isolated components increases due to wetting cap or pendular cluster formation. On the other hand, the number of circular holes within the WP decreases as fewer grains are enclosed by the WP. The distinct characteristic of stable displacement is that the rate of change of , and with saturation is constant, that is, invariant with respect to the position of the compact saturation front. At initial time , WP saturates the entire pore-space such that is equal to the genus of the pore-space and is one. As saturation levels increase both and are observed to increase linearly with saturation (Figure 8(b)). These results are consistent with experiments reported in [34], where the Euler characteristics of compact displacement patterns that emerge during main imbibition were shown to change linearly with saturation. We therefore argue that 2D simulations of compact displacement patterns exhibit topological properties applicable to 3D compact patterns as well.
(a)
(b)
For viscous fingering, on the other hand, 2D simulations are not topologically representative of 3D processes. In 3D, dynamic wetting films are considered to extend over the surface of the entire pore-network and remain connected to the WP reservoir such that should remain constant during viscous fingering in 3D unless discrete WP clusters snap off. On the contrary, dynamic wetting films that emerge during viscous fingering in 2D become reservoir-disconnected as interfaces coalesce and thus form isolated objects of topological genus one (homeomorphic to a circle). As increases with the number of isolated films, the total number of holes , however, is expected to remain unchanged since grains remain wetted by WP. As a result, we expect to change with the same rate as that of as NWP saturation increases. While our results qualitatively reproduce the expected trends (Figure 8(a)), is found to increase with a slightly higher rate as that of . This implies that slightly decreases which can be attributed to the unstable breakup of dynamic wetting films, that is, isolated wetting films of genus one breakup into multiple clusters of genus zero. For an isolated wetting film to be stable in equilibrium, the pressure inside the film must equal the capillary pressure difference across the curved film meniscus when disjoining pressure effects are neglected. As this is generally not the case since film isolation occurred during viscous entrainment, films tend to break up at later stages of the simulation. Both the isolation and subsequent breakup of dynamic wetting films are 2D phenomena. Nonetheless, the impacts of films on specific interfacial area and viscous coupling as discussed in Section 3.1 are expected to apply regardless of above topological issues.
4. Summary and Open Questions
While pore-scale imaging methods, for example, fast X-ray tomography methods, to study the distribution and evolution of fluid interfaces for capillary-dominated flow become increasingly mature, their applicability to viscous-dominated flow is limited by temporal resolution. We consider hydrodynamic direct numerical simulations a suitable complementary approach for pore-scale modeling of viscous-dominated two-phase flow. Using a Smoothed-Particle Hydrodynamics model, we performed pore-scale simulations of primary drainage in partially wettable 2D porous media of particulate microstructure at large . We have demonstrated the use of the mesh-free model to study hydrodynamic entrainment of thick wetting films during viscous fingering as well as the dispersion of stable saturation fronts.
While the impact of wetting films on the dynamics of viscous fingering was outlined, the identification of capillary numbers where the onset of wetting film formation takes place requires further efforts. As the modeling of thin wetting films requires disjoining pressure forces to be accounted for, the present model was only applied to thick hydrodynamic wetting films well within the viscous fingering regime. While the relative dominance of interfacial viscous coupling forces is observed to monotonically increase during viscous fingering, their impact on relative permeability functions remains open. While mesh-free methods were shown to be a viable approach to study viscous fingering phenomena, significant computational resources and algorithmic optimizations are currently required to render 3D mesh-free simulations feasible.
In contrast to viscous fingering, the 2D topological properties of the compact wetting phase during stable displacement were found qualitatively equivalent to those in 3D. The observation that the width of capillary dispersion zones increases with microstructure heterogeneity is thus considered to apply to 3D processes as well. The latter may thus motivate the regularization of Buckley-Leverett models to account for capillary dispersion. Mobilization, fragmentation and coalescence of disconnected wetting phase clusters were found as processes entirely local to the capillary dispersion zone. Hence, the dispersion zone decisively impacts residual wetting phase distribution in form of disconnected wetting caps and pendular bridges. This implies that the region of interest to study wetting phase entrapment during stable displacement using pore-scale imaging methods can be restricted to the capillary dispersion zone. For an accurate numerical prediction of threshold sizes for wetting phase cluster mobilization in reservoir systems, however, more realistic models that take into account contact angle hysteresis are considered necessary.
Appendix
Hydraulic and Morphological Properties of Microstructures
In analogy to what was done in [50, 78], pore-throat size distributions of computational unit cells are computed on the basis of a Delaunay triangulation of grain center points (Figure 9). The pore-throat size between two neighboring grains is considered equal to the length of the connecting Delaunay edge upon subtracting respective grain radii. Delaunay edges that intersect the system boundaries are excluded during generation of pore-throat size histograms. Resulting pore-throat size histograms and grain diameter histograms (Table 1) are subsequently fitted to the Gaussian distributions and , respectively.
(a)
(b)
(c)
Data Availability
Simulation results and data used to generate figures are available upon request via e-mail to the first author.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Acknowledgments
Holger Steeb and Rakulan Sivanesapillai acknowledge the financial support from the German Science Foundation (DFG) through the Project STE 969/14-1, within the SPP 2005 “Opus Fluidum Futurum-Rheology of Reactive, Multiscale, Multiphase Construction Materials”.
Supplementary Materials
Animations showing WP and NWP distribution, normalized pressure, and velocity magnitude fields during the discussed numerical primary drainage experiments are available as online resources and accompany the electronic version of this manuscript. Animations show capillary fingering at in microstructures A–C (Online Resources 1–3), viscous fingering at in microstructures A–C (Online Resources 4–6), and stable displacement at in microstructures A–C (Online Resources 7–9). (Supplementary Materials)