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

Authors: Rakulan Sivanesapillai; Holger Steeb
Journal: Geofluids
 (2018)

---

## Abstract

We perform pore-scale resolved direct numerical simulations of immiscible two-phase flow in porous media to study the evolution of fluid interfaces. Using a Smoothed-Particle Hydrodynamics approach, we simulate saturation-controlled primary drainage in heterogeneous, partially wettable 2D porous microstructures. While imaging the evolution of fluid interfaces near capillary equilibrium becomes more feasible as fast X-ray tomography techniques mature, imaging methods with suitable temporal resolution for viscous-dominated flow have only recently emerged. In this work, we study viscous fingering and stable displacement processes. During viscous fingering, pore-scale flow fields are reminiscent of Bretherton annular flow, that is, the less viscous phase percolates through the core of a pore-throat forming a hydrodynamic wetting film. Even in simple microstructures wetting films have major impact on the evolution of fluid interfacial area and are observed to give rise to nonnegligible interfacial viscous coupling. Although macroscopically appearing flat, saturation fronts during stable displacement extend over the length of the capillary dispersion zone. While far from the dispersion zone fluid permeation obeys Darcy’s law, the interplay of viscous and capillary forces is found to render fluid flow within complex. Here we show that the characteristic length scale of the capillary dispersion zone increases with the heterogeneity of the microstructure.


---

## Body

## 1. Introduction

Assessing the stability and evolution of saturation fronts, or, from a pore-scale point of view, interfaces between immiscible bulk fluid phases, is key with respect to understanding a multitude of subsurface processes. Examples include sequestration of carbon dioxide in geological media, groundwater contamination remediation, or enhanced oil recovery. Depending on the governing capillary number (Ca), viscosity ratio (M), properties of the porous microstructure, and boundary conditions, primary drainage results in flow regimes as diverse as viscous fingering, stable displacement, or capillary finger branching [1]. In traditional macroscopic continuum models, a phenomenological extension of Darcy’s law is assumed to be applicable with relative permeability and capillary pressure functions representing constitutive model inputs [2]. Calibration of these constitutive relations in light of a particular flow regime typically renders them nonlinear and hysteretic [3–5]. In an attempt to face the latter, contemporary models acknowledge the role of interfacial areas in hysteresis [6–10] or explicitly account for mass-exchange between hydraulically reservoir-connected and reservoir-disconnected subphases [11, 12].Considerable effort has been devoted to studying two-phase flow patterns at the length scale of pore-networks, both experimentally [13–18] and numerically [19–25], providing a reliable set of data for the Ca-M phase diagram of drainage displacement patterns as introduced by [1]. However, the pore-scale dynamics of fluid interfaces and the mechanisms by means of which they evolve remain poorly understood.To this end, recent advances in pore-scale imaging-based characterization methods (see review [26]) that enable the fast visualization of two-phase flow at pore-scale resolution, most notably microscopy imaging of thin micromodels [18, 27–29], X-ray computed tomography [30–34], and confocal microscopy [35, 36], have provided valuable insights into the interplay of viscous, capillary, gravitational, and inertial forces constituting the complexity of interface dynamics at the pore-scale. For instance, free-energy driven Haines jumps have been confirmed as dominant displacement mechanism for flow at small capillary numbers [29, 32]. Clearly, these observations deviate from the assumptions underlying generalized Darcy flow. Besides experimental approaches, we consider direct numerical simulations (DNS) to be an important complementary tool for quantitative characterization of multiphase flow in porous media [26, 37].Our method of choice is a quasi-incompressible Smoothed-Particle Hydrodynamics (SPH) model [38–41] which incorporates the Navier-Stokes equations together with the Continuum Surface Stress method [42, 43] to account for the interfacial balance equations. Advantages of SPH in the context of two-phase flow in porous media include its mesh-free nature. The reproducing kernel interpolation of SPH does not require nodal integration points (particles) to be distributed on grids or meshes. The latter renders spatial discretization of complex pore spaces less computationally expensive as compared to traditional grid or mesh-based methods. Typical SPH methods use an updated Lagrangian approach; that is, particles are advected in space according to their local advection velocity. Hence, nonlinear convection terms are not required to be modeled. The latter further implies that phase indicator fields are advected through particle motion. As a result, implementation of the Continuum Surface Stress method does not require interface-tracking such that SPH constitutes an attractive approach to model formation, coalescence, and fragmentation of fluid interfaces in complex pore spaces. Moreover, the updated Lagrangian formulation simplifies modeling of locally large Reynolds numbers [44, 45]. Disadvantages of SPH include its high computational costs associated with repetitive use of neighbor searching algorithms. In the context of using explicit time integration schemes, time stepping stability criteria such as the CFL-condition may be rather restrictive. Moreover, since uncorrected SPH interpolation stencils lack the Kronecker delta property, the application of essential boundary conditions remains nontrivial [46, 47].In this work, we present pore-scale resolved DNS of two-phase flow in 2D partially wettable porous media of particulate microstructure. We perform numerical experiments of saturation-controlled primary drainage for various magnitudes ofCa and M. While much effort has been devoted to studying the pore-scale distribution of fluid interfaces near capillary equilibrium, pore-scale numerical studies for viscous-dominated flow remain comparatively sparse. Rather than studying the emerging macroscopic displacement patterns we discuss pore-scale flow fields. For viscous fingering we discuss the formation of hydrodynamic wetting films and their impact on the evolution of specific interfacial area and interfacial viscous coupling effects. For stable displacement we discuss the microscopic dispersion of stable saturation fronts (capillary dispersion zone) and the role of microstructure heterogeneity. In an attempt to elucidate the trade-off between modeling assumptions, topological constraints of 2D simulations, and physical meaningfulness, we critically discuss our numerical approach.


## 2. Methods

### 2.1. Direct Numerical Simulation Approach

We model isothermal flow of a Newtonian, quasi-incompressible wetting (Ω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 Ωα∈{Ωn,Ωw} read(1)ϱ˙α=-ϱαdivuα,(2)ϱαu˙α=divTEα-gradpα,respectively. The Newtonian viscous extra stress tensor reads TEα=μα(graduα+gradTuα) and uα denotes local velocity. Dynamic viscosity μα is considered constant. A superimposed dot denotes the material time derivative. Fluid pressure pα is related to bulk mass density ϱα by a stiff barotropic equation of state pα=f(ϱα)+p0 such that density fluctuations are small implying quasi-incompressibility. A constant, positive backpressure p0 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(3)ϱu-uαβ·nαβαβ=0,where uαβ denotes interfacial velocity, imply kinematic coupling along the interfacial normal vector nαβ. Interfacial balances of linear momentum read(4)TE-pI·nαβαβ=-divαβΠαβ,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(5)TEn-TEw·nαβ+pw-pn+σwnκwnnwn=0,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(6)ϱfu˙f=divTEf-gradpf+∑αβdivΠαβδαβ·Ifs,where the solid surface identity tensor Ifs≔I-nfs⊗nfs is introduced to remove the otherwise unbalanced contact line stress σwnsin⁡Θ 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(7)miu˙i=∑jFijV-∑jFijP+∑αβFiαβfor a fluid particle Pi∈Ωf with lumped mass mi. Upon discretization, internal bulk forces due to viscous diffusion and pressure gradient result in discrete interaction forces (FijV and FijP) that act between Pi and its neighboring particles Pj.


### 2.2. Simulation Setup and Microstructures

The total computation domain (Figure1) 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 LW=20 mm and LH=15 mm in direction of the unit vectors e1 and e2, respectively. The microstructure exhibits LH-periodicity in direction of e2 such that periodic boundary conditions are applied with respect to e2. 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 UP=UPe1 of pistons that displace the reservoir fluids relative to Ωcuc. As a result of choosing UP constant, saturation rates remain constant during drainage. Simulations are halted as soon as the NWP penetrates into the WP reservoir, that is, at breakthrough.Figure 1
Schematic diagram of simulation setup at initial timet=0. Black bars represent displacing pistons.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 mT and vT for pore-throat size distribution and mG and vG for grain diameter distribution, respectively. Microstructures are constructed such that they only differ by the degree of heterogeneity as measured by vT and vG, while mean values mT and mG, intrinsic permeabilities, and porosities are in close agreement (Table 1). Mean pore-throat size mT serves as characteristic length scale for the reference capillary pressure PCref≔2cos⁡ΘσwnmT-1. Given its small size, the computational unit cell Ωcuc does not constitute a representative unit cell. However, clipping the sample width LW by 50% was observed to yield changes in vT and vG 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.Table 1
Overview of microstructure parameters. Intrinsic permeabilitiesk were calculated using single-phase flow simulations.
Microstructure A B C Pore-throat size distribution Mean valuemT [mm] 0.42 0.41 0.41 Standard deviationvT [mm] 0.14 0.15 0.21 Grain diameter distribution Mean valuemG [mm] 1.63 1.61 1.60 Standard deviationvG [mm] 0.20 0.34 0.44 Miscellaneous Porosityϕcuc [%] 45.5 46.0 44.3 Intrinsic permeabilitykcuc [10−9 m2] 4.41 4.55 4.55 Dimensions [mm2] 20 × 15 Grain count 88 85 86Varying the piston velocityUP, which also represents the prescribed specific discharge, we control the capillary number(8)Ca≔μnUPσwn.Primary drainage processes are studied forCa∈{10-2,10-3,10-4,10-5}. Fluid viscosities are varied to prescribe the viscosity ratio:(9)M≔μnμw.The simulated set of viscosity ratios isM∈{10-1,100,101} by considering the viscosity 2-tuples log10(μ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=σwncos⁡Θ 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 mT≈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 mT, time stepping constraints become increasingly restrictive.SPH particles are initially placed on the vertices of a Cartesian grid with grid spacingdx. Our choice of numerical resolution dx depends on mean pore-throat size mT as well as Ca and M. 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 dx=m¯T/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 dx=m¯T/14=29.5 μm. It was previously shown in [45, 50] that the incorporated choices of dx can reproduce local pressure, menisci, and velocity profiles with reasonable accuracy. The total number of SPH particles is approximately 6·105. Average computation time for a single time step was 0.6 s using a moderately optimized code running 6 threads in parallel on an Intel® Xeon® CPU E5-1650. The total number of time steps varied between 1.1·106 for small and 1.6·105 for large capillary numbers.Animated simulation results showing phase, pressure, and velocity magnitude distributions during primary drainage are available as supplementary information (MoviesS01–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 M>1, and viscous fingering at sufficiently large Ca and M<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 PC and linear growth of specific fluid interfacial area awn (Figure 2(b)) until breakthrough occurs. Macroscopic capillary pressure is computed as the difference in reservoir-averaged fluid pressures PC≡〈pn〉res-〈pw〉res whereas awn is the ratio of total fluid-fluid interfacial area to total unit cell volume. While constant PC 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].Figure 2
(a) NWP distributions (black) at breakthrough for microstructure C in(Ca,M) feature plane. Solid phase and WP are colored gray and white, respectively. (b) Normalized macroscopic capillary pressure PC(Sw)/PCref and specific fluid interfacial area awn(Sw) for capillary-dominated primary drainage in microstructure C at log10(Ca,M)=(-5,0). During the preentry regime, PC increases to the entry capillary pressure level while awn(Sw) decreases as the meniscus comes in contact with the solid phase. During the postentry regime, PC exhibits a constant mean value with fluctuations related to spontaneous Haines jumps.
(a) (b)In an attempt to highlight methodological limitations, we discuss the most crucial restrictions of the present simulations. (a) The present model is only applicable to ideal solid surfaces with absent chemical imperfections, surface roughness, or dust particles. In the presence of inhomogeneous solid surfaces a phenomenon referred to as contact line hysteresis has to be taken into account. (b) Given the computational costs associated with pore-scale resolved SPH simulations, we restrict ourselves to two-dimensional problems. The latter implies topological constraints related to connectivity in 2D microstructures as well as a lack of representing depth-related effects, for example, due to curvature of menisci in the out-of-plane direction. (c) When a gas bubble invades a capillary tube that is initially saturated with a WP of nonnegligible viscosity, a residual wetting film is hydrodynamically entrained between gas bubble and channel wall regardless of solid wettability characteristics [56]. The latter constitutes a special case of the Landau-Levich-Derjaguin problem [57]. When the thickness of the wetting film δ~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.


## 2.1. 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 Ωα∈{Ωn,Ωw} read(1)ϱ˙α=-ϱαdivuα,(2)ϱαu˙α=divTEα-gradpα,respectively. The Newtonian viscous extra stress tensor reads TEα=μα(graduα+gradTuα) and uα denotes local velocity. Dynamic viscosity μα is considered constant. A superimposed dot denotes the material time derivative. Fluid pressure pα is related to bulk mass density ϱα by a stiff barotropic equation of state pα=f(ϱα)+p0 such that density fluctuations are small implying quasi-incompressibility. A constant, positive backpressure p0 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(3)ϱu-uαβ·nαβαβ=0,where uαβ denotes interfacial velocity, imply kinematic coupling along the interfacial normal vector nαβ. Interfacial balances of linear momentum read(4)TE-pI·nαβαβ=-divαβΠαβ,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(5)TEn-TEw·nαβ+pw-pn+σwnκwnnwn=0,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(6)ϱfu˙f=divTEf-gradpf+∑αβdivΠαβδαβ·Ifs,where the solid surface identity tensor Ifs≔I-nfs⊗nfs is introduced to remove the otherwise unbalanced contact line stress σwnsin⁡Θ 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(7)miu˙i=∑jFijV-∑jFijP+∑αβFiαβfor a fluid particle Pi∈Ωf with lumped mass mi. Upon discretization, internal bulk forces due to viscous diffusion and pressure gradient result in discrete interaction forces (FijV and FijP) that act between Pi and its neighboring particles Pj.


## 2.2. Simulation Setup and Microstructures

The total computation domain (Figure1) 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 LW=20 mm and LH=15 mm in direction of the unit vectors e1 and e2, respectively. The microstructure exhibits LH-periodicity in direction of e2 such that periodic boundary conditions are applied with respect to e2. 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 UP=UPe1 of pistons that displace the reservoir fluids relative to Ωcuc. As a result of choosing UP constant, saturation rates remain constant during drainage. Simulations are halted as soon as the NWP penetrates into the WP reservoir, that is, at breakthrough.Figure 1
Schematic diagram of simulation setup at initial timet=0. Black bars represent displacing pistons.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 mT and vT for pore-throat size distribution and mG and vG for grain diameter distribution, respectively. Microstructures are constructed such that they only differ by the degree of heterogeneity as measured by vT and vG, while mean values mT and mG, intrinsic permeabilities, and porosities are in close agreement (Table 1). Mean pore-throat size mT serves as characteristic length scale for the reference capillary pressure PCref≔2cos⁡ΘσwnmT-1. Given its small size, the computational unit cell Ωcuc does not constitute a representative unit cell. However, clipping the sample width LW by 50% was observed to yield changes in vT and vG 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.Table 1
Overview of microstructure parameters. Intrinsic permeabilitiesk were calculated using single-phase flow simulations.
Microstructure A B C Pore-throat size distribution Mean valuemT [mm] 0.42 0.41 0.41 Standard deviationvT [mm] 0.14 0.15 0.21 Grain diameter distribution Mean valuemG [mm] 1.63 1.61 1.60 Standard deviationvG [mm] 0.20 0.34 0.44 Miscellaneous Porosityϕcuc [%] 45.5 46.0 44.3 Intrinsic permeabilitykcuc [10−9 m2] 4.41 4.55 4.55 Dimensions [mm2] 20 × 15 Grain count 88 85 86Varying the piston velocityUP, which also represents the prescribed specific discharge, we control the capillary number(8)Ca≔μnUPσwn.Primary drainage processes are studied forCa∈{10-2,10-3,10-4,10-5}. Fluid viscosities are varied to prescribe the viscosity ratio:(9)M≔μnμw.The simulated set of viscosity ratios isM∈{10-1,100,101} by considering the viscosity 2-tuples log10(μ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=σwncos⁡Θ 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 mT≈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 mT, time stepping constraints become increasingly restrictive.SPH particles are initially placed on the vertices of a Cartesian grid with grid spacingdx. Our choice of numerical resolution dx depends on mean pore-throat size mT as well as Ca and M. 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 dx=m¯T/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 dx=m¯T/14=29.5 μm. It was previously shown in [45, 50] that the incorporated choices of dx can reproduce local pressure, menisci, and velocity profiles with reasonable accuracy. The total number of SPH particles is approximately 6·105. Average computation time for a single time step was 0.6 s using a moderately optimized code running 6 threads in parallel on an Intel® Xeon® CPU E5-1650. The total number of time steps varied between 1.1·106 for small and 1.6·105 for large capillary numbers.Animated simulation results showing phase, pressure, and velocity magnitude distributions during primary drainage are available as supplementary information (MoviesS01–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 M>1, and viscous fingering at sufficiently large Ca and M<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 PC and linear growth of specific fluid interfacial area awn (Figure 2(b)) until breakthrough occurs. Macroscopic capillary pressure is computed as the difference in reservoir-averaged fluid pressures PC≡〈pn〉res-〈pw〉res whereas awn is the ratio of total fluid-fluid interfacial area to total unit cell volume. While constant PC 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].Figure 2
(a) NWP distributions (black) at breakthrough for microstructure C in(Ca,M) feature plane. Solid phase and WP are colored gray and white, respectively. (b) Normalized macroscopic capillary pressure PC(Sw)/PCref and specific fluid interfacial area awn(Sw) for capillary-dominated primary drainage in microstructure C at log10(Ca,M)=(-5,0). During the preentry regime, PC increases to the entry capillary pressure level while awn(Sw) decreases as the meniscus comes in contact with the solid phase. During the postentry regime, PC exhibits a constant mean value with fluctuations related to spontaneous Haines jumps.
(a) (b)In an attempt to highlight methodological limitations, we discuss the most crucial restrictions of the present simulations. (a) The present model is only applicable to ideal solid surfaces with absent chemical imperfections, surface roughness, or dust particles. In the presence of inhomogeneous solid surfaces a phenomenon referred to as contact line hysteresis has to be taken into account. (b) Given the computational costs associated with pore-scale resolved SPH simulations, we restrict ourselves to two-dimensional problems. The latter implies topological constraints related to connectivity in 2D microstructures as well as a lack of representing depth-related effects, for example, due to curvature of menisci in the out-of-plane direction. (c) When a gas bubble invades a capillary tube that is initially saturated with a WP of nonnegligible viscosity, a residual wetting film is hydrodynamically entrained between gas bubble and channel wall regardless of solid wettability characteristics [56]. The latter constitutes a special case of the Landau-Levich-Derjaguin problem [57]. When the thickness of the wetting film δ~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.


## 3. Results and Discussion

### 3.1. Viscous Fingering Regime: Hydrodynamic Wetting Films

Viscous fingering, i.e., the Saffman-Taylor instability [60], occurs when a less viscous fluid displaces a more viscous fluid (M<1) at sufficiently large capillary numbers. While invasion percolation is considered a suitable stochastic model for capillary fingering, viscous fingering is stochastically reminiscent of diffusion-limited aggregation (DLA) [61–63]. In this section, we empirically study the annular pore-scale flow profiles during viscous fingering and the impact of wetting films on the evolution of specific interfacial area and related dynamic effects. Hydrodynamic entrainment of wetting films at large 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 δ~mTCam2/3, where the microscopic capillary number Cam incorporates a characteristic microscopic velocity rather than UP. 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 (Figure3). Transition from spherical tip to wetting film is observed less obvious when the finger percolates through a pore-body (Figure 3, t0+0.6ms). 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, t0+0.9ms). 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, t0+0.9ms).Figure 3
Hydrodynamic wetting film entrainment in subdomain of microstructure A forlog10(Ca,M)=(-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.Following [66], we consider our choice of viscosity parameters for M=10-1 a rough representation of a carbon dioxide- (CO2-) water system in deep sedimentary formations. In particular, the modeled NWP has small but nonnegligible viscosity (μn=0.1 mPa s). Annular flow may thus constitute a relevant transport mechanism during supercritical CO2 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 atCa=10-2 and M=10-1, the rate of change of awn is observed to be significantly larger as compared to stable displacement and capillary fingering (Figure 4(a)). The magnitude of awn 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.Figure 4
(a) Evolution of specific fluid-fluid interfacial areaawn 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 SnB denotes NWP saturation at breakthrough.
(a) (b)In an attempt to quantify the dynamic role of wetting films, we measure the total forces the reservoir pistons have to do work against to displace the NWP into the pore-space. For viscous-dominated flow we expect the total resistance force against drainage to be additively composed of the total skin friction force due to viscous shear between fluid and solid phase, that is, Darcian drag, and the total viscous coupling force acting between WP and NWP, that is, the Yuster effect [67, 68]. Total skin friction and viscous coupling force are evaluated as(10)FfsV≔∮∂ΩsTEf·tfsda,FwnV≔∮ΓwnTEn·twnda,respectively, where tαβ denotes a vector tangent to the interface Γαβ. The ratio FwnV/(FfsV+FwnV) 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, M>100) viscous coupling forces are negligible (Figure 4(b)). The latter does not apply for viscous fingering where the relative dominance of FwnV increases monotonically with Sn and in a manner qualitatively similar to awn(Sn). At breakthrough, nearly one-sixth of total resistance is due to viscous coupling and this trend is expected to be even more pronounced for finer microstructures, larger domain sizes, and higher saturations. These results suggest that viscous coupling terms should be accounted for in macroscopic models for viscous fingering in porous media since otherwise lubrication effects are lumped into relative permeability functions, which, by definition, are related to momentum exchange between solid and fluid phases only.


### 3.2. Stable Displacement Regime: Capillary Dispersion Zone

Viscosity stabilizes interfacial perturbations when a more viscous fluid displaces a less viscous fluid (M>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 M>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–71].Disregarding boundary effects, saturation profilesSw(X1) that arise during stable displacement are sigmoidal in shape and, in a rough approximation, upper and lower asymptotes to Sw(X1) can be considered the saturation limits Swu=1 and Swl≈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.Figure 5
Evolving saturation profiles in flow direction for stable displacement in microstructures A (a) and C (b) atlog10(Ca,M)=(-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 ΔX1=LW/20. Red markers indicate residual saturation profile SwB(X1) at breakthrough.
(a) (b)We define the widthLC(t) of the CDZ as spatial difference between both points where saturation profiles Sw(X1,t) approach the asymptotes. While these asymptotes are observed to be barely sensitive to properties of the microstructures, the width LC(t), 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 LC(t) during drainage, hereafter denoted by L¯C. In order to avoid boundary artifacts, we consider the averaging window for L¯C 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 vG of grain diameter distribution (Figure 6(a)). In particular, L¯C≈14vGmG yields a reasonable fit to the present data. We emphasize that drainage rates, interfacial tensions, and viscosities, which are all expected to affect LC(t), are kept constant meanwhile. Moreover, microstructures A–C only differ by the degree of heterogeneity expressed in terms of vG while porosities, intrinsic permeabilities, and mean values mG and mT 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 (~vG).Figure 6
(a) Normalized average widthL¯C of CDZ during stable displacement as a function of standard deviation vG of grain diameter distribution. Blue error bars correspond to standard deviations of LC(t). (b) Normalized pressure profiles in flow direction for microstructure C at log10(Ca,M)=(-2,1). Solid lines represent average fluid pressure profiles p¯(X1,t) within vertical averaging slices of width ΔX1=LW/50. Dotted and dashed lines correspond to lower and upper bounds to fluid pressure distribution and envelop 87% of pressure data (shaded area).
(a) (b)Within the CDZ, fluid interfaces and hydraulically reservoir-disconnected WP clusters are subject to the interplay of viscous and capillary forces. This is discussed in terms of the temporal evolution of mean pressure profilesp¯(X1,t) (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 (Figure7). 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, t0). 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, t1). 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, t2). 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.Figure 7
Pore-scale interfacial dynamics during stable displacement atlog10(Ca,M)=(-2,1). Visual representation is equivalent to Figure 3. Spherical wetting cap formation due to menisci overlap near the saturation front ((a), t0). Fragmentation of a pendular bridge into multiple wetting caps within the CDZ ((b), t1). Wetting phase cluster mobilization by sequence of fragmentation and coalescence events within CDZ ((c), t2).
(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 χw of the total WP to emphasize topological implications of 2D simulations. In 2D, the Euler characteristic can be expressed as χw=Nw-Lw, where Nw denotes the number of isolated components and Lw 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 componentsNw 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(Sn), Nw(Sn) and Lw(Sn) with saturation is constant, that is, invariant with respect to the position of the compact saturation front. At initial time t0, WP saturates the entire pore-space such that Lw(t0) is equal to the genus of the pore-space and Nw(t0) is one. As saturation levels increase both Nw 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.Figure 8
(a) Number of isolated components (Nw) and Euler characteristic (χw) of WP for viscous fingering at log10(Ca,M)=(-2,-1) in microstructure B. (b) Euler characteristics of WP and reservoir-connected WP (χw,C) as well as number of isolated components for stable displacement at log10(Ca,M)=(-2,1) in microstructure B.
(a) (b)For viscous fingering, on the other hand, 2D simulations are not topologically representative of 3D processes. In 3D, dynamic wetting films are considered to extend over the surface of the entire pore-network and remain connected to the WP reservoir such thatχ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 Nw increases with the number of isolated films, the total number of holes Lw(Sn)=Lw(t0), 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 Nw 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 Nw. This implies that Lw 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 awn and viscous coupling FwnV as discussed in Section 3.1 are expected to apply regardless of above topological issues.


## 3.1. Viscous Fingering Regime: Hydrodynamic Wetting Films

Viscous fingering, i.e., the Saffman-Taylor instability [60], occurs when a less viscous fluid displaces a more viscous fluid (M<1) at sufficiently large capillary numbers. While invasion percolation is considered a suitable stochastic model for capillary fingering, viscous fingering is stochastically reminiscent of diffusion-limited aggregation (DLA) [61–63]. In this section, we empirically study the annular pore-scale flow profiles during viscous fingering and the impact of wetting films on the evolution of specific interfacial area and related dynamic effects. Hydrodynamic entrainment of wetting films at large 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 δ~mTCam2/3, where the microscopic capillary number Cam incorporates a characteristic microscopic velocity rather than UP. 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 (Figure3). Transition from spherical tip to wetting film is observed less obvious when the finger percolates through a pore-body (Figure 3, t0+0.6ms). 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, t0+0.9ms). 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, t0+0.9ms).Figure 3
Hydrodynamic wetting film entrainment in subdomain of microstructure A forlog10(Ca,M)=(-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.Following [66], we consider our choice of viscosity parameters for M=10-1 a rough representation of a carbon dioxide- (CO2-) water system in deep sedimentary formations. In particular, the modeled NWP has small but nonnegligible viscosity (μn=0.1 mPa s). Annular flow may thus constitute a relevant transport mechanism during supercritical CO2 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 atCa=10-2 and M=10-1, the rate of change of awn is observed to be significantly larger as compared to stable displacement and capillary fingering (Figure 4(a)). The magnitude of awn 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.Figure 4
(a) Evolution of specific fluid-fluid interfacial areaawn 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 SnB denotes NWP saturation at breakthrough.
(a) (b)In an attempt to quantify the dynamic role of wetting films, we measure the total forces the reservoir pistons have to do work against to displace the NWP into the pore-space. For viscous-dominated flow we expect the total resistance force against drainage to be additively composed of the total skin friction force due to viscous shear between fluid and solid phase, that is, Darcian drag, and the total viscous coupling force acting between WP and NWP, that is, the Yuster effect [67, 68]. Total skin friction and viscous coupling force are evaluated as(10)FfsV≔∮∂ΩsTEf·tfsda,FwnV≔∮ΓwnTEn·twnda,respectively, where tαβ denotes a vector tangent to the interface Γαβ. The ratio FwnV/(FfsV+FwnV) 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, M>100) viscous coupling forces are negligible (Figure 4(b)). The latter does not apply for viscous fingering where the relative dominance of FwnV increases monotonically with Sn and in a manner qualitatively similar to awn(Sn). At breakthrough, nearly one-sixth of total resistance is due to viscous coupling and this trend is expected to be even more pronounced for finer microstructures, larger domain sizes, and higher saturations. These results suggest that viscous coupling terms should be accounted for in macroscopic models for viscous fingering in porous media since otherwise lubrication effects are lumped into relative permeability functions, which, by definition, are related to momentum exchange between solid and fluid phases only.


## 3.2. Stable Displacement Regime: Capillary Dispersion Zone

Viscosity stabilizes interfacial perturbations when a more viscous fluid displaces a less viscous fluid (M>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 M>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–71].Disregarding boundary effects, saturation profilesSw(X1) that arise during stable displacement are sigmoidal in shape and, in a rough approximation, upper and lower asymptotes to Sw(X1) can be considered the saturation limits Swu=1 and Swl≈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.Figure 5
Evolving saturation profiles in flow direction for stable displacement in microstructures A (a) and C (b) atlog10(Ca,M)=(-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 ΔX1=LW/20. Red markers indicate residual saturation profile SwB(X1) at breakthrough.
(a) (b)We define the widthLC(t) of the CDZ as spatial difference between both points where saturation profiles Sw(X1,t) approach the asymptotes. While these asymptotes are observed to be barely sensitive to properties of the microstructures, the width LC(t), 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 LC(t) during drainage, hereafter denoted by L¯C. In order to avoid boundary artifacts, we consider the averaging window for L¯C 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 vG of grain diameter distribution (Figure 6(a)). In particular, L¯C≈14vGmG yields a reasonable fit to the present data. We emphasize that drainage rates, interfacial tensions, and viscosities, which are all expected to affect LC(t), are kept constant meanwhile. Moreover, microstructures A–C only differ by the degree of heterogeneity expressed in terms of vG while porosities, intrinsic permeabilities, and mean values mG and mT 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 (~vG).Figure 6
(a) Normalized average widthL¯C of CDZ during stable displacement as a function of standard deviation vG of grain diameter distribution. Blue error bars correspond to standard deviations of LC(t). (b) Normalized pressure profiles in flow direction for microstructure C at log10(Ca,M)=(-2,1). Solid lines represent average fluid pressure profiles p¯(X1,t) within vertical averaging slices of width ΔX1=LW/50. Dotted and dashed lines correspond to lower and upper bounds to fluid pressure distribution and envelop 87% of pressure data (shaded area).
(a) (b)Within the CDZ, fluid interfaces and hydraulically reservoir-disconnected WP clusters are subject to the interplay of viscous and capillary forces. This is discussed in terms of the temporal evolution of mean pressure profilesp¯(X1,t) (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 (Figure7). 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, t0). 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, t1). 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, t2). 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.Figure 7
Pore-scale interfacial dynamics during stable displacement atlog10(Ca,M)=(-2,1). Visual representation is equivalent to Figure 3. Spherical wetting cap formation due to menisci overlap near the saturation front ((a), t0). Fragmentation of a pendular bridge into multiple wetting caps within the CDZ ((b), t1). Wetting phase cluster mobilization by sequence of fragmentation and coalescence events within CDZ ((c), t2).
(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 χw of the total WP to emphasize topological implications of 2D simulations. In 2D, the Euler characteristic can be expressed as χw=Nw-Lw, where Nw denotes the number of isolated components and Lw 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 componentsNw 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(Sn), Nw(Sn) and Lw(Sn) with saturation is constant, that is, invariant with respect to the position of the compact saturation front. At initial time t0, WP saturates the entire pore-space such that Lw(t0) is equal to the genus of the pore-space and Nw(t0) is one. As saturation levels increase both Nw 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.Figure 8
(a) Number of isolated components (Nw) and Euler characteristic (χw) of WP for viscous fingering at log10(Ca,M)=(-2,-1) in microstructure B. (b) Euler characteristics of WP and reservoir-connected WP (χw,C) as well as number of isolated components for stable displacement at log10(Ca,M)=(-2,1) in microstructure B.
(a) (b)For viscous fingering, on the other hand, 2D simulations are not topologically representative of 3D processes. In 3D, dynamic wetting films are considered to extend over the surface of the entire pore-network and remain connected to the WP reservoir such thatχ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 Nw increases with the number of isolated films, the total number of holes Lw(Sn)=Lw(t0), 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 Nw 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 Nw. This implies that Lw 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 awn and viscous coupling FwnV as discussed in Section 3.1 are expected to apply regardless of above topological issues.


## 4. Summary and Open Questions

While pore-scale imaging methods, for example, fast X-ray tomography methods, to study the distribution and evolution of fluid interfaces for capillary-dominated flow become increasingly mature, their applicability to viscous-dominated flow is limited by temporal resolution. We consider hydrodynamic direct numerical simulations a suitable complementary approach for pore-scale modeling of viscous-dominated two-phase flow. Using a Smoothed-Particle Hydrodynamics model, we performed pore-scale simulations of primary drainage in partially wettable 2D porous media of particulate microstructure at largeCa. We have demonstrated the use of the mesh-free model to study hydrodynamic entrainment of thick wetting films during viscous fingering as well as the dispersion of stable saturation fronts.While the impact of wetting films on the dynamics of viscous fingering was outlined, the identification of capillary numbers where the onset of wetting film formation takes place requires further efforts. As the modeling of thin wetting films requires disjoining pressure forces to be accounted for, the present model was only applied to thick hydrodynamic wetting films well within the viscous fingering regime. While the relative dominance of interfacial viscous coupling forces is observed to monotonically increase during viscous fingering, their impact on relative permeability functions remains open. While mesh-free methods were shown to be a viable approach to study viscous fingering phenomena, significant computational resources and algorithmic optimizations are currently required to render 3D mesh-free simulations feasible.In contrast to viscous fingering, the 2D topological properties of the compact wetting phase during stable displacement were found qualitatively equivalent to those in 3D. The observation that the width of capillary dispersion zones increases with microstructure heterogeneity is thus considered to apply to 3D processes as well. The latter may thus motivate the regularization of Buckley-Leverett models to account for capillary dispersion. Mobilization, fragmentation and coalescence of disconnected wetting phase clusters were found as processes entirely local to the capillary dispersion zone. Hence, the dispersion zone decisively impacts residual wetting phase distribution in form of disconnected wetting caps and pendular bridges. This implies that the region of interest to study wetting phase entrapment during stable displacement using pore-scale imaging methods can be restricted to the capillary dispersion zone. For an accurate numerical prediction of threshold sizes for wetting phase cluster mobilization in reservoir systems, however, more realistic models that take into account contact angle hysteresis are considered necessary.



---

*Source: 8269645-2018-06-14.xml*