Dynamic pore-scale modeling of residual fluid configurations in disordered porous media

. Fluid trapping in porous media is important in many subsurface flow processes such as enhanced oil recovery and geological sequestration of carbon dioxide. To achieve optimal performance in such applications, a fundamental understanding of residual trapping mechanisms at the pore scale is necessary. In this work, we present a computational study of fluid trapping behaviors in natural porous media under different flow regimes by employing a dynamic pore-network modeling approach. The model incorporates many advanced features that have not been collectively used in previous dynamic platforms. For instance, it rigorously solves for fluid pressure fields from two-phase mass balance equations in each pore element, incorporates a detailed description of pore-scale fluid displacement dynamics of piston-like advance, snap-off, and pore-body filling, and explicitly accounts for flow through wetting layers forming in corners and rough surfaces of pore spaces. Moreover, we extend the ability of our model by including contact angle hysteresis, which is often neglected in existing dynamic models. A heavily-parallelized implementation of this platform is further advanced to achieve an efficient computational performance. We first conduct primary drainage and imbibition simulations in pore networks representing Bentheimer and Berea sandstones. We show that the predicted two-phase relative permeability curves agree well with their experimental counterparts reported in the literature. Afterwards, the validated model is used to systematically probe fluid trapping behaviors in a core-sized pore network that is constructed from high-resolution micro-computed tomography images of a Berea sandstone core sample. The effects of dynamic flow conditions and fluid properties on core-scale two-phase displacement pattern, residual-fluid configuration, and residual oil saturations are examined in detail. Fluid trapping properties such as maximum and average residual-fluid cluster size and capillary-controlled invasion selectivity at the pore scale are analyzed under both capillary-and viscous-dominated flow regimes.


Introduction
Residual trapping of the non-wetting fluid in porous media is important in many science and engineering applications such as enhanced oil recovery (EOR), geological sequestration of carbon dioxide, and contaminant removal from groundwater [1][2][3].To achieve the optimal performance in these applications, it is crucial to accurately model the flow behaviors and predict the residual fluid configurations.However, multiphase flow systems are exceedingly complex due to, for instance, the complicated pore space geometries and topologies, fluid-fluid capillary interaction, and wettability effects [4][5][6][7][8][9][10][11].Among many other efforts to better understand the residual trapping behaviors, extensive numerical simulation studies have been conducted at multiple scales.In the conventional continuum models, nevertheless, most of the complex pore-scale processes are significantly simplified and taken into consideration using a series of input functions, e.g., relative permeability, with little physical basis.Yet these properties are expensive and time-consuming to obtain experimentally.Therefore, physics-based models that can predict fluid trapping properties and bridge the gap between the pore and continuum scales flow processes are desperately needed.Over the years, with the advances in imaging techniques, e.g., micro-computed tomography (micro-CT), different pore-scale simulation techniques have been proposed and investigated to fulfill different purposes.In general, porescale models can be categorized into direct pore-level models and pore-scale network models.In the first group, computational fluid dynamics (CFD) methods, such as levelset [12], lattice Boltzmann [13,14], and volume-of-fluid [15], are utilized to solve the governing flow equations in the pore space obtained through the imaging techniques.The direct models preserve the detailed pore space geometries.However, the complex pore space geometries of naturallyoccurring rocks make their computations expensive.The complicated pore-scale fluid dynamics of two-phase system carry the modeling works into a further level of complexity.In the second group of pore-scale models, fluid flow and solute transport are solved in a network of inter-connected pores and throats with idealized geometries.Through this approximation, the fluid interfacial movements can be described by a set of displacement rules with closed-form mathematical expressions to substantially reduce the computational costs.In 1956, Fatt [16] made the pioneering work of modeling drainage process in porous media by using a two-dimensional (2D) network of cylindrical tubes.Since then, tremendous progress has been made in developing different network models, which mainly include quasi-static and dynamic models.In quasi-static models, different fluids are assumed to be at capillary equilibrium and solely the capillary pressure determines the fluid displacement sequence [17][18][19].Therefore, the quasi-static models are only applicable to slow flow processes where the capillary number is approximately in the order of 10 −7 or lower.To account for the viscous effect, the dynamic models solve the fluid pressure field by writing mass conservation equations for all volume-occupying phases in every pore spaces.In this regard, the computational demand of the dynamic models is orders of magnitude higher than that of the quasi-static simulation.Koplik et al. [20], Dias and Payatakes [21], and Lenormand et al. [22] are among the first to develop the dynamic algorithm for two-phase flow in porous media.They wrote the Washburn equation at the fluid-fluid invasion front to account for the capillary and viscous forces.In these early studies, many pore-scale fluid dynamics such as fluid flow through the wetting layers are neglected, even though they are crucial in maintaining the connectivity of the wetting phase and affecting the residual trapping behaviors of the non-wetting fluid during two-phase flow [23,24].Later, Singh and Mohanty [24] proposed a rule-based scheme to implement layer flow for primary drainage displacement.They calculated the wetting phase pressure in layers using a spatial extrapolation approximation and updated the wetting layer volume following a heuristic rule that the amount of wetting phase removed from a layer was proportional to its local capillary pressure drop.Nguyen et al. [25] studied the effect of flow rate on imbibition displacement properties assuming that no viscous pressure drop was associated with the non-wetting phase and the pressure of the wetting phase was constant behind the displacement front.With a similar approach, Wang et al. [26] simulated two-phase imbibition in a regular 2D network consisting of 800 pores and 1540 throats.Thompson [27] introduced a fully dynamic model, in which the local capillary pressure is approximated as a function of the local water saturation in regular cubic pore spaces.To extend the Thompson's approach, Sheng and Thompson [28] conducted steady-state flow simulations to predict two-phase relative permeability; Joekar-Niasar and Hassanizadeh [29] investigated the non-equilibrium capillary pressure during two-phase flow.More recently, Li et al. [30] studied the effect of viscosity ratio on two-phase imbibition inside a regular 2D pore network consisting of 100 × 50 pores and a regular three-dimensional (3D) pore network with 20 × 20 × 20 pores.Boujelben et al. [31] presented a new dynamic model to study low-salinity water injection, where the wetting layer was treated as volumeless.Qin et al. [32] extended the model of Joekar-Niasar and Hassanizadeh [29] and studied air-water drainage flow with phase change in two pore networks with 1,163 and 10,619 pores, respectively.Yelkhovsky and Pinczewski [33] modeled two-phase capillary rise in both round and polygonal capillary tubes by solving an averaged Navier-Stokes equation.While the geometries of the pore spaces of naturally occurring porous media can be complex [34,35], to the best of our knowledge, existing fully dynamic models always assume the pore elements have square cross sections.Moreover, factors such as neglecting wettability hysteresis, ignoring throat volumes, and neglecting snap-off displacements in pore bodies further limit the predictability of the previous developed models.This study aims at developing a new, generalized dynamic pore-network modeling platform to study fluid trapping behaviors in natural porous media under different flow regimes.We start with an introduction of the representative pore networks used in this work.This is followed by a detailed description of the new dynamic pore-scale modeling algorithm and its parallel implementation.Then, quantitative and qualitative validations of the modeling platform are presented.Afterwards, the model is used to examine the effects of dynamic flow conditions on two-phase displacement patterns, residual oil cluster configurations, and final residual oil saturation during different imbibition displacements.In addition, the impacts of fluid properties including viscosity ratio and wettability are discussed in detail.Finally, the main findings of this work are summarized in the conclusion section.
2 Pore network constructions In this work, we utilize four networks for different purposes.For the purpose of validation, we first generate a 2D lattice pore network that is statistically equivalent to the micromodel reported in the experimental study of Lenormand et al. [22].This 2D pore network consists of 100 × 100 pores with cross-sectional shape of square.Then, we utilize two pore networks being representative of the Berea and Bentheimer sandstones to predict their two-phase relative permeabilities.The Berea network (see Fig. 1a), having a dimension of 3 × 3 × 3 mm 3 , is constructed from a process-based technique [36] and has been successfully tested in previous modeling studies [37][38][39].The Bentheimer network (see Fig. 1b), having a length of 5.00 mm and a cross-sectional diameter of 4.97 mm, is extracted from the high-resolution (scanning resolution 2.6 µm) micro-CT images of a Bentheimer sandstone rock sample [40].The Bentheimer rock sample has a measured porosity of 23.42% and an absolute permeability of 2568 mD.To investigate two-phase flow behaviors and residual fluid configurations in disordered porous media, we further extract a longer, cylindrical pore network directly from highresolution (scanning resolution 2.49 µm) micro-CT images of a Berea sandstone sample with a measured porosity of 21.69% and an absolute permeability of 623 mD [40].The network (see Fig. 1c) has a length of 10.374 mm and a crosssectional diameter of 3.276 mm.Hereinafter, we will refer to this network as the large Berea network to distinguish it from the smaller network constructed using the process-based technique.Pores and throats of the three pore networks can have different cross-sectional shapes of square, equilateral and scalene triangles.
3 Model description

Model assumptions
Fluids are assumed to be Newtonian, incompressible, and immiscible.The flow process is isothermal and there is no interaction between the fluids and the rock material.Also, clay content is assumed to be fully saturated with immobile water.It is worth noting that in our model, both pore and throat elements have their own volume and conductivity.This is distinct from the existing dynamic models, where zero volume or zero conductance are often imposed on either throat elements or pore bodies to simplify the computations.Hereinafter, we will use the terms "water" and "oil" to represent the wetting and non-wetting phases, respectively, and "pore elements" to denote both pore bodies and pore throats.

Pore-level physics of two-phase flow
Based on the experimental observations of two-phase flow in porous media [41][42][43], three fundamental pore-level displacing mechanisms are included in this dynamic model: piston-like, cooperative pore-body filling, and snap-off.Their detailed descriptions are presented in our previous work [44].However, it is worth mentioning that snap-off corresponds to the displacement of the non-wetting phase in the center of a pore element by arc menisci (AMs) formed in the corners and on the rough surfaces of the pore space.When local capillary pressure decreases, the relevant wetting corners gradually swell.Eventually, at a critical threshold capillary pressure, oil at the center is snapped off and the element is spontaneous filled by water.This displacement mechanism significantly contributes to oil entrapment in porous media [45,46].

Phase hydraulic conductance
The hydraulic conductance of a pore element filled with a single phase is calculated based on its geometries [18,47]: where R, A, and G are the element's inscribed radius, crosssectional area, and shape factor, respectively; and µ is the fluid viscosity.After a local drainage, it is impossible for oil to completely drain water out of an angular pore element as a certain amount of water remains in the corners.In such fluid configurations, the flow area open to the corner water and the center oil can be computed as: (2) where α is the corner half angle of the pore element; θow is the contact angle of the oil-water-solid system; r is the curvature radius of the AMs; σow is the oil-water interfacial tension; pc is the local capillary pressure; and o and w are the phase indexes denoting oil and water phases.The center oil conductance is then evaluated using Ao and Eq. 1, while the calculation of the corner water conductance can be found in [44,48,49].
One should note that at pore scale, the oil conductivity of an element could instantly change from a finite value to zero immediately after the oil phase is snapped off or displaced by water.

Governing equations
The volume conservation equation of each fluid phase is applied at every pore and throat: where Vi and Ni are the respective total volume and coordination number of element i; Sα and pα are the saturation and pressure of phase α.The equivalent phase conductance , through two neighboring elements, i and j, is considered as the harmonic mean of the phase conductance of i and j.
In addition, we know that for any element containing two phases, the pressures of oil and water are related by the local capillary pressure, as: and the summation of phase saturations in i follows: Combining Eqs.5-7, we obtain the following pressure equation, as: We choose pw as the principal unknown because water is hydraulically connected across the entire pore network.

Local displacement
At pore scale, the occurrence of the local displacement is controlled by the threshold capillary pressure,  .In this work, the threshold capillary pressure for piston-like displacement is calculated using the Mayer-Stowe-Princen (MSP) theory [50,51].Interested readers are referred to [37,52] to find the detailed computations of the  for pistonlike, cooperative pore-body filling, and snap-off displacements.However, it is worth mentioning that for a given pore element, the threshold capillary pressure for piston-like displacement is always larger than that of the snap-off event.Therefore, piston-like displacement is more favored when the fluid configuration is allowed.However, snap-off does not require the presence of a main terminal meniscus (MTM) to occur, meaning that it can take place anywhere in the porous media as long as the fluids are hydraulically connected.In this work, for a partially invaded element, we assume its local capillary pressure equals to the threshold capillary pressure of the local displacement (i.e., local drainage or imbibition) until the element is fully invaded.By the end of a local imbibition, oil is fully displaced and the local capillary pressure disappears.In contrast, after a local drainage, water layers form in the corners of the angular pore space.In this case, we define a critical saturation corresponding to the point when the entire central pore space is occupied by oil under the threshold drainage capillary pressure of  , , .Denoted as  , , , this critical saturation can be obtained using Eqs.
2-4 and  , , .Snap-off handling.Immediately after a water element is invaded by oil, we compute its threshold capillary pressure for snap-off,  , , , and its critical saturation  , , corresponding to  , , using Eqs.2-4.In our simulations, as fluid saturations are updated after every time step, a snap-off checking process is performed for all elements containing two phases; if the updated water saturation in an element i exceeds  , , , spontaneous snapoff will occur and the oil phase will break into n parts, where n is the coordination number of i.Note that during this process, new MTMs will form, resulting in an instant increase in the local capillary pressure.This, in turn, accelerates the retreat of oil from i to its neighbors [53].Accordingly, we update the local capillary pressure in element i based on the curvature of the newly formed MTMs: where Ri and  , are the inscribed radius and advancing contact angle of i.

Contact angle hysteresis and corner interface handling
We assign each pore element a receding ( ) and advancing ( ) contact angle for the local drainage and imbibition displacement, respectively.Moreover, within each element, different parts of the solid surface can exhibit different wettabilities.In specific, each element is initially full of water and is considered as strongly water-wet.After it is invaded by oil, the water-filled corners will remain strongly water wet, while the rest of the solid surfaces contacted by oil will become less water-wet (i.e., have higher contact angles).
When a new AM forms, its initial location can be uniquely determined by the threshold capillary pressure, contact angle, and the corner half angle of the pore element through a geometrical analysis [37].As fluid flow continues, the configuration of the AM is then subject to the dynamic change of the local capillary pressure.It can stay pinned with a fixed apex-meniscus distance b as the contact angle, , hinges between  and  or it can move at a constant  or  .The detailed calculation of the corner interface configuration can be found in [44].It is evident that such wettability hysteresis behavior results in a non-smooth local capillary pressure-water saturation function.

Local capillary pressure
We note that for a pore element containing two phases, the fluid interface configuration and the phase saturation can be uniquely determined under a given local capillary pressure, as illustrated in Appendix C in [37].In addition, at a new time step n+1, the local capillary pressure,  , , will be either in the range of [ , ,  , ] (or [ , , , ,  , ]) for a local drainage (or imbibition) displacement, where  , is known from the last time step n and  , is a local maximum.In view of these, we propose a bi-sectional iterative method to find the  , in any pore element using the newly updated local saturation,  , .The detailed calculation procedure is illustrated in Fig. 2. It is important to note that this iterative numerical scheme allows our dynamic pore-network modeling platform to work with pore elements with irregular and arbitrary cross-sectional shapes.

Boundary and initial conditions
At the inlet of the pore network, a Neumann boundary condition is employed to accommodate a constant fluid injection rate.Specifically, the following two equations are added to the linear system of the pressure equations: where Nin is the total number of inlet-connected throats; and Qα,in and pα,in are the fluid injection rate and injection pressure, respectively.At the outlet of the pore network, a constant production pressure pw,out is enforced for the water phase.When oil is also present at the outlet, we adopt the assumption proposed by Joekar-Niasar and Hassanizadeh [29] and set the local capillary pressure gradient in the outlet-connected throats as zero, namely, = 0.A no-flow condition is applied at the side boundaries of the network.
In our simulations, the pore networks are initially full of water and are gradually drained by oil during the primary drainage displacement.After a steady-state condition (or a targeted initial oil saturation) is reached, the simulation proceeds to the imbibition displacement, during which water is injected at the inlet of the pore network under a constant volumetric flow rate.

Numerical implementation
At a new time step n + 1, we first update the fluid pressure field implicitly by listing the pressure equation (Eq.8) for all connected pores and throats, as: A linear system of pressure equations is then obtained and solved with an algebraic multigrid preconditioned conjugated gradient method using PETSc parallel solver package [54].
Then, we explicitly update the fluid saturations in pores and throats by discretizing the volume balance equation (Eq.5), as: To keep the numerical stability, the time step size (∆t) is chosen in a way that at most one pore element can be fully invaded at every time step.Moreover, after every certain number of invasions, we adjust the pore-scale capillary pressures in accordance with the macroscopic capillary pressure distribution to avoid possible numerical oscillations in local pc.The inherent assumption here is that in a porous medium, the variation of the capillary pressure along the principal flow direction is smooth and continuous.To this end, the macroscopic pc distribution along the length of the pore network is first calculated as the averages of the porescale capillary pressures over the cross-sections of the network.Then, the capillary pressures in different pore elements are updated according to the element's distance from the inlet of the network and the macroscopic pc distribution.

Domain decomposition and parallelization
For a fully dynamic pore network model, the computation of phase pressure updating, local capillary pressure iteration, and fluid interface tracking can become significant, particularly for multiphase displacements in a large pore network.In this work, the modeling platform is designed to run on massively parallel computing architectures.To this end, the entire large Berea network is divided into a number of sub-domains.The MVAPICH2 library, a high performance message passing interface (MPI) implementation, is used to communicate information between different computing processors.Please refer to our previous works [44] for more details of parallel computing.

Results and discussions
In section 4.1, we first validate our dynamic pore network model.Afterwards, we conduct a series of two-phase displacement simulations to examine the effects of flow rate, viscosity ratio (M), and wettability on dynamic displacement behaviors and residual fluid configurations.Specifically, we conduct five imbibition simulations (simulations 1-5) under elevated water injection rates.Then, another two imbibition simulations (simulations 6 and 7) are performed with enhanced water viscosity to study the viscosity ratio effect.
We define the M as the ratio of the viscosity of the invading phase to that of the defending phase.Lastly, four more imbibition simulations (simulations 8-11) are designed with various oil-water contact angles to study the wettability effect.In all simulations, the imbibition displacements are initiated at an initial oil saturation of Soi = 0.9.Properties of the fluids used in these simulations are listed in Table 1.The detailed simulation conditions and residual oil saturations (Sor) are summarized in Table 2.

Model validation
To validate our model, we first carry out a series of primary drainage simulations in the 2D network over a wide range of capillary numbers and viscosity ratios, following the experimental conditions reported in Lenormand et al. [22].3. Three important two-phase displacement configurations of stable displacement (Fig. 3a), viscous fingering (Fig. 3b), and capillary figuring (Fig. 3c) are clearly identified and they agree well with the experimental observations presented in [22].In the stable displacement scenario with a viscosity ratio of M >> 1 and a large capillary number, the invasion front between the fluids is flat and oil invades nearly all of the pore spaces behind of it.In this case, the pore-scale displacements exclusively take place at the invasion front because the viscous pressure drop associated with the invading phase dominates the displacing process (Fig. 3a).As the viscosity ratio becomes much smaller than one, there is almost no viscous resistance to the flow of the invading phase, while the viscous defending phase consumes nearly all the pressure drop.In consequence, as the invading oil creates a small number of fingers with relatively high conductivities at the inlet, a few of them grow rapidly toward the outlet as the invasions continues, resulting in a quick breakthrough and leaving most of the network uninvaded (Fig. 3b).With the decrease in Nc, the viscous force becomes negligible in both fluids and the pore-scale displacement sequence is largely controlled by the capillary force.In this case, the invading phase selectively takes paths with more favorable capillary pressures and develops capillary fingering, as shown in Fig. 3c.Note that the configurations of the viscous fingering and capillary fingering have distinct differences.The former invasion configuration primarily orient towards the outlet, while the latter one grows anomalously inside the network.This is also in line with the experimental observations reported in [22].We further conduct primary drainage and imbibition simulations in the small Berea and the Bentheimer networks to predict their two-phase relative permeability curves at a low flow rate.During the simulations, the phase relative permeability is computed as: where the subscripts s and t represent the single-and twophase flow condition, respectively, and ∆pα is the [33]pressure drop of phase α across the network.A receding contact angle of 0° is used for the primary drainage displacements, whereas advancing contact angles in the ranges of 60−75° and 30−50° are assigned for the Berea and the Bentheimer networks, respectively, following the works of [37] and [18].In Figs. 4 and 5, the computed oil and water relative permeabilities versus water saturation for the small Berea and the Bentheimer networks are compared quantitatively against their corresponding measured data collected from the work of Oak [55] and Oren et al. [18].As shown in these two figures, our predicted relative permeability curves agree closely with the experimental data for both primary drainage and imbibition displacements.Moreover, the dynamic model accurately predicts the residual oil saturations at the end of the imbibition displacements in both sandstones.We also note that at the early stage of the primary drainage, the predicted water relative permeability underestimates the experimental data.This discrepancy is caused by the fact that phase saturations along the network are not uniform before oil breakthrough.Overall, the simulation results illustrate that our dynamic pore network model is capable of capturing twophase flow properties under a wide range of flow conditions.

Fluid displacement patterns
In this section, simulations 1-5 are conducted with elevated water injection flow rates to achieve different numbers, Nc, ranging from 2.0 × 10 -7 to 2.0 × 10 -3 .Fig. 6 shows the fluid pressure distribution along the Berea network at early stage of the imbibition process in simulations 1-5.
As shown in Fig. 6, the oil pressure is higher than the water pressure in all five simulations due to the positive capillary pressure under the water-wet condition.However, at low capillary numbers (Nc ≤ 2.0 × 10 -5 ), the capillary pressure difference (∆pcap) between oil and water is much higher than the viscous pressure drop (∆pvis).In this case, pore-scale displacement is driven by the capillary force, meaning that water prefers to fill the elements with more favorable threshold capillary pressures.Moreover, the low flow rate allows sufficient time for water to penetrate through the wetting layers and do displacement in elements that are ahead of the invasion front.The low capillary number and the hydraulically connected water channels jointly lead to a capillary-dominated imbibition pattern: local invasion can take place throughout the network and they proceeds largely in the order controlled by the size and geometry of the pores, regardless of their locations along the flow direction.As a result, water saturation profile shows a uniform increase across the network in Fig. 7b.In contrast, at high capillary numbers of Nc ≥ 2.0 × 10 -4 , the viscous pressure drop surpasses the capillary pressure difference.In particular, simulation 5 has a pronounced viscous pressure drop of ∆pvis = 172 kPa compared to a capillary pressure difference about ∆pcap = 10 kPa, as depicted in Fig. 6e.This significant viscous pressure gradient favors frontal displacement (piston-like and pore-body filling) close to the water invasion front.Meanwhile, the rapid water injection rate leaves limited time for water flowing through the wetting layers with restricted flow conductance.This further suppresses water invasions ahead of the invasion front.Therefore, a piston-like invasion front develops and gradually moves towards the outlet of the network as the injection continues, as shown in Fig. 7c.

Trapped oil clusters and residual oil saturations
During imbibition in a water-wet rock, bypassing and snapoff displacement can lead to entrapment of oil and formation of oil clusters with complex structures and a wide range of sizes [23,56].This section presents characterizations of the trapped oil clusters, including the volumetric cluster size distribution and the final residual oil saturation, under different capillary numbers.At the end of each imbibition displacement, we calculate a total amount of the trapped oil clusters ( , ) and a cumulative volumetric distribution to the  , from clusters that are smaller than a volume size of  , .Fig. 8 plots the computed cumulative volumetric distribution of trapped oil clusters after the imbibition processes in simulations 1-5.Fig. 9 summarizes the final residual oil saturation of simulations 1-5.Under the capillary-dominated flow conditions, the injected water has sufficient time to flow through the wetting layers and to swell in the corners of the angular porous space across the network, causing snap-offs in small pore elements.This results in a large amount of hydraulically disconnected oil regions.From Fig. 8, we find that the contributions of oil clusters with sizes up to 0.1 mm 3 are similar in simulations 1-3.Accordingly, Fig. 9 shows that up to a critical capillary number about 2 × 10 -5 , there is a minor drop in the Sor with the increase of the capillary number.As Nc increases, the viscous pressure gradient is sufficient to restrain snap-off throughout the imbibition process and favors more piston-like advances, which produce many small-to-medium sized blobs in the vicinity of the invasion front.As a result, in simulations 4 and 5 conducted under the viscous-dominated flow regime, the overall oil cluster sizes significantly decrease (see Fig. 8) and the final residual oil saturation drops sharply (see Fig. 9).The same tendency was also described in experimental studies [57,58].Note that all five simulations create many blobs with volumes in the range of 1.0 × 10 -9 ∼ 1.0 × 10 -8 mm 3 , which is even smaller than the order of element size.The formation of the sub-pore-sized blobs during imbibition was also observed in previous experimental work [58,59].Such isolated sub-pore-sized or pore-sized blobs creates local obstacles for water flow through the narrow pore or throat spaces across the network.We further plot the maximum trapped oil cluster size,  , , and average trapped oil cluster size,  , , of simulations 1-5 in Fig. 10.It is found that simulation 1 has a maximum trapped oil cluster volume of  , = 1.022 mm 3 , which is the largest among the five simulations.As the capillary number increases, both the maximum and the average trapped oil cluster size decreases as the high viscous pressure drop breaks the larger oil clusters.However, it is interesting to note that the variations of the maximum and the average trapped oil cluster sizes with the capillary number are different from the Sor-Nc curve shown in Fig. 9; as capillary number increases, both  , and  , show a continuous discernible decrease, although the drop in Sor is minor until the capillary number reaches about 2 × 10 -5 .This suggests that even under the capillary-dominated flow regime, the configurations of the residual oil clusters experience evident changes as the viscous pressure gradually builds up, although the magnitude of the viscous force is insufficient to push the oil clusters out of the porous media.Previous studies [58,60] have suggested that the driving force needed to mobilize a trapped oil cluster is inversely proportional to the cluster length.In this respect, the knowledge of the cluster size distribution of the residual oil is an important aid for optimizing the design of EOR processes.

Fluid property effect
In this section, we first investigate the effects of viscosity ratio on the residual oil saturation by performing simulations 6 and 7.Both simulations are conducted with the same fluid injection rate as that of simulation 3.However, the viscosity of water is increased to10.10 and 101.0 cP in simulations 6 and 7, respectively.In simulations 8-11, different oil-water contact angles are conducted to examine the wettability effect.In these four simulations, a low capillary number of 2.0 × 10 -7 is used to keep them under the capillary-dominated flow regime.The detailed fluid properties, simulation conditions, and the final Sor are tabulated in Tables 1 and 2.

Higher capillary number
The cumulative volumetric contribution of trapped oil clusters, and the maximum and the average trapped oil cluster sizes for simulations 6 and 7 are included in Figs. 8 and 10.As the viscosity ratio increases, more contributions to the final Sor are from the small-to-medium-sized oil clusters and both the maximum and the average oil cluster sizes are decreased significantly.A closer examination of the spatial distribution map of the trapped oil clusters along the length of the network in simulations 3, 6, and 7 is presented in Fig. 11, where each dot represents the volume-weighted centroid location of an oil cluster.We find that the oil globules distribute uniformly across the network and the majority of the globules' size is in the range of δ = 10 -6 ∼ 10 -3 mm 3 .As the viscosity ratio becomes higher, there is an evident narrowing in the range of δ because of the magnified viscous pressure gradient of water.Notably, the population of the subpore-and pore-sized oil globules reduces as δ becomes narrower.This can be attributed to the fact that partial invasions are more likely to be completed as the invading water becomes more viscous.The above trend also agrees with that observed from simulations 3-5 as we increase the water injection flow rate.Nevertheless, even under a same capillary number, simulations 6 and 7 exhibit lower final Sor and smaller  , and  , than simulations 4 and 5 (see Figs. 9 and 10).This is because the higher water-oil viscosity ratio further confines the occurrence of water invasion within the vicinity of the invasion front and thereby, creates a more stable displacement that leaves less residual oil by the end of the imbibition.

Wettability effect
The residual oil saturations of simulations 8-11 are summarized in Fig. 9.As indicated in Fig. 9, when porous medium is strongly water-wet ( ≤ 45°), there is a clear reduction in the residual oil saturation as the contact angle increases.When the contact angles are lower than 15° (simulation 8), a large number of snap-off events take place across the network under the capillary-dominated flow condition.Spontaneous snap-offs and water bypassing [23,61] jointly leave a significant amount of residual oil.As the contact angle increases, it becomes more difficult for snapoff invasion to occur.This is particularly the case for weakly water-wet porous media (simulations 10 and 11), where the majority of oil entrapment solely results from water bypassing.Consequently, much less residual oil is trapped in these two simulations.Similar results have also been reported in the experimental studies of [62,63].To further explore the rules of wettability on the residual trapping behavior, we categorize the overall pore elements into three groups based on their inscribed radius (R) and analyze the water saturation of each group.More specifically, water saturations of small (R < 10 μm), medium (10 ≤ R ≤ 45 μm), and large (R > 45 μm) pore elements are calculated at the end of drainage and imbibition displacements.The results are plotted and compared in Fig. 12.In general, the magnitude of the threshold capillary pressure is lower for larger pore elements assuming uniform contact angles and interfacial tension.Therefore, it is easier for the non-wetting phase to invade into large pores and throats.Accordingly, in Fig. 12, we find that nearly all the mediumto-large elements are invaded by oil while the water saturation is still as high as 77% in small elements after primary drainage.On the other hand, in the following imbibition process, the wetting water phase preferably fill the small elements due to their higher threshold capillary pressures.As shown in Fig. 12, water occupies almost all of the small elements, a fair amount of the medium-sized elements, and a limited amount of the large elements at the end of the imbibition in these five simulations.This selective pore-scale displacement leads to capillary fingering (see section 4.1), resulting in water bypassing and oil entrapment.From Fig. 12, we find that after imbibition displacement in water-wet porous media, the water saturations are always high in small elements.However, as the wettability condition is altered towards neutral wetness and the absolute difference in threshold invasion capillary pressures of pore elements across the pore network shrinks, we see a clear increase in the water saturation (i.e., reduced residual oil trapping) in medium-sized pore elements.In other words, the suppression of the capillary invasion selectivity is more pronounced in the medium-sized pores than in the small elements.For the large pores, however, improved oil recovery is only observed when the porous media have large contact angles of  = 60−75°.Note that only a limited number of piston-like invasions can occur in large pores under the strongly water-wet conditions.There is even a slight drop in the water saturation as the contact angle increases from 0−15° to 30−45°, which can be caused by the changes in the water-corner configurations.In summary, the residual oil saturation becomes lower as the fluids-rock system become more neutral wet.The additional oil recovery is mainly attributed to restrained snap-off invasions and suppressed capillary invasion selectivity, particularly in medium-to-large-sized elements.

Conclusions
In this work, we developed a new dynamic pore network model that is designed for simulating two-phase flow in disordered porous media.The model has several advantages in comparison with the existing dynamic models: it iteratively solves for local capillary pressure in pore elements with irregular cross-sectional shapes, explicitly accounts for contact angle hysteresis, and rigorously considers pore-scale two-phase flow dynamics including frontal displacement, snap-off, and flow through wetting layers in both pore bodies and throat elements.A heavily-parallelized computing scheme of this model was implemented to achieve computational efficiency.We performed a series of two-phase primary drainage and imbibition simulations in different pore networks representing two-dimensional (2D) micro-model, threedimensional (3D) Bentheimer and Berea sandstones.We illustrate that our model is capable of capturing fluid invasion behaviors of stable displacement, viscous fingering, and capillary fingering over a wide range of fluid-fluid viscosity ratios and capillary numbers.The predicted primary drainage and imbibition relative permeability curves of Bentheimer and Berea sandstones show satisfactory matches with their experimental counterparts.The validated model is then applied to systematically study the effects of dynamic flow conditions and fluid properties on two-phase displacement patterns, residual oil configurations, and final residual oil saturations.A series of dynamic imbibition simulations are conducted in a core-sized Berea sandstone pore network that is constructed from highresolution micro-computed tomography (micro-CT) images of a Berea rock sample.Results show that after capillarydominated imbibition displacement, a large amount of oil is trapped in the form of large-sized oil clusters because of the abundant snap-offs and bypassing.As the water injection flow rate increases up to a critical capillary number of about Nc = 2 × 10 -5 , there is a minor change in the final residual oil saturation.However, both the maximum and the average residual-oil cluster sizes reduce considerably.With a further increase in the water injection flow rate, a clear piston-like invasion front is formed at the core scale.Under such viscousdominated flow regime, the significant viscous pressure drop favors frontal displacement in the vicinity of the invasion front, resulting in less amount of large-sized residual oil clusters and a much lower residual oil saturation.Finally, the effects of fluid-fluid viscosity ratio and wettability on residual trapping are investigated.It is found that the overall residual oil cluster sizes are further decreased by enhancing the water-to-oil viscosity ratio, owing to a strengthened frontal displacement.The amount of residual trapping can also be reduced by altering the wettability condition toward neutral wetness, where the pore-scale capillary invasion selectivity is suppressed, particularly in medium-to-largesized pore elements.In summary, our parallelized dynamic pore network model provides novel insights into residual trapping behaviors in disordered porous media under a wide range of flow conditions and fluid properties.This is vital for optimizing subsurface applications such as enhanced oil recovery, carbon sequestration, and soil remediation.
We gratefully acknowledge the financial support of Thermo Fisher Scientific and the School of Energy Resources at the University of Wyoming.

Figure 1 .
Figure 1.Visualization of (a) the small Berea pore network, (b) the Bentheimer pore network, and (c) the large Berea pore network.Orange and gray colors represent the pores and throats, respectively.

Figure 2 .
Figure 2. Flow chart of the iterative local capillary pressure finding algorithm.

Figure 4 .
Figure 4. Predicted relative permeabilities (kr) versus water saturation (Sw) compared against experimental data collected from the literature [55] for (a) primary drainage and (b) imbibition displacements.Both simulations are conducted in the small Berea network.

Figure 5 .
Figure 5. Predicted relative permeabilities (kr) versus water saturation (Sw) compared against experimental data from the literature [18] for the (a) primary drainage and (b) imbibition displacements.Both simulations are conducted with the Bentheimer network.

Figure 7 .
Figure 7. Snapshots of fluid distributions in the large Berea network (a) after the primary drainage process and (b, c) during the imbibition processes with the capillary number of Nc = 2.0 × 10 −7 (simulation 1) and Nc = 2.0 × 10 −3 (simulation 5), respectively.Oil elements are depicted in red and water elements are depicted in blue.Partially invaded elements are shown in orange.

Figure 8 .
Figure 8. Cumulative volumetric distribution of trapped oil blobs to the overall trapped oil volume after the imbibition displacements in the large Berea pore network in simulations 1-5.

Figure 10 .
Figure 10.Calculated maximum and average trapped oil cluster sizes after imbibition in the large Berea network in simulations 1-7.

Figure 11 .
Figure 11.(a) Distribution map of the trapped oil clusters along the length of the large Berea network in simulations 3 (viscosity ratio of M = 1.22), 6 (M = 12.2), and 7 (M = 122).Each dot represents the location of the volume-weighted centroid of an oil cluster.The three vertical lines from left to right show the average volumes of the pore throats, pore elements, and pore bodies, respectively.Stacked histograms of the size distribution of the trapped oil clusters are shown in (b).

Figure 12 .
Figure 12.Comparison of the calculated water saturations in various size ranges of elements at the end of the primary drainage and imbibition displacements in simulations 1 and 8-10, all of which are conducted under a same capillary number of Nc = 2.0 × 10 -7 .

Table 1 .
Properties of the fluids used in simulations 1-11.