Contribution of PoreScale Approach to Macroscale Geofluids Modelling in Porous Media
View this Special IssueResearch Article  Open Access
Fluid Interfaces during ViscousDominated Primary Drainage in 2D Micromodels Using PoreScale SPH Simulations
Abstract
We perform porescale resolved direct numerical simulations of immiscible twophase flow in porous media to study the evolution of fluid interfaces. Using a SmoothedParticle Hydrodynamics approach, we simulate saturationcontrolled primary drainage in heterogeneous, partially wettable 2D porous microstructures. While imaging the evolution of fluid interfaces near capillary equilibrium becomes more feasible as fast Xray tomography techniques mature, imaging methods with suitable temporal resolution for viscousdominated flow have only recently emerged. In this work, we study viscous fingering and stable displacement processes. During viscous fingering, porescale flow fields are reminiscent of Bretherton annular flow, that is, the less viscous phase percolates through the core of a porethroat 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 porescale 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 massexchange between hydraulically reservoirconnected and reservoirdisconnected subphases [11, 12].
Considerable effort has been devoted to studying twophase flow patterns at the length scale of porenetworks, 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 porescale dynamics of fluid interfaces and the mechanisms by means of which they evolve remain poorly understood.
To this end, recent advances in porescale imagingbased characterization methods (see review [26]) that enable the fast visualization of twophase flow at porescale resolution, most notably microscopy imaging of thin micromodels [18, 27–29], Xray 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 porescale. For instance, freeenergy 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 quasiincompressible SmoothedParticle Hydrodynamics (SPH) model [38–41] which incorporates the NavierStokes equations together with the Continuum Surface Stress method [42, 43] to account for the interfacial balance equations. Advantages of SPH in the context of twophase flow in porous media include its meshfree 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 meshbased 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 interfacetracking 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 CFLcondition 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 porescale resolved DNS of twophase flow in 2D partially wettable porous media of particulate microstructure. We perform numerical experiments of saturationcontrolled primary drainage for various magnitudes of and . While much effort has been devoted to studying the porescale distribution of fluid interfaces near capillary equilibrium, porescale numerical studies for viscousdominated flow remain comparatively sparse. Rather than studying the emerging macroscopic displacement patterns we discuss porescale 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 tradeoff 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, quasiincompressible wetting () and nonwetting () fluid phase through the porespace () 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 quasiincompressibility. 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 RankineHugoniot 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 surfaceactive agents. Solid wettability properties are prescribed by setting parameters , , and appropriately.
An alternative expression of the interfacial balance of linear momentum (4) for the fluidfluid 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 YoungLaplace 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 wholefluid 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 wholefluid 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 eventdriven particle dynamics algorithm [51, 52]. Porethroat sizes and grain diameters are normally distributed. Mean values and standard deviations are denoted by and for porethroat 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 porethroat 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 porescale effects related to the length scale of microstructure heterogeneity. Nonetheless, a suitable representative size of the computation domain needs be studied using an indepth 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 2tuples Pa s, respectively. The fluidfluid interfacial tension is chosen N/mm whereas Young’s equilibrium contact angle . The latter is achieved by setting fluidsolid 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 porethroat 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 porethroat size as well as and . For viscous fingering, the porescale 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 E51650. 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 capillarydominated flow, in the following we focus on viscousdominated 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 capillarydominated primary drainage is found to occur at nearconstant 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 reservoiraveraged fluid pressures whereas is the ratio of total fluidfluid 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 porescale resolved SPH simulations, we restrict ourselves to twodimensional problems. The latter implies topological constraints related to connectivity in 2D microstructures as well as a lack of representing depthrelated effects, for example, due to curvature of menisci in the outofplane 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 LandauLevichDerjaguin 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 SaffmanTaylor 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 diffusionlimited aggregation (DLA) [61–63]. In this section, we empirically study the annular porescale 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 coreannular gasliquid 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 closedform expressions for wetting film thicknesses in curved capillary tubes as well [64], the complexity of porenetworks 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 porebody (Figure 3, ). The mean curvature flow assumption clearly does not apply for viscous fingering. In 2D microstructures, dynamic wetting films become hydraulically reservoirdisconnected (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 porenetwork 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 noslip 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 porespace. For viscousdominated 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 onesixth 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 antiDLA 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 porescale 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 BuckleyLeverett 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 reservoirdisconnected 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 porescale 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 porescale 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, reservoirconnected 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 porespace such that is equal to the genus of the porespace 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 porenetwork 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 reservoirdisconnected 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 porescale imaging methods, for example, fast Xray tomography methods, to study the distribution and evolution of fluid interfaces for capillarydominated flow become increasingly mature, their applicability to viscousdominated flow is limited by temporal resolution. We consider hydrodynamic direct numerical simulations a suitable complementary approach for porescale modeling of viscousdominated twophase flow. Using a SmoothedParticle Hydrodynamics model, we performed porescale simulations of primary drainage in partially wettable 2D porous media of particulate microstructure at large . We have demonstrated the use of the meshfree 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 meshfree 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 meshfree 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 BuckleyLeverett 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 porescale 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], porethroat size distributions of computational unit cells are computed on the basis of a Delaunay triangulation of grain center points (Figure 9). The porethroat 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 porethroat size histograms. Resulting porethroat 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 email 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/141, within the SPP 2005 “Opus Fluidum FuturumRheology 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)
References
 R. Lenormand, E. Touboul, and C. Zarcone, “Numerical models and experiments on immiscible displacements in porous media,” Journal of Fluid Mechanics, vol. 189, pp. 165–187, 1988. View at: Publisher Site  Google Scholar
 J. Bear, “Dynamics of Fluids in Porous Media,” Soil Science, vol. 120, no. 2, pp. 162163, 1975. View at: Publisher Site  Google Scholar
 G. R. Jerauld and S. J. Salter, “The effect of porestructure on hysteresis in relative permeability and capillary pressure: porelevel modeling,” Transport in Porous Media, vol. 5, no. 2, pp. 103–151, 1990. View at: Publisher Site  Google Scholar
 O. Van Genabeek and D. H. Rothman, “Macroscopic manifestations of microscopic flows through porous media: Phenomenology from simulation,” Annual Review of Earth and Planetary Sciences, vol. 24, pp. 63–87, 1996. View at: Publisher Site  Google Scholar
 J. C. Muccino, W. G. Gray, and L. A. Ferrand, “Toward an improved understanding of multiphase flow in porous media,” Reviews of Geophysics, vol. 36, no. 3, pp. 401–422, 1998. View at: Publisher Site  Google Scholar
 W. G. Gray and S. M. Hassanizadeh, “Averaging theorems and averaged equations for transport of interface properties in multiphase systems,” International Journal of Multiphase Flow, vol. 15, no. 1, pp. 81–95, 1989. View at: Publisher Site  Google Scholar
 S. M. Hassanizadeh and W. G. Gray, “Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries,” Advances in Water Resources, vol. 13, no. 4, pp. 169–186, 1990. View at: Publisher Site  Google Scholar
 S. Majid Hassanizadeh and W. G. Gray, “Toward an improved description of the physics of twophase flow,” Advances in Water Resources, vol. 16, no. 1, pp. 53–67, 1993. View at: Publisher Site  Google Scholar
 S. M. Hassanizadeh and W. G. Gray, “Thermodynamic basis of capillary pressure in porous media,” Water Resources Research, vol. 29, no. 10, pp. 3389–3405, 1993. View at: Publisher Site  Google Scholar
 V. JoekarNiasar, S. M. Hassanizadeh, and A. Leijnse, “Insights into the relationships among capillary pressure, saturation, interfacial area and relative permeability using porenetwork modeling,” Transport in Porous Media, vol. 74, no. 2, pp. 201–219, 2008. View at: Publisher Site  Google Scholar
 R. Hilfer, “Macroscopic equations of motion for twophase flow in porous media,” Physical Review E: Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics, vol. 58, no. 2, pp. 2090–2096, 1998. View at: Publisher Site  Google Scholar
 R. Hilfer and H. Besserer, “Macroscopic twophase flow in porous media,” Physica B: Condensed Matter, vol. 279, no. 13, pp. 125–129, 2000. View at: Publisher Site  Google Scholar
 M. Fleischmann, D. J. Tildesley, and R. C. Ball, Fractals in the Natural Sciences, Princeton University Press, Princeton, 1990. View at: Publisher Site
 O. I. Frette, K. J. Måløy, J. Schmittbuhl, and A. Hansen, “Immiscible displacement of viscositymatched fluids in twodimensional porous media,” Physical Review E: Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics, vol. 55, no. 3, pp. 2969–2975, 1996. View at: Publisher Site  Google Scholar
 M. Ferer, C. Ji, G. S. Bromhal, J. Cook, G. Ahmadi, and D. H. Smith, “Crossover from capillary fingering to viscous fingering for immiscible unstable flow:Experiment and modeling,” Physical Review E: Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics, vol. 70, no. 1, p. 7, 2004. View at: Publisher Site  Google Scholar
 M. A. Theodoropoulou, V. Sygouni, V. Karoutsos, and C. D. Tsakiroglou, “Relative permeability and capillary pressure functions of porous media as related to the displacement growth pattern,” International Journal of Multiphase Flow, vol. 31, no. 1011, pp. 1155–1180, 2005. View at: Publisher Site  Google Scholar
 V. Berejnov, N. Djilali, and D. Sinton, “Labonchip methodologies for the study of transport in porous media: Energy applications,” Lab on a Chip , vol. 8, no. 5, pp. 689–693, 2008. View at: Publisher Site  Google Scholar
 C. Zhang, M. Oostrom, T. W. Wietsma, J. W. Grate, and M. G. Warner, “Influence of viscous and capillary forces on immiscible fluid displacement: Porescale experimental study in a waterwet micromodel demonstrating viscous and capillary fingering,” ENERGY & FUELS, vol. 25, no. 8, pp. 3493–3505, 2011. View at: Publisher Site  Google Scholar
 M. M. Dias and A. C. Payatakes, “Network models for twophase flow in porous media Part 1. Immiscible microdisplacement of nonwetting fluids,” Journal of Fluid Mechanics, vol. 164, no. 2, pp. 305–336, 1986. View at: Publisher Site  Google Scholar
 M. J. Blunt and H. Scher, “Porelevel modeling of wetting,” Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 52, no. 6, pp. 6387–6403, 1995. View at: Publisher Site  Google Scholar
 J. F. Chau and D. Or, “Linking drainage front morphology with gaseous diffusion in unsaturated porous media: A lattice Boltzmann study,” Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 74, no. 5, Article ID 056304, 2006. View at: Publisher Site  Google Scholar
 A. Ferrari and I. Lunati, “Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy,” Advances in Water Resources, vol. 57, no. 9, pp. 19–31, 2013. View at: Publisher Site  Google Scholar
 H. Liu, Y. Zhang, and A. J. Valocchi, “Lattice boltzmann simulation of immiscible fluid displacement in porous media: Homogeneous versus heterogeneous pore network,” Physics of Fluids, vol. 27, no. 5, Article ID 052103, 2015. View at: Publisher Site  Google Scholar
 A. Ferrari, J. JimenezMartinez, T. L. Borgne, Y. Méheust, and I. Lunati, “Challenges in modeling unstable twophase flow experiments in porous micromodels,” Water Resources Research, vol. 51, no. 3, pp. 1381–1400, 2015. View at: Publisher Site  Google Scholar
 P. Kunz, I. M. Zarikos, N. K. Karadimitriou, M. Huber, U. Nieken, and S. M. Hassanizadeh, “Study of Multiphase Flow in Porous Media: Comparison of SPH Simulations with Micromodel Experiments,” Transport in Porous Media, vol. 114, no. 2, pp. 581–600, 2016. View at: Publisher Site  Google Scholar
 T. Bultreys, W. De Boever, and V. Cnudde, “Imaging and imagebased fluid transport modeling at the pore scale in geological materials: a practical introduction to the current stateoftheart,” EarthScience Reviews, vol. 155, pp. 93–128, 2016. View at: Publisher Site  Google Scholar
 N. K. Karadimitriou and S. M. Hassanizadeh, “A review of micromodels and their use in twophase flow studies,” Vadose Zone Journal, vol. 11, no. 3, 2012. View at: Publisher Site  Google Scholar
 F. Moebius and D. Or, “Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries,” Journal of Colloid and Interface Science, vol. 377, no. 1, pp. 406–415, 2012. View at: Publisher Site  Google Scholar
 R. T. Armstrong and S. Berg, “Interfacial velocities and capillary pressure gradients during Haines jumps,” Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 88, no. 4, 2013. View at: Publisher Site  Google Scholar
 D. Wildenschild and A. P. Sheppard, “Xray imaging and analysis techniques for quantifying porescale structure and processes in subsurface porous medium systems,” Advances in Water Resources, vol. 51, pp. 217–246, 2013. View at: Publisher Site  Google Scholar
 M. J. Blunt, B. Bijeljic, H. Dong et al., “Porescale imaging and modelling,” Advances in Water Resources, vol. 51, pp. 197–216, 2013. View at: Publisher Site  Google Scholar
 S. Berg, H. Ott, S. A. Klapp et al., “Realtime 3D imaging of Haines jumps in porous media flow,” Proceedings of the National Acadamy of Sciences of the United States of America, vol. 110, no. 10, pp. 3755–3759, 2013. View at: Publisher Site  Google Scholar
 T. Bultreys, M. A. Boone, M. N. Boone et al., “Realtime visualization of Haines jumps in sandstone with laboratorybased microcomputed tomography,” Water Resources Research, vol. 51, no. 10, pp. 8668–8676, 2015. View at: Publisher Site  Google Scholar
 S. Schlüter, S. Berg, M. Rücker et al., “Porescale displacement mechanisms as a source of hysteresis for twophase flow in porous media,” Water Resources Research, vol. 52, no. 3, pp. 2194–2205, 2016. View at: Publisher Site  Google Scholar
 A. T. Krummel, S. S. Datta, S. Münster, and D. A. Weitz, “Visualizing multiphase flow and trapped fluid configurations in a model threedimensional porous medium,” AIChE Journal, vol. 59, no. 3, pp. 1022–1029, 2013. View at: Publisher Site  Google Scholar
 S. S. Datta, T. S. Ramakrishnan, and D. A. Weitz, “Mobilization of a trapped nonwetting fluid from a threedimensional porous medium,” Physics of Fluids, vol. 26, no. 2, Article ID 022002, 2014. View at: Publisher Site  Google Scholar
 P. Meakin and A. M. Tartakovsky, “Modeling and simulation of porescale multiphase fluid flow and reactive transport in fractured and porous media,” Reviews of Geophysics, vol. 47, no. 3, Article ID RG3002, 2009. View at: Publisher Site  Google Scholar
 J. J. Monaghan, “Smoothed particle hydrodynamics,” Reports on Progress in Physics, vol. 68, no. 8, pp. 1703–1759, 2005. View at: Publisher Site  Google Scholar  MathSciNet
 A. M. Tartakovsky and P. Meakin, “A smoothed particle hydrodynamics model for miscible flow in threedimensional fractures and the twodimensional RayleighTaylor instability,” Journal of Computational Physics, vol. 207, no. 2, pp. 610–624, 2005. View at: Publisher Site  Google Scholar
 A. M. Tartakovsky and P. Meakin, “Pore scale modeling of immiscible and miscible fluid flows using smoothed particle hydrodynamics,” Advances in Water Resources, vol. 29, no. 10, pp. 1464–1478, 2006. View at: Publisher Site  Google Scholar
 X. Y. Hu and N. A. Adams, “A multiphase SPH method for macroscopic and mesoscopic flows,” Journal of Computational Physics, vol. 213, no. 2, pp. 844–861, 2006. View at: Publisher Site  Google Scholar  MathSciNet
 J. U. Brackbill, D. B. Kothe, and C. A. Zemach, “A continuum method for modeling surface tension,” Journal of Computational Physics, vol. 100, no. 2, pp. 335–354, 1992. View at: Publisher Site  Google Scholar  MathSciNet
 B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti, “Modelling merging and fragmentation in multiphase flows with SURFER,” Journal of Computational Physics, vol. 113, no. 1, pp. 134–147, 1994. View at: Publisher Site  Google Scholar  MathSciNet
 R. Sivanesapillai, H. Steeb, and A. Hartmaier, “Transition of effective hydraulic properties from low to high Reynolds number flow in porous media,” Geophysical Research Letters, vol. 41, no. 14, pp. 4920–4928, 2014. View at: Publisher Site  Google Scholar
 R. Sivanesapillai, Porescale study of nondarcian fluid flow in porous media using smoothedparticle hydrodynamics. Doctoral thesis [Doctoral, thesis], 2016.
 J. Bonet and S. Kulasegaram, “Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations,” Int. J. Numer. Meth. Eng, vol. 47, pp. 1189–1214, 2000. View at: Publisher Site  Google Scholar
 L. CuetoFelgueroso, I. Colominas, G. Mosqueira, F. Navarrina, and M. Casteleiro, “On the Galerkin formulation of the smoothed particle hydrodynamics method,” International Journal for Numerical Methods in Engineering, vol. 60, no. 9, pp. 1475–1512, 2004. View at: Publisher Site  Google Scholar
 J. W. Swegle, D. L. Hicks, and S. W. Attaway, “Smoothed particle hydrodynamics stability analysis,” Journal of Computational Physics, vol. 116, no. 1, pp. 123–134, 1995. View at: Publisher Site  Google Scholar  MathSciNet
 P. S. Kurzeja and H. Steeb, “Variational formulation of oscillating fluid clusters and oscillatorlike classification. I. Theory,” Physics of Fluids, vol. 26, no. 4, Article ID 042106, 2014. View at: Publisher Site  Google Scholar
 R. Sivanesapillai, N. Falkner, A. Hartmaier, and H. Steeb, “A CSFSPH method for simulating drainage and imbibition at porescale resolution while tracking interfacial areas,” Advances in Water Resources, vol. 95, pp. 212–234, 2016. View at: Publisher Site  Google Scholar
 A. Donev, S. Torquato, and F. H. Stillinger, “Neighbor list collisiondriven molecular dynamics simulation for nonspherical hard particles. I. Algorithmic details,” Journal of Computational Physics, vol. 202, no. 2, pp. 737–764, 2005. View at: Publisher Site  Google Scholar
 A. Donev, S. Torquato, and F. H. Stillinger, “Neighbor list collisiondriven molecular dynamics simulation for nonspherical hard particles,” Journal of Computational Physics, vol. 202, no. 2, pp. 765–793, 2005. View at: Publisher Site  Google Scholar
 L. Chen and T. C. G. Kibbey, “Measurement of airwater interfacial area for multiple hysteretic drainage curves in an unsaturated fine sand,” Langmuir, vol. 22, no. 16, pp. 6874–6880, 2006. View at: Publisher Site  Google Scholar
 C. E. Schaefer, D. A. DiCarlo, and M. J. Blunt, “Experimental measurement of airwater interfacial area during gravity drainage and secondary imbibition in porous media,” Water Resources Research, vol. 36, no. 4, pp. 885–890, 2000. View at: Publisher Site  Google Scholar
 M. L. Porter, D. Wildenschild, G. Grant, and J. I. Gerhard, “Measurement and prediction of the relationship between capillary pressure, saturation, and interfacial area in a NAPLwaterglass bead system,” Water Resources Research, vol. 46, no. 8, Article ID W08512, 2010. View at: Publisher Site  Google Scholar
 F. P. Bretherton, “The motion of long bubbles in tubes,” Journal of Fluid Mechanics, vol. 10, pp. 166–188, 1961. View at: Publisher Site  Google Scholar  MathSciNet
 B. Widom, Physics Today, vol. 57, no. 12, pp. 6667, 2004. View at: Publisher Site
 V. Starov, M. Velarde, and C. Radke, Wetting and Spreading Dynamics. Surfactant Science Series, CRC Press, Boca Raton, FL, USA, 2007.
 G. F. Teletzke, H. T. Davis, and L. Scriven, “Wetting hydrodynamics,” Revue de Physique Appliquée, vol. 23, no. 6, pp. 989–1007, 1988. View at: Publisher Site  Google Scholar
 P. G. Saffman and G. Taylor, “The penetration of a fluid into a porous medium or HeleShaw cell containing a more viscous liquid,” Proceedings of the Royal Society A Mathematical, Physical and Engineering Sciences, vol. 245, no. 1242, pp. 312–329, 1958. View at: Publisher Site  Google Scholar  MathSciNet
 P. Meakin, “Diffusioncontrolled cluster formation in 2 6dimensional space,” Physical Review A: Atomic, Molecular and Optical Physics, vol. 27, no. 3, pp. 1495–1507, 1983. View at: Publisher Site  Google Scholar
 L. Paterson, “Diffusionlimited aggregation and twofluid displacements in porous media,” Physical Review Letters, vol. 52, no. 18, pp. 1621–1624, 1984. View at: Publisher Site  Google Scholar
 R. Lenormand, “Liquids in porous media,” Journal of Physics: Condensed Matter, vol. 2, no. S, pp. SA79–SA88, 1990. View at: Publisher Site  Google Scholar
 M. Muradoglu and H. A. Stone, “Motion of large bubbles in curved channels,” Journal of Fluid Mechanics, vol. 570, pp. 455–466, 2007. View at: Publisher Site  Google Scholar
 T. C. Ransohoff and C. J. Radke, “Laminar flow of a wetting liquid along the corners of a predominantly gasoccupied noncircular pore,” Journal of Colloid and Interface Science, vol. 121, no. 2, pp. 392–401, 1988. View at: Publisher Site  Google Scholar
 J. M. Nordbotten, M. A. Celia, and S. Bachu, “Injection and storage of CO2 in deep saline aquifers: Analytical solution for CO2 plume evolution during injection,” Transport in Porous Media, vol. 58, no. 3, pp. 339–360, 2005. View at: Publisher Site  Google Scholar
 S. Yuster, “Some Theoretical Considerations of Tilted Water Tables,” Journal of Petroleum Technology, vol. 5, no. 05, pp. 149–156, 2013. View at: Publisher Site  Google Scholar
 M. Ayub and R. G. Bentsen, “Interfacial viscous coupling: a myth or reality?” Journal of Petroleum Science and Engineering, vol. 23, no. 1, pp. 13–26, 1999. View at: Publisher Site  Google Scholar
 G. Jerauld, H. Davis, and L. Scriven, “Stability Fronts of Permanent Form in Immiscible Displacement,” in Proceedings of the SPE Annual Technical Conference and Exhibition, Houston, Texas. View at: Publisher Site  Google Scholar
 A. Riaz and H. A. Tchelepi, “Linear stability analysis of immiscible twophase flow in porous media with capillary dispersion and density variation,” Physics of Fluids, vol. 16, no. 12, pp. 4727–4737, 2004. View at: Publisher Site  Google Scholar
 M. Rücker, S. Berg, R. T. Armstrong et al., “From connected pathway flow to ganglion dynamics,” Geophysical Research Letters, vol. 42, no. 10, pp. 3888–3894, 2015. View at: Publisher Site  Google Scholar
 S. E. Buckley and M. C. Leverett, “Mechanism of fluid displacement in sands,” Transactions of the AIME, vol. 146, no. 1, pp. 107–116, 2013. View at: Publisher Site  Google Scholar
 S. Berg and H. Ott, “Stability of CO2brine immiscible displacement,” International Journal of Greenhouse Gas Control, vol. 11, pp. 188–203, 2012. View at: Publisher Site  Google Scholar
 S. Motealleh, M. Ashouripashaki, D. di Carlo, and S. Bryant, “Mechanisms of capillarycontrolled immiscible fluid flow in fractionally wet porous media,” Vadose Zone Journal, vol. 9, no. 3, pp. 610–623, 2010. View at: Publisher Site  Google Scholar
 R. Holtzman and E. Segre, “Wettability Stabilizes Fluid Invasion into Porous Media via Nonlocal, Cooperative Pore Filling,” Physical Review Letters, vol. 115, no. 16, Article ID 164501, 2015. View at: Publisher Site  Google Scholar
 K. R. Mecke and D. Stoyan, Statistical Physics and Spatial Statistics, vol. 554, Springer Berlin Heidelberg, Berlin, Heidelberg, 2000. View at: Publisher Site
 K. Mecke and C. H. Arns, “Fluids in porous media: a morphometric approach,” Journal of Physics: Condensed Matter, vol. 17, no. 9, pp. S503–S534, 2005. View at: Publisher Site  Google Scholar
 M. Moura, E.A. Fiorentino, G. Schäfer, and R. Toussaint, “Impact of sample geometry on the measurement of pressuresaturation curves: Experiments and simulations,” Water Resources Research, vol. 51, no. 11, pp. 8900–8926, 2015. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Rakulan Sivanesapillai and Holger Steeb. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.