Pore-scale direct numerical simulation of Haines jumps in a porous media model

TL;DR

This study uses pore-scale direct numerical simulations to analyze fluid displacement and Haines jumps in a porous medium, revealing the influence of contact angle and capillary number on displacement patterns and jump dynamics.

Contribution

It demonstrates that contact angle and capillary number jointly influence pore-scale displacement patterns and Haines jump characteristics, expanding understanding beyond traditional capillary number effects.

Findings

Displacement patterns depend on contact angle and capillary number.

Haines jumps involve sudden interface movements with specific time and length scales.

Contact angle significantly affects Haines jump behavior in simple pore geometries.

Abstract

Direct numerical simulations are presented for a porous media model consisting of two immiscible fluids, an invading and defending phase, in a two-dimensional micro-geometry filled with randomly sized and randomly distributed cylinders. First, interface instability and penetration modes are studied when varying the wetting features of a single pore in the porous medium. It is found that the displacement patterns not only change with the capillary number, as previously observed, but also are a function of the contact angle, even for a viscosity ratio of unity. This is an important conclusion suggesting that capillary number and viscosity ratio alone cannot completely describe the pore-scale displacement. Second, rapid pore-scale displacement is considered, where the displacements are accompanied by sudden interface jumps from one site to another, known as Haines jumps. The characteristic time and length scales of a Haines jump are examined to better understand the transient dynamics of the jump. We then focus on analyzing the Haines jump in a simple pore configuration where cylinders of equal size are placed at the vertices of equilateral triangles. We use this geometry to provide more insight into the effect of the contact angle at which the Haines jump is predicted.

Equations

$$\text{partialderivative}\rho \mathbf{u}t+\text{divergence}\rho \mathbf{uu}=-\text{gradient}p+\text{divergence}\mu \text{quantity}(\text{gradient}\mathbf{u}+\text{gradient}{\mathbf{u}}^{\text{T}})+{\mathbf{f}}_{\sigma }$$

$$\text{partialderivative}\rho t+\text{divergence}\rho \mathbf{u}=0$$

$$\text{partialderivative}\gamma t+\text{divergence}\mathbf{u}\gamma =0$$

$${\mathbf{f}}_{\sigma }=\sigma \kappa \text{gradient}\gamma$$

$$\kappa =\text{divergence}\vec{\text{n}}$$

$$\vec{\text{n}}=-\frac{\text{gradient}\tilde{\gamma }}{\Vert \text{gradient}\tilde{\gamma }\Vert }$$

$$\text{divergence}\mathbf{u}=0$$

$$\Re \text{quantity}(\text{partialderivative}{\lambda }_{\rho }\overrightarrow{{\mathbf{u}}^{\ast }}{t}^{\ast }+{\text{gradient}}^{\ast }\cdot {\lambda }_{\rho }{\mathbf{u}}^{\ast }{\mathbf{u}}^{\ast })=-{\text{gradient}}^{\ast }{p}^{\ast }+{\text{gradient}}^{\ast }\cdot {\lambda }_{\mu }\text{quantity}({\text{gradient}}^{\ast }{\mathbf{u}}^{\ast }+{{\text{gradient}}^{\ast }{\mathbf{u}}^{\ast }}^{\text{T}})+\frac{1}{\mathrm{Ca}}{\kappa }^{\ast }{\text{gradient}}^{\ast }\gamma$$

$$\Re =\frac{\rho \text{absolutevalue}\mathbf{U}\overline{D}}{\mu }$$

$$\mathrm{Ca}=\frac{\mu \text{absolutevalue}\mathbf{U}}{\sigma }$$

Full text

1 Introduction

The flow of immiscible fluids in porous media has applications in subsurface water flows, tar sands oil production, and $\text{CO}_{2}$ sequestration, to name a few (Perazzo et al., 2018). At the micro-scale, where the capillary forces are dominant, the wetting effects can play a significant role in the displacement dynamics. For example, when a fluid is pushed into a porous medium, displacing an immiscible fluid, the distribution of fluids depends on wetting properties of the medium. In the pioneering work of Lenormand, the displacement patterns of a wetting fluid, which is known as the imbibition scenario, and a non-wetting fluid, known as the drainage scenario, are studied, showing that in slow displacement drainage, fingering like patterns emerge (Lenormand, 1990). This condition can lead to the invading phase not being able to completely drain the other (defending) phase, resulting in some of the displaced phase becoming trapped in the porous medium (Blunt and Scher, 1995; Singh et al., 2016, 2019). The trapping can be important as it affects the saturation of the invading liquid which limits the extent to which the defending fluid can be displaced from the medium (Anbari et al., 2018; Zacharoudiou et al., 2018).

The study of the influence of the contact angle on different displacement patterns is important, as it is also one of the key factors in defining the flow behavior, and therefore the resulting fluid distribution within a porous media (Zhao et al., 2016; Singh et al., 2019; Primkulov et. a, 2019). In the pioneering work of Cieplak and Robbins (Cieplak and Robbins, 1988, 1990), a two-dimensional quasistatic numerical model was developed to describe wetting conditions. However, the approach in (Cieplak and Robbins, 1988, 1990) simplifies the interface representation as circular arcs connecting the solid geometries, which in their work are two-dimensional array of disks with random radii. Recent particle-based simulations account for wall wettability in two-dimensional (Holtzman and Segre, 2015) and quasi-two-dimensional (Jung et al., 2016) models. However, particle-based simulations present difficulties in enforcing the incompressibility in very viscous flows. Nonetheless, all the studies above consistently show the transition from stable front displacement at low contact angles to capillary fingering at high contact angles.

Pore-scale models have increasingly become a viable tool for predicting two-phase flows through porous media, mainly due to recent improved developments; see (Zhao et al., 2019) for a comparison of a variety of pore-scale models. Direct numerical simulation techniques based on solving the full Navier-Stokes equations provide a robust and high accuracy solution for pore-scale analysis of multiphase flows in porous media. However, a great challenge in pore-scale direct numerical simulations is the inclusion of wetting effects into the numerical model through the contact angle boundary condition. There exist direct numerical studies on pore-scale modeling of wettability effects, such as the phase-field simulation in (Basirat et al., 2017) and the lattice Boltzmann simulation in (Liu et al., 2014). Here we present the results of a direct numerical method that is based on an interface capturing approach. While the majority of previous (Eulerian) interface capturing methods include wetting on body-fitted meshes (see e.g. (Bussmann et al., 1999; Afkhami et al., 2009) and for a review (Sui et al., 2014)), less progress has been made to develop numerical methods that handle arbitrary geometries. To that end, one attractive approach is the Immersed Boundary Method (IBM) (Peskin, 1972; Mittal and Iaccarino, 2005), because of its ability to model complex geometries, and its potential for modeling systems of solid bodies with arbitrary relative motion (O’Brien and Bussmann, 2020). Here we use an accurate and robust numerical model developed in (O’Brien and Bussmann, 2018) and later extended in O’Brien and Bussmann (2020) for combining wetting dynamics, in a Volume-of-Fluid (VOF) framework, and the IBM.

The result we present here is first of its kind in that it uses a combined VOF/IBM pore-scale direct numerical simulations to quantify the effects of the capillarity, characterized by the capillary number, $\mathrm{Ca}$, and wetting, characterized by the contact angle, $\theta$, on the displacement phenomena in a porous media model. We examine how the wetting controls displacement patterns, and show the crossover from a stable propagating front to a fractal pattern, namely a fingering instability. Our work extends on the literature by identifying the displacement regimes as functions of independently varying $\mathrm{Ca}$ and $\theta$. We illustrate the interface patterns as we change the contact angle, and analyze the instabilities in terms of the growth mechanisms and the rapid movement of interfaces between pores. Based on our results, we describe three main patterns and the crossover from fingering to stable flow as a function of $\mathrm{Ca}$ and $\theta$, at a viscosity ratio of $1$. We show that decreasing the contact angle results in a change of the flow behavior from strong fingering to stabilized front propagation, and that the contact angle effect diminishes as $\mathrm{Ca}$ increases. The threshold for this transition is related to the contact angle amongst other parameters. For wetting invading fluids, the front propagates rather smoothly, draining all of the displaced phase, while a non-wetting invading fluid can lead to residual fluid. Our simulations reveal how post-fingering production flux is different above a critical contact angle from the one below the critical contact angle where the front is stable, for sufficiently small $\mathrm{Ca}$, and that the production flux is similar when $\mathrm{Ca}$ is sufficiently large, irrespective of the contact angle. Our study therefore allows the quantification of the flow mechanisms when including the wetting effects. We note that we do not discuss here the effects of gravity, as the capillary forces dominate the gravitational forces, and so the buoyancy effect is secondary (see e.g. Hilfer and Øren (1996) for more discussion). We also note that an important characteristic of the liquid–solid interaction is the contact angle hysteresis, caused mainly due to surface roughness and heterogeneity. While to avoid further complications, we ignore the effects, it is shown in Helland and Skjæveland (2007) that if contact angle hysteresis is assumed, major displacement mechanisms still remain unaltered.

We also pay particular attention to sharp interfacial jumps during the displacement, known as Haines jumps Haines (1930); Cieplak and Robbins (1988); Lenormand (1990). This phenomena can lead to a highly discontinuous rate at which the front progresses, identified as bursting, when it would otherwise progress slowly and continuously. Here we numerically investigate a two-dimensional invading fluid displacing another fluid in an arrangement of circular disks. We show the influence of the Haines jump on the propagation of the front, which is accompanied by a sharp pressure oscillation. We further investigate the effects of the interfacial tension and the wettability on the jump transient dynamics. Finally, we study the frequency of Haines jumps as a function of $\theta$ and $\mathrm{Ca}$. This information can be useful to predict the residual saturation of the displaced phase, which is a measure of the displacement efficiency.

2 Computational Setup

In this work, the flow equations are solved numerically and the pore geometry is modeled using an Immersed Boundary (IB) method O’Brien and Bussmann (2018) in place of body-fitted meshes. The incompressible Navier-Stokes equations for a two-component fluid system are solved, where the interface is tracked using a Volume-of-Fluid (VOF) method. A contact line boundary condition is imposed on surfaces of the immersed boundaries. This approach has many advantages; in particular it greatly simplifies the inclusion of complex domain geometries, as the need to produce a body-fitted mesh is eliminated.

The Navier-Stokes equations for a two component fluid include the momentum equation,

[TABLE]

and the mass conservation equation,

[TABLE]

where $\rho$ is the fluid density, $\mu$ is the fluid viscosity, $\mathbf{u}$ is the fluid velocity, $p$ is the pressure, and $\mathbf{f}_{\sigma}$ is a surface tension body force.

For an incompressible flow consisting of two immiscible phases, the density is expressed as $\rho(\gamma)=\rho_{1}+\gamma\quantity(\rho_{2}-\rho_{1})$, while the viscosity is computed from harmonic averaging, $\frac{\rho(\gamma)}{\mu(\gamma)}=\frac{\rho_{1}}{\mu_{1}}+\gamma\quantity(\frac{\rho_{2}}{\mu_{2}}-\frac{\rho_{1}}{\mu_{1}})$, where $\gamma$ is the fluid volume fraction governed by the VOF advection equation,

[TABLE]

which is solved using a suitable numerical scheme that limits numerical diffusion in order to maintain a sharp interface between the two fluids. In this work, the Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) scheme of Ubbink and Issa Ubbink and Issa (1999) is used to solve the discretized form of Eq. (3). The computation of $\mathbf{f}_{\sigma}$ is carried out using the Continuum Surface Force (CSF) method of Brackbill et al.Brackbill et al. (1992), which models surface tension as

[TABLE]

where $\sigma$ is the constant surface tension coefficient, and $\kappa$ is the interface curvature, computed as

[TABLE]

where $\vec{\textbf{n}}$ is the interface unit normal vector, defined as

[TABLE]

where $\tilde{\gamma}$ is the smoothed volume fraction field. The contact angle is imposed through the computation of the interface curvature near the IB; see O’Brien and Bussmann (2020) for details. Briefly, in order to compute the interface curvature in cells neighboring an IB cell, the orientation of the interface normal at the contact line, determined from the prescribed contact angle, is used for the calculation in Eq. (5), instead of a value of $\vec{\textbf{n}}$ calculated from Eq. (6). To compute $\gradient\gamma$ at each point where the numerical stencil intersects the IB, ghost-cell values of $\gamma$ are evaluated using a bilinear interpolation and a least-squares stencil to reflect the contact angle boundary condition.

Finally, to close the system of equations, the divergence-free constraint is imposed on the velocity field, which has the form

[TABLE]

Eqs. (1) and (7) are solved together using a fractional-step projection method Chorin (1967), based on the balanced-force algorithm of Francois et al.Francois et al. (2006). The details of the discretization are omitted here, and can be found in O’Brien and Bussmann (2018). For the simulations reported in this work, Eq. (1) is non-dimensionalized as follows,

[TABLE]

where the superscript ∗ denotes a non-dimensional quantity, $\lambda_{\rho}$ is the density ratio, $\lambda_{\mu}$ is the viscosity ratio, and the Reynolds number $\real$ and Capillary number $\mathrm{Ca}$ are defined as follows,

[TABLE]

[TABLE]

where $\absolutevalue{\mathbf{U}}$ is the characteristic velocity, taken here to be the bulk flow velocity. A no-slip boundary condition is employed on the IB surfaces. It is noted that when using a no-slip boundary condition, a stress singularity is present at the contact line which can lead to results being dependent on the grid size. This phenomenon is discussed in detail by Afkhami et al.in Afkhami et al. (2009, 2018). Although the dynamic contact line models developed in Afkhami et al. (2009, 2018) are not used here, we believe that the results presented for the front instability and penetration modes are not strongly altered by this effect. Nevertheless, the effect of the stress singularity warrants further study.

In this work, density and viscosity ratios are set to unity. The characteristic velocity, $\absolutevalue{\mathbf{U}}$, corresponds to $Q/A$, where $Q$ is the flow rate through the domain and $A$ is the cross-sectional area of the domain. $\overline{D}$ is the mean cylinder diameter, which is sampled from a normal distribution. We note that a normal distribution assumption may not be a representative of more realistic geometries, but nonetheless we use it for a better control of the generated sizes that can be accommodated in our finite computational domain. More representative distributions will be considered for future study.

3 Results and Discussions

The results are organized as follows. In Section 3.1, the primary mechanism of pore-filling, the Haines jump, as well as its two modes, are investigated. In section 3.2, the different displacement regimes are characterized as functions of $\mathrm{Ca}$ and $\theta$. In Section 3.3, the effects of the contact angle and capillary number on the shape of the propagation front of the invading phase are examined.

3.1 Haines Jump Modes

A simplified setup, depicted in Fig. 1, is used in order to isolate the primary mechanisms of pore-filling. Three different $\mathrm{Ca}$ are compared, and the results are illustrated in Figs. 2 to 4. Singh et al.Singh et al. (2017) identified two different geometry-dependent Haines jump modes for cylinders and spheres. In 2D, these modes effectively correspond to the ‘touching’ and ‘coalescence’ modes described in Cieplak and Robbins (1988, 1990). We use the constructed geometry in Fig. 1 to qualitatively demonstrate both modes. We find that both modes of Haines jump will only occur at certain values of $\mathrm{Ca}$. As the interfaces corresponding to \theta=$$$, \theta=, and $\theta= pass between the two left-most cylinders, they begin to ‘bulge’ until they make contact with the center cylinder (Fig. 2c). This formation of a new contact line on the center cylinder initiates a rapid reconfiguration of the interface due to unbalanced surface tension forces. This represents the first mode of the Haines jump. As $\mathrm{Ca}$ is increased, and as the interfaces approach the throats formed by the right-most cylinders, they merge with a secondary interface passing around a cylinder (Figs. 3d and 4d). This causes a rapid release of the capillary pressure gradient acting against the flow, which causes the invading fluid to rapidly fill the pore and move through the throat. This corresponds to the second mode of the Haines jump. The effect of the geometry, injection rate and capillary number on the Haines jumps were studied in detail by Armstrong et al.Armstrong et al. (2015), showing that the occurrence of a Haines jump could be predicted from simple geometrical arguments. Here we also show that at low invading phase contact angles, and when the media is not densely packed, it is then possible that neither mode of Haines jump will occur. For the first mode, the contact lines on the adjacent sides of an obstacle will simply merge, and the protruding of the interface from the throat necessary for the second mode to occur will not be observed. In this scenario, complete displacement of the defending phase will be observed. When a Haines jump occurs, a pathway forms for the fluid to travel, in which no resistance from capillary pressure is present. As a result, pockets of defending fluid may remain in place, since the invading fluid will simply flow around them through sites in which Haines jumps have occurred.

As seen in Fig. 4, a viscous boundary layer forms initially, preventing a contact line to form on the middle cylinder before the contact lines traversing the surface of the left-most cylinders merge. At low $\mathrm{Ca}$, as seen in Fig. 2, the driving pressure may be too small to overcome the capillary forces at the throat, and the interface stalls, effectively blocking flow through the throat. At higher $\mathrm{Ca}$ (see Figs. 3 and 4), the liquid film trapped along the solid boundary collapses into a droplet upon the formation of contact lines for lower invading phase contact angles. Decreasing $\mathrm{Ca}$ results in a more stable contact line, with a $\theta$ that corresponds more closely to the prescribed equilibrium value.

We also note that the resolution needs to be sufficient to appropriately resolve the small features noticed above. Our developed numerical schemes provide a robust and high accuracy solution to capture these small structures within the pores. References (O’Brien and Bussmann, 2018) and (O’Brien and Bussmann, 2020) detail the numerical implementation in this work that allows to accurately capture such features. However, we have not investigated the required resolution to resolve these small length scales, as we do not believe that the phenomena of interest here, which is the way that the large length scale front invades the pore under different contact angles, is influenced significantly by this kind of small scale structure.

3.2 Displacement Regimes

The primary motivation behind this work is to investigate the effect that varying invading phase contact angle has on the fraction of the defending fluid that has been displaced from the domain in a porous media model. In order to investigate this, a random array of randomly-sized cylinders was generated, with different prescribed contact angles ranging from \theta=$$$ to \theta=. The secondary phase is injected into the cylinder array at a constant rate, at $\mathrm{Re}=5$, where the characteristic length is defined as the mean cylinder diameter. The displacement regimes include viscous fingering (high $\mathrm{Ca}$), capillary fingering (low $\mathrm{Ca}$, $\theta>) and stable displacement (low $\mathrm{Ca}$, \theta<$$$). The existence of these regimes has been discussed by previous authors Zhang et al. ([2011](#bib.bib34)); Ferrari and Lunati ([2013](#bib.bib12)); Mei et al. ([2012](#bib.bib21)); Zacharoudiou et al. ([2018](#bib.bib33)); Lenormand ([1990](#bib.bib19)); Jung et al. ([2016](#bib.bib18)) for different independent parameters. Injection patterns as a function of \mathrm{Ca}$and$\thetaare shown in Fig. [5](#S3.F5). Viscous fingering can be identified by the formation of long, thin fingers of the invading phase which often do not touch the obstacles as the relatively weaker surface tension results in the fingers simply following the velocity streamlines. Capillary fingers form contact lines with the obstacles, and are much rounder in shape due to the relatively strong surface tension forces present with high curvatures. In stable displacement, the front is much flatter and the interface much less distorted. Based on the observed data, transitions between the three known displacement regions can be identified, in a similar manner to how Zhang et al.Zhang et al. ([2011](#bib.bib34)) identified the existence of different displacement regions for a range of\mathrm{Ca}and viscosity ratios. Fig. [5](#S3.F5) depicts the transition between different regions in terms of\mathrm{Ca}$and$\theta$.

As seen in Fig. 5, as the viscosity ratio is held constant at one, all cases of $\mathrm{Ca}>0.1$ represent examples of viscous fingering, with little difference observed between media at different contact angles. As $\mathrm{Ca}$ drops and $\theta$ begins to increase, capillary fingering behavior is observed. Stable fronts are generally characterized by low $\mathrm{Ca}$ and low invading phase contact angles. These regimes were not observed to change much in our numerical experiments when slightly varying $\real$; nevertheless, for the applications we consider in this work, we keep $\real$ small. However, investigating the effects of varying $\real$ on the results is of interest and since it is beyond the scope of this study, we leave it for future work.

3.3 Displacement Patterns

Here we demonstrate the progression of the front through the porous media described in the previous section, for $\mathrm{Ca}=10^{-3}$. Figs. 6, 7 and 8 depict snapshots of the fluid injection into the porous media for contact angles of \theta=$$$, \theta= and $\theta= respectively.

The transition from stable to unstable displacement is clearly visible. At small contact angles, the front advances through nearly all available channels in part due to capillary imbibition action. As the porous media become increasingly non-wetting, i.e.  high invading phase contact angles, the adverse capillary pressure gradient causes fingers to form, as well as pockets of defending fluid which are never displaced. The presence of Haines jumps can be detected from the pressure differential across the porous media, since occurrences are accompanied by sudden rapid decreases in local capillary pressure. The pressure signals for various $\mathrm{Ca}$ and $\theta$ are shown in Fig. 9.

At high $\mathrm{Ca}$, the dynamics of the flow are dominated primarily by viscous forces. As a result, the average pressure gradient remains relatively constant, and varies only slightly as the flow establishes a steady state. In the viscous regime, the dynamics of the interface are nearly independent of contact angle. As $\mathrm{Ca}$ is lowered, two main things are observed to occur. First, the pressure differential needed to drive the flow increases with $\theta$. For a non-wetting porous media, a high mean curvature interface forms as it travels through a throat, resulting in a large capillary force acting in the direction opposite to the flow. For a more wetting porous media, i.e. low invading phase contact angles, the direction of the curvature is reversed, which acts to expedite the progression of the interface through the throat. The second observation is that due to the presence of Haines jumps and capillary imbibition, the intensity of pressure fluctuations increases; as we discuss next, when capillary fingering is present, the fraction of the defending fluid that has been displaced from the domain decreases in this case.

The degree to which capillary forces play a role in the flow dynamics can also be depicted by examining the front propagation. In Fig. 10, the front displacement, defined as the farthest longitudinal location of the interface at each time, for different $\mathrm{Ca}$ at various $\theta$ is shown. For any contact angle, the front propagation is relatively unaffected when Ca is large. Decreasing the Ca results in a greater variance over time in the front propagation, due to imbibition at low Ca, and capillary fingering at higher Ca.

The dominance of capillary forces at low $\mathrm{Ca}$ manifests itself primarily through the unsteady propagation of the front. This is due to high-frequency front oscillations resulting from both capillary rise and Haines jumps, and lower frequency oscillations that are the result of the formation of fingers. Viscous regimes are characterized by a relatively smooth propagation of the front when compared to capillary-driven regimes.

Examining the displacement fraction, $S$, defined as the ratio between the invading fluid volume and the total fluid field volume, can yield useful information regarding the different displacement regimes. Figs. 11 shows the displacement fraction $S$ as a function time for various regimes regime; as shown, a gradual increase in $S$ occurs up to $S_{\text{max}}$, where $S_{\text{max}}$ denotes the displacement fraction after steady state is obtained. In the viscous regime, the shear forces cause the interface to form long fingers which primarily follow the local velocity field. As a result, the interface does not immediately make contact with the surfaces of obstacles and form contact lines, which act to fill the pores normal to the flow direction through capillary action. As a result, the time scale on which the throats and pores normal to the flow direction are filled is larger when a contact line does not form, as the interface will only fill pores at a rate governed by the flux through pores themselves.

As $\mathrm{Ca}$ decreases, increasing the dominance of surface tension, and therefore capillary pressure, the interface is forced to expand in the transverse direction and fill more pores regardless of the pore orientation with respect to the flow. This effect can be seen in Figs. 11b and 11c, where at lower $\mathrm{Ca}$, a more pronounced initial linear displacement phase is observed, and in contrast to Fig. 11a, the linear displacement phase can be seen to end abruptly. This is due to the sudden formation of pockets of the defending fluid resulting from Haines jumps. Only in the case of \theta=$$$ is full displacement observed, since this regime represents a fully stable front which displaces all of the defending fluid. We also note that here we present an integrated quantity S$, which is meant to illustrate the amount of fluid displaced from the porous media by the invading phase, as opposed to identifying an individual Haines jump. The displacement is most efficient for a stable front propagation compared to the extremes where viscous or capillary fingering is present.

4 Conclusions

This work presents an extensive computational study of the dynamics of a two-phase flow in a porous media model. Direct numerical simulations provide an enhanced understanding of the subtle and complex competition between capillarity and wetting effects in a porous media consisting of two immiscible fluids displacing in a two-dimensional micro-geometry filled with randomly sized and randomly distributed cylinders. We show the transition from a stable front to a patterned displacement, with a focus on the wetting features of the porous media. The study provides detailed insight into various scenarios as we change $\mathrm{Ca}$ and $\theta$, while keeping the viscosity ratio fixed. We leave the study of the viscosity ratio effect to future work. We also study the Haines jump, namely rapid pore-scale displacement. Our simulations reveal the characteristic time and length scales of the Haines jump. Our numerical analysis of the Haines jump in a simple pore configuration, where cylinders of equal size are placed at the vertices of equilateral triangles, provide further insight into the effect of the contact angle at which the Haines jump is predicted. Finally, we do not address issues related to the moving contact line problem Afkhami et al. (2018), which is out of the scope of this study. An extension of this work should consider the effects of slip, and validate the results with experimental measurements. One of the main advantages of the direct numerical VOF/IBM method presented in this work is its ability to model complex geometries. Future work will consider the representation of three-dimensional domains to demonstrate the computational capabilities of our numerical methods with respect to other existing methods.