Fluid Interfaces during Viscous-Dominated Primary Drainage in 2D Micromodels Using Pore-Scale SPH Simulations

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 equilibriumbecomesmore feasible as fast X-ray tomography techniquesmature, imagingmethodswith 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.


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][4][5].In an attempt to face the latter, contemporary models acknowledge the role of interfacial areas in hysteresis [6][7][8][9][10] or explicitly account for mass-exchange between hydraulically reservoir-connected and reservoirdisconnected subphases [11,12].
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][28][29], X-ray computed tomography [30][31][32][33][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][39][40][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 twophase flow in 2D partially wettable porous media of particulate microstructure.We perform numerical experiments of saturation-controlled primary drainage for various magnitudes of Ca 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.

Direct Numerical Simulation
Approach.We model isothermal flow of a Newtonian, quasi-incompressible wetting (Ω w ) and nonwetting (Ω n ) fluid phase through the pore-space (Ω f ) of a rigid solid matrix (Ω s ).Bulk fluid balance equations for mass density and linear momentum in phase domains respectively.The Newtonian viscous extra stress tensor reads T   =   (grad u  + grad  u  ) and u  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   = (  )+ 0 such that density fluctuations are small implying quasi-incompressibility.A constant, positive backpressure  0 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 Γ  ∈ {Γ wn , Γ ns , Γ ws }.Interfacial balance equations are expressed in terms of jump conditions using the Rankine-Hugoniot jump operator ⟦ ⋅ ⟧  .The interfacial mass balances where u  denotes interfacial velocity, imply kinematic coupling along the interfacial normal vector n  .Interfacial balances of linear momentum read where div  denotes the surface divergence operator and the interfacial Cauchy stress tensor Π  =   (I − n  ⊗ n  ).Interfacial tensions   are considered constant which is considered reasonable in the absence of surface-active agents.Solid wettability properties are prescribed by setting parameters  wn ,  ns , and  ws appropriately.An alternative expression of the interfacial balance of linear momentum (4) for the fluid-fluid interface Γ wn reads where   = −div  n  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  Ω f = Ω n ∪ Ω w .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 as where the solid surface identity tensor I fs fl I − n fs ⊗ n fs is introduced to remove the otherwise unbalanced contact line stress  wn sin Θ 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 equation for a fluid particle P  ∈ Ω f with lumped mass   .Upon discretization, internal bulk forces due to viscous diffusion and pressure gradient result in discrete interaction forces (F   and F   ) that act between P  and its neighboring particles P  .

Simulation Setup and Microstructures.
The total computation domain (Figure 1) is comprised of wetting phase (WP) and nonwetting phase (NWP) reservoirs denoted by Ω w,res and Ω n,res , respectively, as well as the porous computational unit cell Ω cuc .The unit cell has side lengths   = 20 mm and   = 15 mm in direction of the unit vectors e 1 and e 2 , respectively.The microstructure exhibits   -periodicity in direction of e 2 such that periodic boundary conditions are applied with respect to e 2 .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 U  =   e 1 of pistons that displace the reservoir fluids relative to Ω cuc .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 V  for pore-throat size distribution and   and V  for grain diameter distribution, respectively.Microstructures are constructed such that they only differ by the degree of heterogeneity as measured by V  and V  , 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  ref  fl 2 cos Θ wn  −1  .Given its small size, the computational unit cell Ω cuc does not constitute a representative unit cell.However, clipping the sample width   by 50% was observed to yield changes in V  and V  of less than 10%.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 Ca ∈ {10 −2 , 10 −3 , 10 −4 , 10 −5 }.Fluid viscosities are varied to prescribe the viscosity ratio: The simulated set of viscosity ratios is  ∈ {10 −1 , 10 0 , 10 1 } by considering the viscosity 2-tuples log 10 ( n ,  w ) = (−4, −3), (−3, −3), (−3, −4) Pa s, respectively.The fluid-fluid interfacial tension is chosen  wn = 5 N/mm whereas Young's equilibrium contact angle Θ = 30 ∘ .The latter is achieved by setting fluid-solid interfacial tensions  ws = 0 and  ns =  wn cos Θ 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  f =  n =  w .Neither the mean pore-throat size   ≈ 0.41 mm nor the interfacial tension  wn is chosen with regard to a particular reservoir system but rather in an attempt to optimize computational costs.Both, for larger values of  wn 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 d.Our choice of numerical resolution d depends on mean pore-throat size   as well as Ca 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 d =   /19 = 21.7 m which leads to wetting films being represented by approximately 4 particles in direction of film thickness.For all remaining cases, we choose d =   /14 = 29.5 m.It was previously shown in [45,50] that the incorporated choices of d can reproduce local pressure, menisci, and velocity profiles with reasonable accuracy.The total number of SPH particles is approximately 6 ⋅ 10 5 .Average computation time for a single time step was 0.6 s using a moderately optimized code running 6 threads in parallel on an Intel5 Xeon5 CPU E5-1650.The total number of time steps varied between 1.1 ⋅ 10 6 for small and 1.6 ⋅ 10 5 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 Ca, stable displacement at sufficiently large Ca and  > 1, and viscous fingering at sufficiently large Ca and  < 1); 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  wn (Figure 2(b)) until breakthrough occurs.Macroscopic capillary pressure is computed as the difference in reservoir-averaged fluid pressures   ≡ ⟨ n ⟩ res − ⟨ w ⟩ res whereas  wn 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][54][55].
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 t 0 (2.4 mm × 1.9 mm) t 0 + 0.3 ms t 0 + 0.6 ms t 0 + 0.9 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  ∼ O (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.
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 Ca 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 (Ca ∼ 10 −2 ) 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  ∼   Ca 2/3 m , where the microscopic capillary number Ca m 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,  0 + 0.6 ms).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,  0 + 0.9 ms).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,  0 + 0.9 ms).
Following [66], we consider our choice of viscosity parameters for  = 10 −1 a rough representation of a carbon dioxide-(CO 2 -) water system in deep sedimentary formations.In particular, the modeled NWP has small but nonnegligible viscosity ( n = 0.1mPa s).Annular flow may thus constitute a relevant transport mechanism during supercritical CO 2 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 Ca = 10 −2 and  = 10 −1 , the rate of change of  wn is observed to be significantly larger as compared to stable displacement and capillary fingering (Figure 4(a)).The magnitude of  wn 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.
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 as respectively, where t  denotes a vector tangent to the interface Γ  .The ratio ‖F  wn ‖/(‖F  fs ‖ + ‖F  wn ‖) is a measure for the relative dominance of viscous coupling forces in total resistance.Both for capillary fingering (Ca ≤ 10 −4 ) and for stable displacement (Ca ≥ 10 −2 ,  > 10 0 ) viscous coupling forces are negligible (Figure 4(b)).The latter does not apply for viscous fingering where the relative dominance of   wn increases monotonically with  n and in a manner qualitatively similar to  wn ( n ).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.

Stable Displacement Regime: Capillary Dispersion Zone.
Viscosity stabilizes interfacial perturbations when a more viscous fluid displaces a less viscous fluid ( > 1) 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  > 1 and Ca > 10 −3 (Figure 2(a)).While large Ca 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][70][71].
Disregarding boundary effects, saturation profiles  w ( 1 ) that arise during stable displacement are sigmoidal in shape and, in a rough approximation, upper and lower asymptotes to  w ( 1 ) can be considered the saturation limits   w = 1 and   w ≈ 0, 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.
We define the width   () of the CDZ as spatial difference between both points where saturation profiles  w ( 1 , ) 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 V  of grain diameter distribution (Figure 6(a)).In particular,   ≈ 14V    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 V  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 (∼ V  ).
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 ( 1 , ) (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,  0 ).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,  1 ).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,  2 ).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.

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  w of the total WP to emphasize topological implications of 2D simulations.In 2D, the Euler characteristic can be expressed as  w = N w − L w , where N w denotes the number of isolated components and L w 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 N w 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  w ( n ), N w ( n ) and L w ( n ) with saturation is constant, that is, invariant with respect to the position of the compact saturation front.At initial time  0 , WP saturates the entire pore-space such that L w ( 0 ) is equal to the genus of the porespace and N w ( 0 ) is one.As saturation levels increase both N w and  w 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.
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  w 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 N w increases with the number of isolated films, the total number of holes L w ( n ) = L w ( 0 ), however, is expected to remain unchanged since grains remain wetted by WP.As a result, we expect  w to change with the same rate as that of N w as NWP saturation increases.While our results qualitatively reproduce the expected trends (Figure 8(a)),  w is found to increase with a slightly higher rate as that of N w .This implies that L w 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  wn and viscous coupling F  wn as discussed in Section 3.1 are expected to apply regardless of above topological issues.

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 Ca.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 .Black lines in microstructure plots correspond to Delaunay edges that have been used to generate pore-throat size histograms (Table 1).
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.

Figure 1 :
Figure 1: Schematic diagram of simulation setup at initial time  = 0. Black bars represent displacing pistons.

Figure 2 :
Figure 2: (a) NWP distributions (black) at breakthrough for microstructure C in (Ca, ) feature plane.Solid phase and WP are colored gray and white, respectively.(b) Normalized macroscopic capillary pressure   ( w )/ ref  and specific fluid interfacial area  wn ( w ) for capillarydominated primary drainage in microstructure C at log 10 (Ca, ) = (−5, 0).During the preentry regime,   increases to the entry capillary pressure level while  wn ( w ) decreases as the meniscus comes in contact with the solid phase.During the postentry regime,   exhibits a constant mean value with fluctuations related to spontaneous Haines jumps.

Figure 3 :
Figure 3: Hydrodynamic wetting film entrainment in subdomain of microstructure A for log 10 (Ca, ) = (−2, −1).Solid phase domain is colored gray whereas WP SPH particles are represented by black markers.NWP SPH particles are omitted for improved visibility.Streamlines are colored according to the magnitude of local fluid velocity.

Figure 4 :
Figure 4: (a) Evolution of specific fluid-fluid interfacial area  wn in microstructure C. (b) Ratio of total interfacial viscous force to total resistance force in microstructure C. Data are presented as a function of normalized NWP saturation where   n denotes NWP saturation at breakthrough.

Figure 5 :
Figure 5: Evolving saturation profiles in flow direction for stable displacement in microstructures A (a) and C (b) at log 10 (Ca, ) = (−2, 1).Among the tested microstructures, microstructure A exhibits the narrowest grain diameter distribution and microstructure C exhibits the widest one.Width of vertical slices for computing saturation profiles was chosen Δ 1 =   /20.Red markers indicate residual saturation profile   w ( 1 ) at breakthrough.

Figure 6 :
Figure 6: (a) Normalized average width   of CDZ during stable displacement as a function of standard deviation V  of grain diameter distribution.Blue error bars correspond to standard deviations of   ().(b) Normalized pressure profiles in flow direction for microstructure C at log 10 (Ca, ) = (−2, 1).Solid lines represent average fluid pressure profiles ( 1 , ) within vertical averaging slices of width Δ 1 =   /50.Dotted and dashed lines correspond to lower and upper bounds to fluid pressure distribution and envelop 87% of pressure data (shaded area).

Figure 9 :
Figure 9: Computational unit cells A, B, and C (from (a) to (c)).Black lines in microstructure plots correspond to Delaunay edges that have been used to generate pore-throat size histograms (Table1).

Table 1 :
Overview of microstructure parameters.Intrinsic permeabilities  were calculated using single-phase flow simulations.