Direct Numerical Simulation of Pore-Scale Trapping Events During Capillary-Dominated Two-Phase Flow in Porous Media

This study focuses on direct numerical simulation of imbibition, displacement of the non-wetting phase by the wetting phase, through water-wet carbonate rocks. We simulate multiphase flow in a limestone and compare our results with high-resolution synchrotron X-ray images of displacement previously published in the literature by Singh et al. (Sci Rep 7:5192, 2017). We use the results to interpret the observed displacement events that cannot be described using conventional metrics such as pore-to-throat aspect ratio. We show that the complex geometry of porous media can dictate a curvature balance that prevents snap-off from happening in spite of favourable large aspect ratios. We also show that pinned fluid-fluid-solid contact lines can lead to snap-off of small ganglia on pore walls; we propose that this pinning is caused by sub-resolution roughness on scales of less than a micron. Our numerical results show that even in water-wet porous media, we need to allow pinned contacts in place to reproduce experimental results.


Introduction
Two-phase flow through geological porous media occurs during enhanced oil recovery (Lake 1989), water infiltration in soils (Hill and Parlange 1972), remediation of contaminated aquifers (Mulligan et al. 2001) and geological storage of carbon dioxide (Marini 2006;Yamasaki 2003). These flows are dictated by the behaviour at the pore scale (Blunt 2017). The dynamics of pore-scale events such as snap-off can significantly influence the macro-scale efficiency of these processes. Capillary trapping of the non-wetting phase helps to maximize the trapping efficiency for carbon dioxide storage, while for oil recovery and subsurface contaminant clean-up processes, minimal trapping is desirable.

3
The prediction of immiscible two-phase flow displacement events, however, can be a significant challenge because of the complicated interaction between the fluid-fluid interface and the solid structure. Moreover, different degrees of affinity for fluids to cover the rock surface, i.e., different wettability, surface roughness and heterogeneity can influence capillary forces that can play a significant role in pore-scale displacement processes (Butt 2008;Singh et al. 2017). A combination of corner flow, piston-like displacement and snapoff events occur during capillary driven imbibition (Blunt 2017).
Snap-off is one of the mechanisms of trapping non-wetting fluid in pores (larger spaces) connected by throats (smaller restrictions connecting pores). During imbibition, wetting fluid may accumulate in roughness and corners of the pore space (Yu and Wardlaw 1986b). When the wetting phase loses contact with the solid, the interface between the two fluid phases may become unstable resulting in rupture of the non-wetting phase breaking its connection between adjacent pores as the wetting phase fills the throat between them. This event is called snap-off and is a function of interrelated parameters such as capillary pressure, wetting layer thickness, pore geometry and wettability (Yu and Wardlaw 1986a, b). Snap-off can play a major role in trapping continuous CO 2 in saline aquifers or waterwet reservoir rocks during water-flooding (Pickell et al. 1966;Roof 1970;Lenormand et al. 1983;Altman et al. 2014;Blunt 2017). Therefore, understanding and predicting the dynamics of snap-off events can be of great importance to determine the amount and location of the trapped non-wetting phase.
Several experimental and theoretical studies have investigated the mechanism of snapoff. Yu and Wardlaw (1986b) performed quasi-static two-phase flow experiments through a single pore-throat pair micro-model. They concluded that the capillary pressure, disjoining pressure, thin-film thickness, and wettability (contact angle) control snap-off mechanisms in a two-phase flow single-throat system. They quantified critical aspect ratios (pore-tothroat inscribed radius ratio at which snap-off occurs) for different ranges of advancing contact angles. In another study, Yu and Wardlaw (1986a) provided information on topological aspects of different snap-off scenarios by performing two-phase flow experiments on two-dimensional transparent micro-models with various interconnected pores. They found that, in addition to wettability, topology plays a major role in determining the location of snap-off (throat, pore, or pore-throat-junction). Beresnev et al. (2009) formulated a geometric criterion for snap-off in capillary-driven imbibition by analysing capillary-pressure distributions in capillary channels with sinusoidal geometry. The main limitation of this analysis, however, is that it is based on static or quasi-static criteria and, therefore, they do not take into account the dynamics of two-phase flow.
These dynamics are typically characterized by the capillary number, the ratio of the viscous force to the capillary force. A number of experimental and computational studies have been performed to investigate the impact of two-phase flow dynamics characterized by capillary number on snap-off (Payatakes 1982;Mogensen and Stenby 1998;Raeini et al. 2014a;Deng et al. 2015). However, due to experimental and computational limitations, these studies only considered simple geometries that do not represent the complexity of geological porous media.
Recent advances in micro-CT scanning can provide high-resolution three-dimensional images of pore-scale displacement mechanisms in geological porous media (Wildenschild et al. 2002;Blunt et al. 2013;Andrew et al. 2014;Arns et al. 2007;Pak et al. 2015;Singh et al. 2016). In addition, high-resolution fast synchrotron imaging has made it possible to visualize and track the time evolution of interfaces during snap-off events in drainage and imbibition processes (Coles et al. 1998;Armstrong et al. 2014;Singh et al. 2017;Andrew et al. 2015). This information can be employed to analyse the dynamics of snap-off within the complex geometry of permeable media, and also to quantify the related parameters such as the critical capillary pressure and curvature, which are a function of geometry and wettability. Moreover, the availability of micro-CT imaging has given new impetus to pore-scale modelling. This type of experiments can serve as a basis to validate pore-scale models including direct numerical and pore-network models (Alpak et al. 2019;Bultreys et al. 2018). On the other hand, direct numerical simulations can complement the experimental analyses by simulating and probing phenomena that classical rules fail to describe (Rabbani et al. 2017).
The goal of this study is to further understand the impact of realistic pore structures on the occurrence of snap-off in a capillary dominated flow regime, where capillary forces dominate over viscous forces. We accomplish this goal by predicting snap-off events within the porespace of a Ketton limestone rock sample via two-phase direct numerical simulation (DNS), and performing one-to-one comparisons against the experimental data (Singh et al. 2017(Singh et al. , 2018. This is in contrast to previous studies which have primarily treated snap-off events as a static (or quasi-static) phenomena, or have modelled snap-off evolutions within synthetic single-pore geometries.
This work is organized as follows. First, we present the direct numerical simulation method used to model immiscible two-phase flow directly on micro-CT images (Shams et al. 2018a;Raeini et al. 2012). Then, we briefly describe the two-phase experimental data that we use as the base for our numerical simulations. Results are then presented in two parts. First, we demonstrate the fidelity of our two-phase flow DNS model in predicting snap-off events by simulating the imbibition process on cropped image samples where snap-off has been observed in the experiments. Then, we employ DNS to simulate and describe the experimental snap-off events that do not follow classic snap-off criteria established in the literature. We simulate a case for which snap-off does not happen in spite of a large aspect ratio, which is favourable for snap-off, and also two snap-off cases where the interface ruptures unexpectedly and consequently the non-wetting phase gets trapped. We evaluate the influence of pore geometry on the interface curvature that may impose a piston-like displacement where snap-off is expected. Moreover, we investigate whether and how the roughness of solid boundaries can promote the pinning of the interface causing the trapping of small ganglia of the non-wetting phase.

Numerical Method
The two fluids are modelled as one single-fluid continuum, and the following single set of the Navier-Stokes equations is numerically solved to describe isothermal, incompressible, immiscible two-phase flow through porous media, where is the velocity vector, D Dt ( ) = t ( ) + ∇ ⋅ ( ) is the inertial term, p d is the dynamic pressure, = (∇ + (∇ ) T ) is the viscous stress tensor, is the local density, is the local dynamic viscosity, and c is the capillary force computed as a volume force based on a contour-level surface force model (see Shams et al. (2018a) for more details); (1) p c , is a potential field called the microscopic capillary pressure, which is defined as (Raeini et al. 2012) Total (physical) pressure, p, then can be computed as (see Raeini et al. (2012) for more details), The interface reconstruction based surface force model in conjunction with some numerical filtering techniques (Shams et al. 2018a;Raeini et al. 2012), is employed to efficiently reduce spurious currents allowing us to simulate two-phase flow at relatively low capillary numbers, Ca = 10 −5 to Ca = 10 −6 ; see Raeini et al. (2012Raeini et al. ( , 2014a and Shams et al. (2018a, b) for more details.
The fluid-fluid interface is implicitly captured as a discontinuity in fluid properties using an indicator function, , where the subscripts 1 and 2 denote the first and second fluid, respectively. In the volumeof-fluid approach, the indicator function physically represents the volume fraction of one of the two fluids that varies smoothly between 0 and 1. The interface can be defined as the surface where the value of is equal to 0.5.
The kinematics of the interface is modelled using the following advection equation, The finite volume method, with a collocated arrangement of variables at centres of discrete non-overlapping control-volumes, is used to discretize the governing equations; see Deshpande et al. (2012) and references therein for more details. The flow algorithm is implemented within OpenFOAM, an open-source computational fluid dynamics package (Weller et al. 1998;OpenFOAM. 2016). On solid boundaries, the wetting effects are taken into account using a prescribed contact angle, (Shams et al. 2018a). The contact angle is an input parameter set equal to the equilibrium contact angle.

Two-Phase Flow Experiment
The two-phase flow experiment that we simulate in this study, was performed by Singh et al. (2017). They employed fast synchrotron X-ray micro-tomography to obtain threedimensional high-resolution (with a voxel size of 3.28 m) dynamic images of twophase fluid flow during brine injection (imbibition) through a water-wet carbonate rock (a 3.8-mm-diameter and 10-mm-long sample) at ambient temperature (20 • C). The flow regime was capillary dominated with a capillary number of approximately 1.26 × 10 −9 . This capillary number is a global (or bulk) capillary number which uses the total Darcy velocity of the invading fluid (wetting phase). At the pore level, however, the local flow velocity of the invading phase within a single pore can be orders of magnitude higher than the total Darcy velocity of the bulk porous medium (Deng et al. 2015;Berg et al. 2013). Therefore, one can assume that, for this experiment, the local (or pore) capillary number, which is determined by the flow velocity within single pores, is approximately of the order of O (10 −7 ) or smaller. Direct numerical simulation of the imbibition process at the same size-and time-scale as in the experiment, however, can be computationally prohibitive. To address the size-scale issue, we perform simulations on smaller cropped sub-images (maximum 300 3 voxels) around the regions where snap-off occurs. The time-scale for snap-off events in the experiment, swelling of the wetting layer prior to the instantaneous layer rupture, was reported to be of the order of tens of minutes. However, within a feasible computational time, we are able to simulate a dynamic process of the order of only a few seconds. To resolve this issue, we increase the swelling speed of the wetting layers in our simulations by injecting the wetting phase about two to four orders of magnitude faster than in the experiment. This results in higher local capillary numbers in our simulations, Ca = 10 −5 − 10 −6 , compared to global/Darcy capillary number of Ca = 10 −9 in the experiment. However, since the simulation capillary number is still very low, the flow regime remains highly capillary dominated, and therefore, simulated dynamic processes inside pores closely follow the ones in the experiment.
We perform direct numerical simulations on six micro-CT images, which are selected and cropped from the whole micro-CT image obtained from the experiment; see Table 1.
The flow domain is discretized to an unstructured mesh using SnappyHexMesh, a mesh generator from OpenFOAM (2016). It uses a uniform hexahedral cell background mesh and snaps the boundary points to the solid boundary leading to split-hexahedral cells near the grain walls. We use a constant uniform contact angle of 45 • , consistent with the direct in situ measurements inside the rock. The interfacial tension is constant and taken to be = 0.04 N/m. Since the injection velocity is very low, and flow regime is capillary dominated, the impact of viscous and inertial forces are insignificant compared to interfacial forces. Therefore, we set viscosity and density ratios equal to one in all the simulations.
A no-slip boundary condition is applied to the grain walls. We consider a fixed-flux boundary condition, with a uniform value of zero, for the dynamic pressure, p d , at the inlet. The wetting phase is injected at the inlet using a relaxed-fixed injection rate boundary condition designed to gradually increase the injection rate to the target rate, here Ca ≈ 10 −5 , while mimicking a zero-gradient boundary condition for the velocity. The velocity boundary condition at the outlet is zero-gradient for normal flow (positive fluxes) and zero fixedvalue for backward flow (negative fluxes). The dynamic pressure boundary condition acts to balance the velocity boundary condition, i.e., zero fixed-value for normal flow (positive flux) and zero-gradient for backward flow (negative flux) at the outlet.
The sub-domain cropped for the two-phase simulation is in capillary equilibrium with the surrounding domain during the capillary dominated flow experiments. To impose this 1 3 equilibrium condition in our simulations, a zero-gradient boundary condition is used for the microscopic capillary pressure, p c , at all boundaries. For the indicator function, we consider a fixed-value boundary condition, = 0 , at the inlet and a zero-boundary condition for the other boundaries. The phase distribution imaged at the end of the drainage experiment is used as the initial condition for the indicator function in our simulations.

Results and Discussion
The test cases 1-3 serve as validation to check the ability of our numerical method to predict snap-off events. The location and size of the trapped phase for each case is compared to the experimental results. Figure 1 shows these three test cases at their initial state; the micro-CT sub-images obtained from the flow experiment are shown in the first row while their equivalent numerical flow domain and the initial state of the non-wetting phase are shown in the second row. Figure 2 shows a comparison of the phase distributions between the experiment and simulations after snap-off occurs. These visualizations indicate a good qualitative agreement between numerical and experimental results. Table 2 reports the relative error in the size of the trapped ganglia, where v exp and v num are the volume of the trapped ganglion computed from the experimental data and the numerical simulation, respectively. The computational results agree with the experimental data with a maximum relative error of 14%. The discrepancy can be due to the uncertainty in some parameters such as contact-angle, roughness and hysteresis effects, as well as the difference between the capillary number of the experiment and simulations. There are also uncertainties in the segmentation of the experimental images, which leads to an uncertainty in the measured ganglion size; however, since each ganglion Fig. 1 Visualizations of the non-wetting phase distribution during the imbibition experiment used as the initial condition for the indicator function in the validation test cases. The wetting phase and solid are transparent with light shading of the pore space occupies around 100,000 voxels, this error is likely to be small. These simulations exhibit the ability of the numerical model to predict snap-off events. In both the experiments and the simulations, flow is capillary controlled. We assume that the computational sub-domains are sufficiently large that boundary effects do not significantly affect the results. Hence, we are able to make quantitative comparisons between the imaged and predicted fluid distributions. However, the results cannot be used to predict the associated kinetics/time-scale of displacement events such as the swelling rate of wetting-phase layers and snap-off time since our simulations are run at a different capillary number to the experiment. Figure 3 shows a one-to-one comparison between the observed and simulated displacement of the non-wetting phase for case 1. It is clear that even though the model can predict the sequence of displacement, there is a significant time-scale discrepancy between numerical and experimental results. The time-scale of the simulated displacement events, t s , is almost four orders of magnitude smaller than observed time-scale in the experiment, t, reported in Singh et al. (2017).
This discrepancy, as mentioned in Sect. 2, is due to the fact that in our simulations the invaded-phase is injected with a higher velocity, Ca ≈ 10 −5 , compared to the experiment, Ca ≈ 10 −9 . In addition, since our simulations are performed on small sub-samples (limited to a few pores) without enforcing kinematic continuity of the velocity between the subdomain and the whole flow domain, simulation results cannot be relied upon to accurately predict the associated kinetics of the flow.  In the following, we use the computational model as a predictive tool to simulate the last three test cases (cases 4-6) for which the experimental results are not consistent with the expected behaviour regarding snap-off. Figure 4 shows the flow domain and initial phase distributions for the three test cases at the beginning of the imbibition simulations.
The geometric properties of pores and throats can be obtained from a pore-throat network extracted from the image . In case 4, the pore-throat aspect ratios for the throats connected to the pore, are approximately 4.4 and 3.3 . Therefore, snap-off is expected to occur in both throats. However, the experimental results show that the non-wetting fluid passes through the throats in a piston-like fashion. As shown in Fig. 5 (case 4), the numerical simulation predicts the same behaviour just as in the flow experiment.  For further investigation, we have post-processed the simulation data to compute average capillary pressures within three small separate control volumes around the throats and the pore. The capillary pressure results presented in Fig. 6, show that the higher capillary pressure in the pore maintains the stability of piston-like displacement of the interface within the pore. Moreover, we observe that the microscopic capillary pressure dictates the order in which the wetting phase fills the pore and throats; the pore is filled first-see the sudden drop in capillary pressure-followed by the two throats and not the other way round as would be expected based on the aspect ratio. The radius of the flow path is just one of the several interrelated parameters such as corner angle, contact angle, and roughness that determine the interface curvature (or micro-scale capillary pressure). Therefore, for complex media such as porous rocks, pore-throat aspect ratio can not always be relied upon to predict snap-off. As it has been shown in this test case, the ultimate factor in determining the flow pattern is the interface curvature, which itself is dependent on dynamic and geometric factors (Rabbani et al. 2018).
As shown in Fig. 5 for cases 5 and 6, the numerical simulations are not consistent with the experimental results. In these cases, the experimental results show that a small fraction of the non-wetting phase has been left behind sticking to the pore wall, while the numerical simulation predicts a complete pore filling. This kind of trapping mechanism seems to be triggered by pinning of the moving interface at the contact line due to sub-micron resolution roughness at the solid boundary .
We have examined this hypothesis by setting the contact angle around the small region where trapping occurs to a different value higher than 45 • , as shown with green patches in Fig. 7a. This allows us to introduce the pinning mechanism through a phenomenological model, namely an advancing contact angle boundary condition, representing the impact of the sub-micron resolution roughness on the contactline movement. We set the advancing contact angle, A , to 70 • . As shown in Fig. 7c, the simulations with the new contact angle condition have predicted the trapping of  The wetting phase and solid are transparent with light shading of the pore space the non-wetting phase similar to the experimental results. The value of the advancing contact angle is considered as a tuning parameter that controls the size of the trapped ganglia.
In Fig. 8, we visualize the simulation of non-wetting phase displacement process at different times for case 6, with the interface pinning mechanism in place. One can clearly observe that, due to the interface pinning, a small portion of the non-wetting phase remains stuck to the pore wall. The volume of non-wetting phase shrinks with time as the invading wetting phase continues to pass over it (see Fig. 8c-d). Interfacial transverse curvatures developed along the body of the stuck non-wetting phase start to grow and destabilize the interface within a relatively short span of time, which eventually leads to snap-off (see Fig. 8e-g).
For any arbitrary control volume within the flow system, the rate of change of dissipated energy due to viscous forces, Ė , is equal to the negative sum of the rate of change of energies introduced into the system by inertial forces, dynamic pressure gradient, and capillary forces, where Ė I , Ė pd , and Ė pc are the rate of change of energies introduced into the system by inertial forces, dynamic pressure gradient, and capillary forces, respectively.
The mechanical energy-balanced equation for a control volume, V, can be obtained by taking the inner product of the momentum equation, Eq.
(2), with the velocity vector field, , and integrating it over the control volume (Raeini et al. 2014b), These relations can be used to quantify the contribution of the pore-scale forces during displacement in terms of the rate of change of dissipated/introduced energy by each force within the system (Raeini et al. 2014b).
For simulation case 6, in addition to the visualization of the non-wetting phase displacement presented in Fig. 8, we have numerically computed the energy rate terms, Eq. (10). The results, presented in Fig. 9, confirm that the flow is mainly controlled by a competition between viscous and capillary forces. Capillary forces are the major source of energy change during the process, and the energy introduced by the capillary force field is balanced by viscous forces,

Fig. 9
Rate of change of energy due to different forces, Eq. 10, during two-phase flow displacement for simulation case 6. Energy can be introduced into the flow system by inertial forces, Ė I , dynamic pressuregradient forces, Ė pd , and capillary forces, Ė pc , or can get dissipated from the flow system by viscous forces, Ė We can also observe that, apart from the two spikes in the rate of change of energy, the flow system is stable, and phases are in full equilibrium with almost no change in the rate of change of energy during the displacement. The first spike in the rate of change of energy is due to the instability introduced by the rearrangement of fluid phases at the start of the displacement process, see Fig. 8a-c. The second spike signals the onset of the instability introduced by the ever-increasing interfacial transverse curvatures along the body of the stuck non-wetting phase, which leads to the break-up of the fluid-fluid interface, see Fig. 8e-g.
There are two parameters that influence the occurrence and the size of trapped ganglia in cases 5 and 6: (a) the size of the boundary patch for which the contact angle is different, and (b) the advancing contact angle, A , assigned to the patch. We have performed a sensitivity analysis to assess the impact of these two parameters on the trapping mechanism. For this purpose, four patches with different surface areas are selected around the region that we have observed trapping in the experiment in case 6: patches 1-4 with surface areas of, respectively, 0.006 mm 2 (the green patch shown in Fig. 7a), 0.008 mm 2 , 0.01 mm 2 and 0.016 mm 2 . For each patch, the two-phase simulations are performed within a range of A from 45 • to 95 • . Figure 10 depicts the results in terms of the relative error in the size of the trapped ganglion, Eq. (7), as the function of the advancing contact angle for patches of different size. Trapping does not happen for any of the patches when A is less than 70 • . It can be observed that by increasing A , interface pinning becomes more dominant leading to a larger trapped ganglion. It is also clear that by increasing the size of the patch, a larger surface area contributes to the pinning, which results in a larger ganglion.
The wide error range, from 2% to over 250%, indicates that the accurate prediction of the size of the trapped ganglion is dependent on the size of the patch and its contact angle. This implies that to increase the ability of the numerical model to predict such a trapping (11) ∫ V ⋅ (∇ ⋅ )dV ≈ ∫ V ⋅ (∇p c − c )dV.

Fig. 10
The relative error in the size of the trapped ganglion, E r , Eq. (7), for case 6 as a function of the advancing contact angle, A , and the size of the patch on which A is applied mechanism, sub-resolution roughness and its related dynamics need to be incorporated in the model.

Conclusions
Three-dimensional computational fluid dynamics simulations were performed directly on images obtained from an X-ray micro-CT two-phase flow experiment. The ability of the computational model to predict snap-off and trapping events was validated by comparing against experimental observations. The simulation results revealed that the interface configuration of two fluids can be affected by the complex geometry of the flow domain, and the resultant interface curvature is the major factor that controls the pore filling pattern. This implies that the traditional concept of pore-to-throat aspect ratio cannot always predict snap-off events accurately, and it may explain some of the discrepancies in pore occupancy analysis between experimental data and pore network results (Bultreys et al. 2018;Raeini et al. 2019). The results obtained in this study suggest the importance of roughness on the dynamics of moving interfaces and fluid flow patterns. Sub-resolution roughness can pin the contact-line which leads to snap-off and trapping of the non-wetting phase in pores, a process that is not captured in traditional pore network models, for instance. Further work could consider mixed-wet media, where the trapping mechanism due to a pinned contactline is likely to be much more significant and may indeed dominate the behaviour (Rücker et al. 2019(Rücker et al. , 2020Lin et al. 2019).