Accepted Manuscript

Effects of intermediate wettability on entry capillary pressure in angular pores

Harris Sajjad Rabbani, Vahid Joekar-Niasar, Nima Shokri

PII: DOI:

Reference:

S0021-9797(16)30195-3 http://dx.doi.org/10.1016/jotis.2016.03.053 YJCIS 21177

To appear in:

Journal of Colloid and Interface Science

Received Date: 17 February 2016

Accepted Date: 24 March 2016

Please cite this article as: H.S. Rabbani, V. Joekar-Niasar, N. Shokri, Effects of intermediate wettability on entry capillary pressure in angular pores, Journal of Colloid and Interface Science (2016), doi: http://dx.doi.org/10.1016/ j.jcis.2016.03.053

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Effects of intermediate wettability on entry capillary pressure in angular pores

Harris Sajjad Rabbani, Vahid Joekar-Niasar, Nima Shokri*

School of Chemical Engineering and Analytical Science, The University of Manchester,

Manchester, United Kingdom

Corresponding author Dr. Nima Shokri School of Chemical Engineering and Analytical Science Room C26, The

ineering and Ai

ie Mill

/ ersity of ] cville St Tel: 044

The University of Manchester

Sackville Street, Manchester, M13 9PL, UK

l: 0441613063980

Email: nima.shokri@manchester.ac.uk

Abstract

Entry capillary pressure is one of the most important factors controlling drainage and remobilization of the capillary-trapped phases as it is the limiting factor against the two-phase displacement. It is known that the entry capillary pressure is rate dependent such that the inertia forces would enhance entry of the non-wetting phase into the pores. More importantly the entry capillary pressure is wettability dependent. However, while the movement of a meniscus into a strongly water-wet pore is well-defined, the invasion of a meniscus into a weak or intermediate water-wet pore especially in the case of angular pores is ambiguous. In this study using OpenFOAM software, high-resolution direct two-phase flow simulations of movement of a meniscus in a single capillary channel are performed. Interface dynamics in angular pores under

ied. Inter]

nstant flow

drainage conditions have been simulated under constant flow rate boundary condition at different wettability conditions. Our results shows that the relation between the half corner angle of pores and contact angle controls the temporal evolution of capillary pressure during the invasion of a

pore. By deviating from pure water-wet conditions, a dip in the temporal evolution of capillary pressure can be observed which will be pronounced in irregular angular cross sections. That enhances the pore invasion with a smaller differential pressure. The interplay between the

contact angle and pore geometry can have significant implications for enhanced remobilization

of ganglia in intermediate contact angles in real porous media morphologies, where pores are

very heterogeneous with small shape factors.

1. Introduction

Understanding of the pore scale physics of immiscible displacement is imperative for maximizing efficiency of many industrial applications, most notably in petroleum industry for enhanced oil recovery (1; 2; 3) as well as soil remediation practices (4). The capillary forces

acting on the interfaces separating two immiscible fluids play the primary role in subsurface capillary trapping of hydrocarbons and low recovery factors. To mobilize the capillary trapped

ganglia, the mobilization force should overcome entry capillary pressure, which is controlled by wettability conditions and pore size (1). From computational point-of-view, entry capillary pressure is an essential entity required for simulation of two phase flow in pore network models

Delineating the impact of different variables on entry capillary pressure in angular pore geometries is not as straightforward as in cylindrical pore shapes (6). Due to the presence of angular corners, interface dynamics and its mathematical formulation become complicated (7). To calculate the entry capillary pressure, Mayer and Stow (8) devised an ingenious technique that relied on force balance in non-circular pore geometries. The technique was later further developed by Princen (9) and is referred to as Mayer-Stow-Princen (MS-P) approach. The principle of MS-P approach has been applied in many analytical (10; 5; 11) and semi-analytical approaches (12; 13) to compute entry capillary pressure in different polygonal cross sections. However, despite its significant value and its applicability in pore-network modelling, MS-P can

only be employed under equilibrium conditions. ?

Direct two phase flow simulation provides a unique opportunity to investigate interface dynamics under different wettability conditions at different pore geometries. Unlike traditional pore-network modelling approach where phase distribution is controlled by "local filling rules"

(14), the direct two-phase flow simulation relies on energy balance to simulate displacement process. Ferrari and Lunati (15) performed direct two phase flow simulation to study the impact

of capillary number (Ca) and viscosity ratio (^22.) on entry capillary pressure in strongly watered is

wet porous system, where |idis is the viscosity of displacing fluid and |iinv is the viscosity of the invading phase. Capillary number is defined as the ratio of viscous to capillary forces: Ca =

mmvQinv, where ^inv is viscosity of invading fluid, qinv is velocity of invading fluid and a is the

lillary forc fluid and

interfacial tension. The numerical simulation results demonstrated insensitivity of entry capillary pressure to viscosity ratio and capillary number. This is due to the two-dimensionality of the domain, which is not representative for a real porous medium. Further, Raeini et al. (16) formulated an algorithm to simulate two-phase flow in star-shaped pore channels (17), showing increasing entry capillary pressure as the pore body to pore throat contraction ratio increases. Similar to (15), the pore scale study conducted by Raeini et al. (17) was restricted to perfectly water-wet conditions.

None of the former studies have investigated the complex interplay of contact angle and pore geometries (corner angle) at intermediate contact angles under dynamic conditions. Since the conventional methodology for calculating the entry capillary pressure - based on balance of forces - can only be applied to a very limited range of water-wet pores, transient analysis of capillary pressure evolution of fluid-fluid interface displacement is required. We analyzed the dynamic evolution of a 3D meniscus at the entrance of a partially-wet pore. Within this context, we delineated the influence of pore angularity on the behaviour of fluid-fluid interfaces by performing finite-volume based numerical simulations using OpenFOAM (which employs "volume of fluid" method to capture the interface). All simulations were performed at Ca of 10- ; thus the capillary forces were dominating viscous forces. In the following, we first present the

governing equations and the numerical methods. Then the simulation results illustrating the dynamics of fluid-fluid interface displacement at different pore geometries and wettability at pore level are presented followed by the interpretation of the results and conclusions.

2. Numerical Model

2.1 Governing equations

ulate d

We use OpenFOAM (Open Field Operation and Manipulation) to simulate dynamics of two-phase flow in porous media. The code has been developed in C++ programming language and has been successfully implemented in several porous-media related research studies (17; 15).

The equations governing incompressible and immiscible displacement in a pore include mass and momentum balances, as follows:

Continuity equation reads V.u = 0 (1)

Momentum balance equation reads

+ V ■ (puu) = -Vp + V ■ (//(Vu + VitT) + fsa (2)

where i the vis

1* 3 12*

is the velocity vector (m s-), p is the density (kg m-), p is the pressure (kg m- s-), ¡/ is

viscosity (kg m-1 s-1) and fsa is the body force due to the capillary forces acting at the

interface (kg m- s- ). In the above system of equations the influence of gravity on flow domain is eliminated and fluid properties (density and viscosity) will vary depending upon the phase present in the computational cell (15).

For interface tracking, the Volume of Fluid (VOF) method introduced by (19) has been used. VOF is relatively straightforward to implement and can accurately model complex interfacial geometries (20). It involves computation of governing equations using volume indicator function (y), which represents fraction of phases in each grid block as;

(0,1) Interface; ye{ [1] Fluid 1; [0] Fluid 2.

.inction

At the interface the volume weighted fluid properties that includes density and viscosity are coupled to VOF algorithm as;

P = YPi + (1 - Y)P2 ¡1 = YHi + ( 1 - Y)ii2 The volume indicator function follo

______ws the a

advection transport as shown in Eq. (5):

— + V ■ (yü) + V ■ (7(1 - Y)Ur) = 0

where is the density of fluid 1 (kg m-3), p2 is the density of fluid 2 (kg m-3), ^is the viscosity of fluid 1 (kg m-1s-1), ¡i2 is the viscosity of fluid 2 (kg m-1s-1) and ur is the relative velocity between two fluids (m s-1). One of the drawbacks of VOF approach is the numerical diffusion of interface, which causes the interface to spread over many cells. This problem is referred to as smearing (21). Eq. (5) is the modified version of conventional volume of fluid transport equation as it includes additional convection term, known as artificial compression (21) (the last term in Eq. (5)). The motive of introducing this term is to minimize the smearing problem at the fluid-

fluid interface which can otherwise be quite cumbersome to resolve with conventional volume of fluid method (22; 21).

In addition to interface capturing algorithm, the surface force fsa acting as a source term in the momentum equation was described by Continuum Surface Force (CSF) method (23) that can be

represented as;

fsa = (6)

Where ox2 is the interfacial tension between fluid 1 and fluid 2 (kg s- ) and k is the curvature of the interface (m ) computed as;

:e volume nd y (17

Eqs. (1-5) are discretized using finite volume approach and solved numerically at each time step

to obtain primary unknowns p , u and y (17). Furthermore, the coupling between pressure and velocity equation was based on PISO (Pressure Implicit with Splitting of Operator) algorithm

developed by (24).

2.2 Modelling specifications

The numerical domain represents a pore connected to a large reservoir. Different cross sections including square, equilateral triangle and irregular triangle have been considered. The 3D views of pores are shown in Fig. 1(a-c) and the side view of the pore and its connecting reservoir has been shown in Fig. 1(d). Each pore can be described by its shape factor (G) expressed as the ratio of cross-sectional area to the square of perimeter. Furthermore, in all simulations pores are kept much longer than the reservoir. Fluid properties are shown in Table 1.

Fig. 1. (a-c) 3D views of the pores with different cross sections. (d) A side view of the simulation domain model. Invading phase (fluid 2) is displacing the receding phase (fluid 1) at a capillary number of 10-7.

uid 2) is

Table 1. Fluid properties used in the simulations

Parameter

Interfacial tension

Receding phase viscosity

Invading phase viscosity

Receding phase density

Invading phase density

Symbol

0.0001

Kg m-1 s-

Kg m-1 s-Kg m1 s

2.3 Initial and boundary conditions

The boundary conditions implemented in simulations can be seen in Fig. 1(d). Initially the interface is located 0.1 mm away from the junction inside the reservoir. As simulation starts the invading phase (fluid 2) is injected in the z direction at a constant flowrate and displaces the receding phase (fluid 1), while the outlet was kept at atmospheric pressure. The equilibrium contact angle, 0, is defined through receding phase imposed to the solid boundaries. In addition, to eliminate the hysteresis effect the advancing and receding contact angles were kept equal and no-slip boundary condition was imposed on the walls (25). Contact angles of 0 = 10o, 45o and 60o are simulated. At 0 = 45o the corner interface would be flat in square pore geometry and this would happen in equilateral triangle pore with 0 = 60o.

2.4 Spatial and temporal discretization

In order to avoid dependency of the interface thickness on the grid size, numerical simulations at different grid resolutions have been carried out. The results of grid independence tests are presented in Fig. A1, showing the variation of versus length of the simulation domain. The green curves show the grid size where the interface thickness does not change with grid resolution. Number of grid blocks per simulation are shown in Table A1. For a typical run the

numerical time step was made adjustable according to the courant number , where A t is the

time step and Ax is the grid block size. The maximum courant number was kept at 0.15 for square geometry and 0.4 in the case of regular and irregular triangle geometry.

ltel CPU

The simulations were performed using 12-node Intel CPU cluster (Xeon X7542 2.67GHz) with 4GB of memory. Only the drainage process was simulated, where invading fluid was acting as a non-wetting phase and receding fluid as a wetting phase unless stated otherwise.

2.5 Calculation of interface capillary pressure

The post-processing of numerical results was conducted with Paraview; an open-source

visualization application (26). Capillary pressure at each time step can be obtained by calculating the curvature of the interface;

Pc = 012*

where pc is the capillary pressure (kg m- s-).

Following (27), the iso-surface of 7=0.5 is considered as the fluid-fluid interface. Then the built-in Paraview plugin was used to evaluate the distribution of mean curvature over the entire interface. Finally, the median value of curvature distribution was used to calculate the capillary pressure using Eq. (8). The complete post-processing route is shown in Fig. 2.

Fig. 2. Post-processing of simulations to calculate the capillary pressure. (a) Phase distribution with red and blue representing receding phase and invading phase respectively. All other colors indicate the transition zone. (b) The interface extracted from the transition zone corresponding to 7=0.5 (c) Curvatures along the interface with the values represented by the color map.

3. Results and discussions

3.1. Entry capillary pressure

During the drainage while the non-wetting phase enters the pore, wetting phase will move out from the pore. The interface between the fluids which moves normal to the flow direction is referred to as "main terminal meniscus (MTM)" (28). Also some parts of the wetting phase remains in the corners. The fluid-fluid interface in the corners is referred to as "arc meniscus (AM)" as shown in Fig. 3 (29). In cases where 0 is equal or greater than critical contact angle (the contact angle at which AM becomes flat) the AM disappears from the corners, resulting in entrapment of receding phase at the junction in form of pendular rings (see Fig. 3) (6). The critical contact angle for square geometry is 45o and in the case of equilateral triangle it is 60o.

When the contact angle is equal to half corner angle, the interface would be flat, indicating almost zero capillary pressure.

Fig. 3. Configuration of receding phase (red) and invading phase (blue) in pore models with different wetting conditions. Each image is labeled with receding phase saturation at the junction inside the pore. G indicates the actor. The front view of capillary is shown with cross-section of reservoir (bottom) and pore (top). Invading

shape fac phase is d

phase is displacing receding phase at Ca of 10- .

Variations in entry capillary pressure profile as MTM passes the junction are shown in Fig. 4.

Time evolution

% 400 0.

.. .. .... ...

..........

(a) io° contact angle

■ Square (G = 0.062) » Equilateral Triangle (G = 0.048) • Irregular triangle (G = 0.041)

_i_i_i_

_l_|_|_L_

0.0 500 r

0.4 0.6

Receding phase saturation (-)

S 3(10

♦A .

(b) 45° contact angle

■ f «Ai ■

' Square (G = 0.062) Equilateral Triangle (G = 0.048) Irregular triangle (G = 0.041)

0.4 0.6

Receding phase saturation (-)

JS 200

• / /

'■^.'•vV-.,...-. .

jAm.. _:_*

• i . il A i

• V 4 * A - 4 V 4 A A »

» » ■ i- ' ** -1 V4 'i . V-1-

« i i 1 V! » »

» 0.4 * * A it% A 4 M » 0.8

(c) 60° contact angle

- Square (G = 0.062) * Equilateral Triangle (G = 0.048) Receding phase saturation (-) * Irregular triangle (G = 0.041 )

Fig. 4. Evolution of entry capillary pressure at meniscus entering pores with different cross sections and contact angles. X-axis represents the cross sectional saturation of the invading fluid at the pore entrance.

In a water-wet pore (contact angle of 10o) the entry capillary pressure in square cross section increases monotonically until reaching the maximum value as it moves into the pore. However in the case of other two geometries, there is a small dip before the entry capillary pressure increases. This localized reduction in entry capillary pressure may correspond to enhanced loc and temporal unstable fluid meniscus. Instability in the meniscus triggered due to rapid change in the size of confined region, makes it less resistant against the flow and elevates the local velocities as seen in Fig. 5. It is evident from Fig. 5 that in square geometry the rise in velocity occurs quite late and is not steep compared to both triangular geometries; as a result the impact of instability on the dynamics of interface is more significant in triangular pore geometries than in square pore.

Fig. 5. Variation in interface velocity as the interface enters the pore with 0 = 10o. All velocity values are scaled with respect to 1.0 x 10-4 m/s.

To further illustrate the phenomenon, Fig. 6 has been presented which shows vector plot of flow domain as invading fluid enter the pore. The motion of interface can be described in four key steps. First, there is a coflow of immiscible fluids (both invading and receding fluid from the corners of junction flow in same direction), that allows the meniscus to relax against capillary

forces (30) as shown in Fig. 6(a). Second, as meniscus makes contact with the pore wall and inflates, the intensity of coflow process decreases due to change in trajectory of corner flow shown in Fig. 6(b). During this phase the entry capillary pressure increases linearly with decrease in receding fluid saturation. Third, flow of fluid in the corners is directed downward from tip to neck of MTM (31), at this point the entry capillary pressure is maximum as shown in Fig. 6(c). The numerical results presented in this investigation indicates that for a pore of given inscribed radius, the maximum entry capillary pressure is inversely related to the shape factor and contact angle which is in agreement with previous studies (32; 6; 10). Fourth, imbibition initiates at the neck of MTM inaugurating Plateau-Rayleigh instability (Fig. 6(d)); this period is clearly visible in Fig. 4(c) after reaching maximum value, the entry capillary pressure in irregular triangle geometry decreases. Evolution of the interface in different pore cross sections as it advances from reservoir into the pore under different wetting conditions are presented in Fig. A2.

Fig. 6. Vector plot showing flow direction of fluids as advancing phase (blue) penetrates and displaces receding phase (red) in the equilateral triangle pore. The black arrow represents flow direction for fluid residing the corner.

The oscillations exhibited in the entry capillary pressure profile near the junction as wettability changes from perfect to partial wet (0 = 45o and 60o) become more pronounced with decrease of the shape factor in angular pores. Inspection of Fig. 3 indicates that as the interface progress into the partially-wet pore, thermodynamically driven rupturing of receding phase film occurs at the

junction, this along with coflow process declines the capillary pressure trend of all three pore models as shown in Fig. 4(b-c). At 0 = 60o the capillary behavior of interface as a consequence of film rupturing changes from drainage to imbibition in both triangle pore models, whereas in square pore model the change is not observed. This phenomenon can be rationalized by examining the variations in morphology of the interface shown in Fig. A2. In equilateral and irregular triangle pore cross-section, the capillary forces imbibe the corners of interface while the MTM is just at the junction. During this phase MTM is almost flat causing the overall curvature of the interface to be dictated by the AM (change in curvature of AM as seen in triangle pore models; see Fig. 3). However, in square pores the rupturing of film is initiated while the MTM is above the junction, causing the entry capillary pressure to remain positive.

Results reported in Fig. 3 and 4 clearly manifest that in pore network models similar to the pore body and pore throat, it is important to model junction as a separate entity. Ignoring it can lead to unrealistic flow scenarios. Practical implication of our numerical results in terms of remobilization of trapped phase is demonstrated schematically in Fig. 7.

Fig. 7. Schematic of water flooding to mobilize the trapped oil phase. The flow direction is from left to right. The driving force acting on imbibition interface marked blue tends to increase the pressure gradient between oil and water phase at the drainage interface. To remobilize the oil, it is vital for the driving force to overcome the capillary forces that are resisting mobilization of the trapped phase (figure adapted and modified from (33)).

During water flooding, in order to invade the pore throat, the local pressure gradient across the interface should overcome the entry capillary pressure. More importantly under same wettability conditions the entry capillary pressure across different pore geometries varies significantly. For

example according to Fig. 4(c), at 0 = 60 in both triangular pore geometries the entry capillary pressure at the moment of entree can drop sharply, while this behavior is less visible in square pore geometry. This drop of entry capillary pressure may indicate an easier remobilization of the trapped nonwetting phase for intermediate saturations when the shape factor is smaller. It is worthwhile to mention that in real porous media and rock formation the shape factor can vary mostly in the range of triangular pores (32; 34).

3.2. Comparison with analytical solutions

In this section the maximum entry capillary pressure estimated from simulation results will be compared with the one derived from MS-P method and equation proposed in the present study.

The entry capillary pressure for rectangular cross sections can be calculated using the following analytical solution (2):

-(a + b) co s(0) + J(a + b)2cos20 + 4 -6 - a/2 cos (j + e) cosO)

— 6 — a/2 cos g + 0 j cos0)

where pc* is the entry capillary pressure (pascal), a is the length (m) and & is the width (m) of the capillary cross-section. For equilateral triangular geometries, the following equation can be used as reported in (6);

= a12 (cos6 + Jjj=(sin29 + n(l

n/6+t n/2

For irregular geometries, the entry capillary pressure can be calculated using Eq. (11) derived by (35) as;

_ £712(l+2V7rG)cO.

Sdi 1+J1+c4qs2 9 1 1+2V7TG

COS20 4 G

where i s the rad ius of inscribed circle (m) and D = 1 —j-j + 3 s in6cos6 -

Using simple trigonometry relationships, we proposed the following equation to estimate the maximum entry capillary pressure in angular pore shapes.

_ mrj a12cose[sin (TT-2[3)

sin ((3)

e due e due

ally in to the

where is the half corner angle (degrees) and is the perimeter of the cross section (m). The comparison between entry capillary pressures predicted by theoretical equations against simulation results are presented in Fig. 8. Although the simulation results are generally in agreement with both theoretical predictions, there is slight discrepancy which could be following two reasons: Firstly, due to the 2D nature of the analytical solutions, the longitudinal curvature of the interface (in flow direction) has been assumed equal to zero (17). Secondly, our

simulation results indicate that even at Ca = 10- , the dynamic effects are present. Not only the contact angle will be different from the equilibrium one (see Fig. 8 (d)) but also the interface will not have equal mean curvature over the whole interface as shown in Fig. A2. Near three-phase contact line, the viscous forces commensurate with capillary forces that results in contact angle to be different from its actual value (36; 37). Moreover, as the corner of pore gets smaller the dissipation of viscous stress at three-phase contact line escalates; this enhances the deformation of interface caused by the dynamic e ffects (36).

Fig. 8. Comparison between theoretical and simulation results for (a) 10o, (b) 45o and (c) 60o contact angles. (d) Shows variation in interface shape as it advances towards the apex at different contact angle and corner angle.

4. Summary and conclusions

In this paper, direct two phase flow simulation was performed to investigate variations in the interface behaviour in angular pores under a wide range of wetting conditions. Our results offer new and fundamental insights regarding entry capillary pressure into pores; essential information for accurate description of flow in porous media. Under partially-wet conditions interface behaviour at the junction is highly unstable, inducing non-monotonicity in the entry capillary pressure trend, which was found to increase with pore angularity. The numerical results show that at 0 = 60o angular pores with smaller shape factors induce enhancement of the meniscus movement that leads to smaller entry capillary pressure. This is an important phenomenon as it suggests that in natural porous media composed of pores of different shape factors the impact of wettability on entry capillary pressure will be highly non-uniform. One may speculate that the remobilization of trapped oil phase will be easier in triangle or silt type pore geometries rather than circular or square pore throats. Taking into account the non-linear temporal evolution of the entry capillary pressure in the dynamic pore network modelling approach can allow manifestation of multi-phase flow scenarios pertinent to real oil reservoir rocks.

e in the < lti-phase flow sci

Acknowle

We would like to acknowledge the UK Engineering and Physical Sciences Research Council (EPSRC) to provide the PhD studentship for Harris Rabbani. We would also like to acknowledge the assistance given by IT Services and the use of the Computational Shared Facility at The University of Manchester.

Appendix A1. Grid independence

Showing Grid independence test results

Table A1. Number of grid blocks per simulation. The grid block size is scaled with respect to inscribed radius of pore shape.

Pore network Scaled grid block size Grid blocks

Square 0.05 105000

Equilateral triangle 0.10 91500 c

Irregular triangle 0.11 73000

Appendix A2. Iso-Surface map

Fig. A2. Demonstrating variation in morphology of interface as interface moves from the reservoir to pore under different contact angles.

Appendix A3. Derivation of proposed equation

Fig. A3. (a) Representing configuration of forces under hydrostatic equilibrium (b) the viscous forces at one of the corners of angular pore can be resolved in form trapezium.

Using sine rule one can show that

ct line

e, visc

sin (180 - 2(3)

sin ((3) ;

At three-phase contact line, viscous forces and capillary forces are comparable (37) so;

sin (180 - 2(3) 2(3

= sin ((3) ^s0(2nrO(—)

The pressure gradient at the interface (Pnw - pw) that is the total driving force (D) as soon as interface makes contact with solid wall can be written as

D - * (2-0

sin (180—2(3) sin ((3)

takes into account the angularity of pore shape. From force balance equation, one can

conclude that the total driving force equals to the summation of capillary force and viscous force. Thus capillary force can be written as

Capillary force = 51 n Jp) ffi2c OS0(2 un) - S'n p) ffi2c OS0(2 nr{)

sin (p)

sin (180-2(3) sin (p)

sin (p)

a12cos9(2TTri) (l-i^j

sin (180-213) ./„ 2B\

J, . sin(P) '^12^0(2^0(1-—)

Entry capillary pressure = -—--

It is important to note that for irregular pore shapes, smallest corner angle should be used.

smallei

allest cor

References

1. Dullien, F A L. Porous media: Fluid transport and pore structure. s.l. : Academic press, 1979.

2. Effects of flow history on oil entrapment in porous media: An experimental study. Khosravian, H, Joekar-Niasar, V and Shokri, N. 2015, American Institute of Chemica Engineers, Vol. 61, pp. 1385-1390.

3. Experimental study on nonmonotonicity of Capillary Desaturation Curves in a 2-network. Rodriguez de Castro, A; Shokri, N; Karadimitriou , N; Oostrom , Joekar-Niasar , V;. 10, 2015, Water Resources Research, Vol. 51, pp. 8517-8

4. Effect of soil texture on remediation of hydrocarbons-contaminated soil at El-Minia district Upper Egypt. Abdel-Moghny , T; Mohamed , RS; El-Sayed , E; Mohammed Aly , S; Snousy, MG;. 2012, ISRN Chemical Engineering 2012.

5. Simulating drainage and imbibition experiments in a high-porosity micromodel using an unstructured pore network model. Joekar-Niasar, V; Hassanizadeh, S M; Pyrak-Nolte, L J; Berentsen, C;. 2, 2009, Water resources research, Vol. 45.

6. Effect of contact angle on drainage and imbibition in regular polygonal tubes. Ma, S, Mason, G and Morrow, NR. 3, 1996, Colloids and Surfaces A: Physicochemical and Engineering Aspects, Vol. 117, pp. 273-291.

7. Menisci in arrays of cylinders: numerical simulation by finite elements. Orr, F M, Scriven, L E and Rivas, A P. 3, 1975, Journal of Colloid and Interface Science, Vol. 52, pp. 602-610.

8. Mercury porosimetry—breakthrough pressure for penetration between packed spheres. Mayer, RP and Stowe, RA. 8, 1965, Journal of colloid Science , Vol. 20, pp. 893-911.

9. Capillary phenomena in assemblies of parallel cylinders: II. Capillary rise in systems with more than two cylinders. Princen, H M. 3, 1969, Journal of Colloid and Interface Science, Vol. 30, pp. 359-371.

10. Threshold pressure in capillaries with polygonal cross Section. Lago, M and Araujo, M. 1, 2001, Journal of Colloid and Interface Science, Vol. 243, pp. 219-226.

11. Network model investigation of interfacial area, capillary pressure and saturation relationship in granular porous media. Joekar-Niasar, V; Prodanovic, M; Wildenschild, D; Hassanizadeh, SM;. 6, 2010, Water resources research, Vol. 46.

12. Analytical computation of arc menisci configuration under primary drainage in convex capillary cross sections. Held, M. 2, 2010, Computational Geosciences, Vol. 14, pp. 311-317.

13. A semi-analytical model for computation of capillary entry pressures and fluid configurations in uniformly-wet pore spaces from 2D rock images. Frette, O I and Helland, J I. 8, 2010, Advances in Water Resources, Vol. 33, pp. 846-866.

14. Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: A review. Joekar-Niasar, V and Hassanizadeh, SM. 18, 2012, Critical reviews in environmental science and technology, Vol. 42, pp. 1895-1976.

15. Direct simulations of interface dynamics: linking capillary pressure, interfacial area and surface energy. Ferrari, A and Lunati, I. University of Illinois at Urbana-Champaign : s.n., 2012. Proceedings of XIX International Conference on Water Resources.

16. Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method. Raeini, AQ, Blunt, MJ and Bijeljic, B. s.l. : Journal of Computational Physics, 2012, Vol. 231, pp. 5653-5668.

17. Numerical modelling of sub-pore scale events in two-phase flow through porous media. Raeini, AQ, Blunt, MJ and Bijeljic, B. 2, 2014, Transport in porous media, Vol. 101, pp. 191213.

18. Direct simulations of two-phase flow on micro-CT images ofporous media and upscaling of pore-scale forces. Raeini, AQ, Blunt, MJ and Bijeljic, B. s.l. : Advances in Water Resources, 2014, Vol. 74, pp. 116-126.

19. Volume of fluid (VOF) method for the dynamics of free boundaries. Hirt, CW and Nichols, BD. 1, 1981, Journal of computational physics, Vol. 39, pp. 201-225.

20. Volume of fluid methods for immiscible-fluid and free-surface flows. Gopala, VR and van Wachem, BGM. 1, 2008, Chemical Engineering Journal, Vol. 141, pp. 204-221.

21. Márquez Damián, S. Description and utilization of interFoam multiphase solver. [Online] http://infofich.unl.edu.ar/upload/3be0e16065026527477b4b948c4caa7523c8ea52.pdf.

22. Effects ofplasma flow velocity on melt-layer splashing and erosion during plasma instabilities. Miloshevsky, G and Hassanein, A. 2014, Nuclear Fusion, Vol. 54, p. 033008.

23. A Continuum Method for Modeling Surface Tension. Brackbill, J U, Kothe, DB and Zemach, C. 2, 1992, Journal of Computational Physics, Vol. 100, pp. 335-354.

24. Solution of the implicitly discretised fluid flow equations by operator-splitting. Issa, RI. 1, 1986, Journal of computational physics, Vol. 62, pp. 40-65.

25. Contact angle effects on microdroplet deformation using CFD. Rosengarten, G, Harvie, DJ and Cooper-White, J. 10, 2006, Applied Mathematical Modelling, Vol. 30, pp. 1033-1042.

26. Ayachit, U. The ParaView Guide: A Parallel Visualization Application. s.l. : Kitware, 2015. p. 276.

27. A level set method for determining critical curvatures for drainage and imbibition. Prodanovic, M and Steven, BL. 2, 2006, Journal of Colloid and Interface Science, Vol. 304, pp. 442-458.

28. Coexistence of menisci and the influence of neighboring pores on capillary displacement curvatures in sphere packings. Masson, G and Morrow, NR. 2, 1984, Journal of colloid and interface science, Vol. 100, pp. 519-535.

29. Meniscus configuration and curvatures in non-axisymmetric pores of open and closed unifrom cross section. Mason, G and Morrow, N R. s.l. : In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences , 1987. Vol. 414, pp. 111-133.

30. A Molecular Dynamics Simulation of Capillary Imbibition. Martic, G; Gentner, F; Seveno, D; Coulon, D; De Coninck , J; Blake, TD;. 21, 2002, Langmuir, Vol. 18, pp. 7971-7976.

31. Entry flow into a circular tube of slowly varying cross-section. Jayaraman, G, Padmanabhan, N and Mehrotra, R. 1986, Fluid Dynamics Research, Vol. 1.

32. Capillary behavior of a perfectly wetting liquid in irregular triangular tubes. Masson, G and Morrow, NR. 1, 1991, Journal of Colloid and Interface Science, Vol. 141, pp. 262-274.

33. Dynamics of oil ganglia during immiscible displacement in water-wet porous media. Payatakes, A C. 1, 1982, Annual Review of Fluid Mechanics, Vol. 14, pp. 365-393.

34. 3-D pore-scale modelling of sandstones and flow simulations in the pore networks. Bakke, S and 0ren, PE. 02, 1997, Spe Journal, Vol. 2, pp. 136-149.

35. Extending Predictive Capabilities to network models. Oren , P E, Bakke, Stig and Arntzen, O J. 4, 1998, SPE Journal , Vol. 3, pp. 324-336.

36. Moving contact lines: scales, regimes, and dynamical transitions. Snoeijer, JH and Andreotti, B. 2013, Annual review of fluid mechanics, Vol. 45, pp. 269-292.

37. Dynamics of immiscible-fluid displacement in a capillary tube. Zhou , MY and Sheng , P. 8, 1990, Physical review letters, Vol. 64, p. 882.

Graphical abstract