International Journal of Multiphase Flow xxx (2010) xxx–xxx Contents lists available at ScienceDirect International Journal of Multiphase Flow j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j m u l fl o w Effect of fluids properties on non-equilibrium capillarity effects: Dynamic pore-network modeling Vahid Joekar-Niasar ⇑, S. Majid Hassanizadeh Department of Earth Sciences, Utrecht University, PO Box 80021, 3508 TA Utrecht, The Netherlands a r t i c l e i n f o Article history: Received 22 January 2010 Received in revised form 14 September 2010 Accepted 16 September 2010 Available online xxxx Keywords: Porous media Non-equilibrium capillarity effect Two-phase flow Dynamic pore-network model Drainage Imbibition a b s t r a c t We have developed a Dynamic Pore-network model for Simulating Two-phase flow in porous media (DYPOSIT). The model is applicable to both drainage and imbibition processes. Employing improved numerical and geometrical features in the model facilitate a physically-based pore-scale simulator. This computational tool is employed to perform several numerical experiments (primary and main drainage, main imbibition) to investigate the current capillarity theory. Traditional two-phase flow formulations state that the pressure difference between the two phase is equal to the capillary pressure, which is assumed to be a function of saturation only. Many theoretical and experimental studies have shown that this assumption is invalid and the pressure difference between the two fluids is not only equal to the capillary pressure but is also related to the variation of saturation with time in the domain; this is referred to as the non-equilibrium capillarity effect. To date, non-equilibrium capillarity effect has been investigated mainly under drainage. In this study, we analyze the non-equilibrium capillarity theory under drainage and imbibition as a function of saturation, viscosity ratio, and effective viscosity. Other aspects of the dynamics of two-phase flow such as trapping and saturation profile are also studied. Ó 2010 Elsevier Ltd. All rights reserved. area and they proposed the following equation for capillary pressure: 1. Introduction 1.1. Theory of non-equilibrium capillarity effect Classical equations for multiphase flow in porous media are based a central equation – capillary pressure–saturation relationship. This equation, which is written based on the thermodynamically-equilibrium assumption, is commonly written as: Pn —Pw ¼ Pc ðSw Þ ð1Þ in which, Pc is the capillary pressure, Sw is the saturation of the wetting fluid, and Pn and Pw are the nonwetting and wetting fluid pressures, respectively. In fact, there are two major assumptions in this equation: (a) fluid distribution and consequently capillary pressure is a function of wetting fluid saturation only and (b) fluids pressure difference is equal to capillary pressure (at all times and under all conditions). Regarding the first assumption, it is known that there is hysteresis in Pc–Sw curves obtained for drainage and imbibition processes. Using a thermodynamic approach, Hassanizadeh and Gray (1990, 1993b) have suggested that the non-uniqueness in the Pc– Sw relationship is indeed due to the absence of specific interfacial ⇑ Corresponding author. E-mail address: joekar@geo.uu.nl (V. Joekar-Niasar). Pc ¼ P c ðSw ; anw Þ ð2Þ A number of computational and experimental works have shown that Pc–Sw–anw surfaces for drainage and imbibition more or less coincide under a wide range of drainage and imbibition history. This means that inclusion of anw leads to the removal, or significant reduction, of hysteresis in Pc–Sw relationship (e.g. Reeves and Celia, 1996; Held and Celia, 2001; Cheng et al., 2004, 2007; Joekar-Niasar et al., 2008, 2009, 2010b; Porter et al., 2009). Regarding the second assumption underlying Eq. (1), it is now an established fact that Pn–Pw is equal to capillary pressure but only under equilibrium conditions (see Hassanizadeh et al. (2002) for an extended review of experimental evidences). According to Entov (1980) capillary pressure–saturation relationship is not unique and, even though it is obtained under equilibrium conditions, it is a function of the history of fluids movements. In fact, it depends not only on the volume fraction of each phase, but also on their microscale distribution and change of saturation with time. Non-equilibrium effects in capillary pressure can be of significant importance in industrial porous media, such as paper pulp drying process, where gradients of fluids pressure and fluids velocities are large (Lewalle et al., 1994). For non-equilibrium conditions, the following equation for the difference in fluids pressure has 0301-9322/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmultiphaseflow.2010.09.007 Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 2 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx been suggested (Stauffer, 1978; Kalaydjian and Marle, 1987; Hassanizadeh and Gray, 1990): Pn —Pw ¼ Pc  s @Sw @t ð3Þ where s, a non-equilibrium capillarity coefficient, is a material property that may still be a function of saturation, and fluid–fluid properties (e.g. interfacial area, contact line). Potential dependencies of s are discussed in the next paragraph. In this work, we focus on investigating the validity of Eq. (3) by means of numerical experiments and determining the dependencies of the coefficient s. 1.2. Non-equilibrium capillarity coefficient Eq. (3) has been the subject of many studies in recent years, computationally using Darcy-scale models (see e.g. Manthey et al., 2005; Das et al., 2006), and pore-scale models (see e.g. Dahle et al., 2005; Gielen et al., 2005; Joekar-Niasar et al., 2010a) as well as experimentally see e.g. Hassanizadeh et al., 2004; Oung et al., 2005; O’Carroll et al., 2005; Bottero and Hassanizadeh, 2006; Bottero, 2009; Berentsen and Hassanizadeh, 2006; Camps-Roach et al., 2010; O’Carroll et al., 2010). Hassanizadeh et al. (2002) reviewed extensively the experimental works in which nonequilibrium effects have been observed. Up to now, most studies have studied the variation of s with saturation. However, based on pore-scale insights Hassanizadeh et al. (2002) conjectured that s may be related to phase trapping, capillary blockage (and consequently interfacial area and contact lines), and contact angle. Furthermore, in a recent study, O’Carroll et al. (2010) stated that s may be a function of fluid/fluid/solid contact line friction. Dahle et al. (2005) developed a bundle-of-tube model to investigate the variation of s with variance of radii distribution and with saturation under drainage process. They found that s increases with decrease of wetting fluid saturation and with increase of variance. In another study, Joekar-Niasar et al. (2010a) investigated variation of s with viscosity ratio (M) and saturation during drainage using a dynamic pore-network model. Viscosity ratio is defined as M = ln/lw, where ln and lw are the viscosities of nonwetting and wetting fluid, respectively. Joekar-Niasar et al. (2010a) found a similar trend for variation of s with saturation as Dahle et al. (2005) did. In addition, they found that larger values of viscosity ratio results in larger values of s. In another study, Manthey et al. (2005) studied variation of s with permeability, heterogeneity and entry capillary pressure during drainage. They employed a continuum-scale two-phase simulator in which local capillary pressure was defined based on Brooks-Corey relation. Based on volumetric phase averaging, they determined the magnitude and behavior of s. Although some aspects of non-equilibrium capillarity effect has been studied, there is still much room for further studies. There is still no study on the combined effect of saturation, effective viscosity, and viscosity ratio on the variation of s or on the presence of hysteresis on s–Sw relationship. 1.3. Dynamic pore-network modeling We employ pore-network modeling as an upscaling technique from pore scale to macroscale. This approach is computationally not too expensive, compared to other pore-scale simulators such as Lattice Boltzmann or Smoothed Particle Hydrodynamics (see Joekar-Niasar et al., 2010a, for a comparison between pore-network modeling and Lattice Boltzmann). Dynamic pore-network models can simulate transient behavior of flow with time for various capillary numbers and viscosity ratios. Capillary number (Ca) is traditionally defined as the ratio of qinv  viscous forces of the invading fluid to capillary forces linrvnw . The first dynamic pore-network model reported in the literature was developed by Koplik and Lasseter (1985) to simulate imbibition process in a two-dimensional pore network. Later, several dynamic pore-network models were developed for various applications, such as simulating pressure field evolution with time in a two-phase drainage process (see e.g. Aker et al., 1998a; Aker et al., 1998b; Van der Marck et al., 1997; Dahle and Celia, 1999), effect of capillary number on residual water saturation as well as fractional flow behavior during drainage (see e.g. Singh and Mohanty, 2003; Al-Gharbi and Blunt, 2005), effect of flow regime on residual nonwetting fluid in two-phase imbibition processes (see e.g. Nguyen et al., 2006; Mogensen and Stenby, 1998; Koplik and Lasseter, 1985; Hughes and Blunt, 2000; Thompson, 2002), evaporation (e.g. Prat, 2002), three-phase flow (see e.g. Pereira et al., 1996) ganglia movement see e.g. Constantinides and Payatakes, 1996; Dias and Payatakes, 1986a,b) and its contribution to relative permeability (e.g. Avraam and Payatakes, 1995a), and fluid–fluid interface velocity under drainage (Nordhaug et al., 2003). Strong nonlinearity at the pore-scale causes severe numerical stability problems in dynamic pore-network models, such that simulation for unfavorable viscosity ratios or for capillary-dominated flow have often proven troublesome. For instance, Al-Gharbi and Blunt (2005) could not obtain a full consistency between fluid occupancy at equilibrium resulted from quasi-static simulations and dynamic simulations for the same boundary conditions. In order to avoid these problems and also to be able to simulate imbibition as well as drainage under a wide range of viscosity ratios and different flow regimes, a new computational algorithm and a robust pore-network model was developed in Joekar-Niasar et al. (2010a). This new model has the following main features: (a) In contrast to all other dynamic pore-network models (except for Thompson, 2002), we assign two separate pressure fields to the two fluids within pore bodies or pore throats, and thus assign local capillary pressure to pore bodies. (b) Pore-scale flow mechanisms such as snap-off and cooperative pore filling (Lenormand and Zarcone, 1983, 1984) during both drainage and imbibition are included. (c) A semi-implicit approach is used for the saturation update in order to obtain numerical stability in simulations even for capillary-dominated and viscous fingering regimes. Thus, the resulting set of equations for fluid pressures contain both advection-type terms (corresponding to viscous forces) and diffusion-type terms (corresponding to capillary forces). (d) To mimic the geometry of pore space in granular media, pore bodies and pore throats are represented as truncated octahedrons and parallelepipeds, respectively. Local capillary pressure is determined from local interface curvature, which is related to local pore body saturation. Details of the computational algorithm employed in our dynamic pore-network model are given in Joekar-Niasar et al. (2010a). They developed a pore-network model with cubic pore bodies and parallelepiped pore throats to simulate drainage processes under different pressure gradients imposed at boundaries. They analyzed effect of viscosity ratio on non-equilibrium capillarity coefficient for different saturations. In addition, nonequilibrium effects on specific interfacial area were studied. Their results were only for primary drainage. 1.4. Objectives The main objectives of this paper are as follows: Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx  Developing a dynamic pore-network model (called DYPOSIT) for simulating two-phase flow under both drainage and imbibition conditions. The model employs improved numerical algorithms for computing pressure fields and updating the saturation. Using these algorithms, the transient behavior of capillary pressure under different flow regimes and viscosity ratios for drainage and imbibition processes are simulated; such studies under imbibition have not been reported in the literature. To mimic larger pore spaces among the solid grains in a granular medium, the shape of pore bodies has been selected as truncated octahedrons. This is a modification of the geometry compared with Joekar-Niasar et al. (2010a), in which pore bodies were represented as cubes.  Studying the evolution of average fluids pressure difference (Pn– Pw) with time and saturation for M = 0.1, 1, 10 for primary and main drainage as well as main imbibition processes. As a result, variation of s with viscosity ratio (and consequently effective viscosity) and saturation under drainage and imbibition has been analyzed.  Analysis of the validity of Eq. (3) during imbibition as well as drainage. Up to now, behavior of s has been studied only during primary drainage. We are interested in analyzing existence of a unique s under conditions where the local equilibrium at interfaces is still assumed (so, for example, dynamics of contact angle is neglected). 3 (a) (b) 2. Model description (c) 2.1. Model features 2.1.1. Structure and geometry Coordination number is a topological property of the network, which is the number of pore throats connected to a given pore body. Although there are many irregular quasi-static pore networks (e.g. Bakke and Øren, 1997; Joekar-Niasar et al., 2009, 2010b; Raoof and Hassanizadeh, 2010), the majority of dynamic pore networks have a regular structure. We have also developed a regular three-dimensional lattice network with fixed coordination number of six. Pore bodies and pore throats are presented by ‘‘truncated octahedron” and ‘‘parallelepiped”, respectively. This allows simultaneous existence of both fluids in a single pore element. The octahedron pore bodies can be unequally truncated since they are connected to pore throats of different sizes, as shown in Fig. 1a. Truncated parts of the octahedron have the shape of square pyramids with base width of 2rij (where rij is equal to the radius of inscribed circle of the pore throat ij) as shown in Fig. 1a. The radius of the inscribed sphere of a pore body, Ri, and the radius 0 of the inscribed circle of the diagonal cross sectionpof ffiffi pore body, Ri , are shown in Fig. 1b. It should be noted that Ri ¼ 36 R0i . A cross section through the vertices of a pore body and the pore throat connecting them is shown in Fig. 1c. The size distribution of pore bodies is specified by a truncated log-normal distribution, with no spatial correlation, expressed by: f ðRi ; rnd Þ ¼ qffiffiffiffiffiffiffiffiffiffiffi "  2 # R pffiffiffi ln i 2 exp  12 r Rm nd " ! Rmax Rmin 2r 2r !# ln R ln R m m  erf pffiffiffiffiffiffiffi pr2nd Ri erf pffiffiffiffiffiffiffi 2 2 nd nd ð4Þ in which Rmin is the lower range of truncation, Rmax is the upper range of truncation, Rm is the mean of inscribed sphere radii, and rnd is the standard deviation. Radius and length of pore throats connecting the pore bodies are determined based on the size of neighboring pore bodies. Spacings between the pore body layers of the network in x-, y- and z-directions are chosen to be variable. Let Fig. 1. (a) Schematic presentation of a pore body and its connected pore throats. Truncated parts of the pore body have the width of 2rij, which is the inscribed radius of pore throat ij. (b) Cross sections of the pore body along the vertices and through the edges. Radius of inscribed sphere is denoted by R and radius of inscribed circle 0 in the cross section along vertices is denoted by R . (c) Cross section of two pore bodies and connected pore throats. Geometrical configuration for determining the pffiffi pffiffi pore throat radius (rij) based on pore bodies radii, R0i ¼ 26 Ri and R0j ¼ 26 Rj . the spacings between layers k and k + 1 in the three directions be denoted by kx,k, ky,l, and kz,m. Then, designating each pore body by its lattice indices, namely k, l, and m, lattice spacings are given as follows: kx;k ¼ pffiffiffi 2 maxfR0 ðk; l; mÞ þ R0 ðk þ 1; l; mÞ : l ¼ 1; ny ; m ¼ 1; nz g; k ¼ 1; nx ð5aÞ ky;l ¼ pffiffiffi 2 maxfR0 ðk; l; mÞ þ R0 ðk; l þ 1; mÞ : k ¼ 1; nx ; m ¼ 1; nz g; ð5bÞ pffiffiffi 2 maxfR0 ðk; l; mÞ þ R0 ðk; l; m þ 1Þ : k ¼ 1; nx ; l ¼ 1; ny g; ð5cÞ l ¼ 1; ny kz;m ¼ m ¼ 1; nz where nx, ny, and nz denote the total number of pore bodies in x, y, and z directions, respectively. If there is no spacing between vertices Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 4 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx of square cross-sections shown in Fig. pffiffiffi 1c, the distance between centers of pore bodies will be equal to 2 times the summation of cross section inscribed radii. Based on the length of a pore throat and sizes of its neighboring pore bodies, radius of the inscribed circle in a pore throat cross section is determined (for more detailed explanation, see (Joekar-Niasar et al., 2008). Consider two pore bodies i and j, with a center-to-center distance dij (see Fig. 1c), and cross-sections inscribed radii R0i and R0j , respectively. We define the dimensionless radii Rei and Rej as follows: Rei ¼ R0i =dij ; Rej ¼ R0j =dij  n ð6Þ Then, we can calculate dimensionless inscribed radius of the pore throat ij, reij , as follows: reij ¼ .i .j .i .i1=n þ .1=n j ; n>0 ð7Þ e i sinðp=4Þ R ¼ e i cosðp=4ÞÞn ð1  R .j ¼ ð8Þ e j sinðp=4Þ R e j cosðp=4ÞÞn ð1  R ð9Þ where n is a positive number, which can control the ratio between radii of a pore body and pore throat, referred to as ‘‘aspect ratio”. A larger n results in smaller pore throats (larger aspect ratios). 2.1.2. Fluids and network parameters and specifications Table 1 shows fluids and network properties used in the simulations. Viscosity of the nonwetting fluid is kept constant and equal n to 0.001 and viscosity ratio is defined as the M ¼ llw . Statistical properties of pore bodies, pore throats, and aspect ratio distributions are shown in Table 2. Aspect ratio is defined as pore body inscribed sphere radius divided by pore throat radius. Corresponding to Table 2, Fig. 2a and b show the aspect ratio distribution as well as pore body-pore throat size distributions, respectively. Aspect ratio is controlled by the parameter n given in Eq. (9). In this work, we have set n = 1.0. 2.2. Governing equations ð10Þ A flux Q aij is assigned to a pore throat ij for each fluid separately. A volume balance for fluid a in pore body i is written as: X a Dsa Vi i ¼  Q ij ; Dt j2N i a ¼ w; n Specifications Ri (mm) rij (mm) Rasp min max mean st. deviation 0.0077 0.0200 0.0125 0.0028 0.0048 0.0162 0.0084 0.0017 1.55 4.00 2.25 0.38 Q aij ¼ K aij Dpaij ; a ¼ w; n ð12Þ a where K ij is a function of geometry and fluid occupancy of pore throats and Dpaij is the pressure difference in fluid a across the pore throats. Summation of Eq. (11) for the two fluids, which are assumed to be incompressible, yields: X j2Ni Q tot ij ¼  X n Q ij þ Q w ¼0 ij ð13Þ j2Ni Note that compared to one single pressure per body and one effective conductivity per pore throat (see e.g. (Al-Gharbi and Blunt, 2005; Mogensen and Stenby, 1998), assigning different pressures and conductivities to each fluid has major advantages. For example, it allows us to include mechanisms related to the local capillary pressure (such as snap-off, counter-current flow) in simulations. w n Eqs. (10)–(12) form a determinate set to be solved for sw i ; pi , and pi . 2.2.2. Pressure field solver To reduce the computational demand, the general equations are i , defined as the saturareformulated in terms of a total pressure p tion-weighted average of fluid pressures in a pore body: w n n  i ¼ sw p i pi þ si pi ð14Þ Combining Eqs. (10) and (14) and using sni þ sw i ¼ 1, we get the following equations for pressures of wetting and nonwetting fluids: n c  pw i ¼ pi  si pi n w i þ si pci pi ¼ p ð15Þ ð16Þ Substitution of Eqs. (15) and (16) in Eq. (12) and subsequently in Eq. i : (13) results in an equation for p 2.2.1. General equations for two-phase flow The local capillary pressure for pore body i is defined as:  w pci ¼ pni —pw i ¼ f si Table 2 Statistical properties of the radii of inscribed spheres of pore bodies (Ri), and of inscribed circle of pore throats rij, and aspect ratio distribution Rasp. ð11Þ where Ni is the set of all pore throats connected to pore body i, Vi is the volume of pore body i, and sai is the saturation of fluid a in pore body i. The volumetric flux of fluid a in pore throat ij is given by an equation similar to Washburn formula: X j2Ni  X h n   c n w w   Kw K ij sw pi i  K ij 1  si ij þ K ij ðpi  pj Þ ¼  j2Ni     i w þ Kw  K nij sw pcj j ij 1  sj ð17Þ In this equation, the right-hand side as well as the coefficients of the left-hand side depend on local saturation only. This linear system of i by diagonally-scaled biconjugate gradiequations was solved for p ent method using SLATEC mathematical library (Fong et al., 1993). 2.2.3. Saturation update i , pressures of fluids can be back-calculated After calculating p from Eqs. (15) and (16), based on saturation values. Then, Table 1 Fluid and network properties used in the simulations. Parameter Symbol Value Unit Contact angle Interfacial tension Wetting fluid viscosity Nonwetting fluid viscosity Total no. of pore bodies in flow direction Total no. of pore bodies in lateral directions Domain size Permeability h 0.0 0.0725 0.0001, 0.001, 0.01 0.001 45 35 1.9  1.9  2.37 1.43  1012 Degree kg s2 kg m1 s1 kg m1 s1 – – mm3 m2 rnw lw ln nz nx,ny – K Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx (a) 5 6 Frequency (%) 5 4 3 2 1 3 3. 1 3. 2 3. 3 3. 4 3. 5 3. 6 3. 7 3. 8 3. 9 4 2 2. 1 2. 2 2. 3 2. 4 2. 5 2. 6 2. 7 2. 8 2. 9 1. 5 1. 6 1. 7 1. 8 1. 9 0 Aspect Ratio (b) 20 Pore throat 18 Pore Body Frequency (%) 16 14 12 10 8 6 4 2 0. 00 9 0. 01 0. 01 1 0. 01 2 0. 01 3 0. 01 4 0. 01 5 0. 01 6 0. 01 7 0. 01 8 0. 01 9 0. 02 7 0. 00 8 6 00 0. 0. 00 5 4 00 0. 3 00 00 0. 0. 0. 00 2 0 Radius (mm) Fig. 2. Network geometry properties (a) Aspect ratio distribution, (b) pore body and pore throat distributions. commonly, Eq. (12) is used to calculate Q aij , and using Eq. (11) new saturation values can be calculated explicitly. This procedure, however, will result in numerical problems for a capillary-dominated flow regime, as mentioned in Koplik and Lasseter (1985). He found that the explicit saturation update was not numerically stable for very small capillary numbers and he could not successfully simulate a capillary-dominated flow. Therefore, we have developed a semi-implicit approach, analogous to a fractional flow formulation, to control the nonlinearities under such flow conditions. Details of the semi-implicit saturation update are given in Joekar-Niasar et al. (2010a). The resulting discretized equation reads: ! ! n w c c  kþ1 X K nij K w V i X K ij K ij @pij  w kþ1 ij @pij  si þ sw j tot Dt j2N K tot @sw @sw ij ij ij j2Ni K ij i ¼ n V i  w k X K ij tot si þ tot Q ij Dt j2Ni K ij where superscript k denotes time step level, and ð18Þ @pcij @sw ij is calculated pci —sw i relationship of the upstream pore body (given in Section from a 2.3.1). One should note that since Q tot ij and K ij are calculated from time step k, this scheme is not fully implicit. Here also, the diagonally-scaled biconjugate gradient method from SLATEC mathematical library (Fong et al., 1993) was used to solve Eq. (18). to a minimum saturation. On the other hand, if local imbibition occurs, the wetting fluid saturation in a pore body can go back to 1. So, we calculate Dti for all pore bodies, depending on the local process, from the following formulas:  8  w < qVni sw for local drainage; Q ni > 0 i  si;min i Dt i ¼   : Vni 1  sw for local imbibition; Q ni < 0 i q ð19Þ Dtglobal ¼ minfDti g ð20Þ i where the accumulation rate of the nonwetting fluid is defined as P Q ni ¼ j2Ni Q nij . Then, the global time step is chosen to be the minimum of all Dti s. It should be noted that we have imposed a truncation criterion of w w 106 for sw i  si;min and 1  si in order to reduce the number of times the equations are solved. Also, note that in Eq. (19), there is a correspondence between saturation change (numerator) and the accumulation rate of nonwetting fluid (denominator). That is, when local saturation is close to the limits, the accumulation rate of nonwetting fluid is also very small. This means that Dti always remains finite and nonzero. 2.3. Local rules 2.2.4. Time step The time step is determined based on the time of filling of pore bodies by the nonwetting fluid (or wetting fluid), denoted by Dti. The wetting fluid saturation of a pore body varies between 1 and sw i;min , as we assume that a pore body may be drained down only 2.3.1. Capillary pressure curves for pore bodies and pore throats Local capillary pressure within a pore is a function of the curvature of fluid–fluid interface through Young–Laplace equation, regardless of whether drainage or imbibition occurs. However, Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 6 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx during imbibition, topology of fluid–fluid interfaces is much more complicated. Interfaces may be categorized in two different types: ‘‘arc menisci (AM)”, formed along the edges, and ‘‘main terminal menisci (MTM)”, spanning the pore throat cross section during drainage (Mason and Morrow, 1987). Behavior of ‘‘MTM” during drainage and imbibition is different. During drainage, a MTM is formed at the entrance to a pore throat within a pore body filled with the nonwetting fluid. However, during imbibition a MTM can be formed such that it spans over a pore body and connected pore throats; this is referred to as ‘‘cooperative filling” (Lenormand and Zarcone, 1984). In any case, for a given fluid–fluid interface position within a pore body, we can determine corresponding capillary pressure and saturation. Therefore, for a given pore body, a unique relationship between capillary pressure and local saturation can be obtained. Local pci —sw curves for drainage and imbibition are i discussed and derived below. 2.3.1.1. Drainage. During drainage, the invasion of a pore throat by the nonwetting fluid, and consequently the movement of main terminal menisci, is controlled by the entry capillary pressure. In drainage, an interface is located within a pore body. If a pore body is filled with both fluids, the wetting fluid is residing along edges of the pore bodies (see Fig. A.1b and c in Appendix A). The saturation of the pore body (i.e. volume of the wetting fluid divided by the volume of the pore body) depends on the prevailing capillary pressure. For a given capillary pressure, the curvature of interfaces in the edges of the pore body can be calculated and, consequently, the corresponding saturation can be estimated. In Appendix A, details of derivation of the (local) pci —sw i relationship for an octahedron pore body are presented. The resulting pci —sw i relationship in terms of the radius Ri of the inscribed sphere of the pore body i and other geometrical parameters is: pci ¼ 2rnw ji ; ji 8 sw sdr 3:5 > > 1 1 dr i i > þ R1 sw < rij  Ri i P si 1sdr i i ¼  a > > sw dr > : R1 sdri ; a ¼ 2:98sdr13:85 smin < sw i i < si i i i ð21Þ in which, sdr i is the wetting fluid saturation corresponding to the inis the scribed sphere of the pore body, given by Eq. (A.2), and smin i minimum possible saturation in a simulation given by Eq. (22). The resulting curve is shown in Fig. A.2a. Obviously, it is impossible to completely displace the wetting fluid from the corners of a pore body. We assume that each pore body has a minimum saturation sw i;min , which depends on the im(Pcglobal posed global pressure difference defined in Section 3.2) as well as the blockage of the invading fluid. The capillary pressure   blockage of the invading fluid Pceblock is also a global variable, de- 2.3.1.2. Imbibition. During imbibition, locations of interfaces are not only controlled by the geometry of pore bodies and pore throats, but also by local fluid configuration in pore throats connected to a pore body (Lenormand and Zarcone, 1984). One of the mechanisms that controls configuration of fluid–fluid interfaces in a pore body during imbibition is ‘‘cooperative pore filling”. Lenormand and Zarcone (1984) observed in micro-model experiments that, depending on the number of pore throats filled with the nonwetting fluid, the topology of local main terminal interface can change. They defined a so-called I-mechanism for pore filling, where I represents number of pore throats filled with the nonwetting fluid. They showed that with the decrease of I, radius of interface may decrease, resulting in an increase of local capillary pressure as shown schematically in Fig. 3. Some researchers such as Øren et al. (1998), Hughes and Blunt (2000) suggested simple algebraic equations for pore filling mechanism, regardless of the geometry of pore bodies and pore throats. These relations do not apply to all pore spaces. For example, they do not apply to highporosity domains with small aspect ratio, such as micro-model experiments of Pyrak-Nolte (2007), as illustrated in simulations of the imbibition process by Joekar-Niasar et al. (2009). We define local capillary pressure–saturation relationship as a function of saturation as well as pore throats filling. Details of derivation of local pci —sw i relationship during imbibition are given in Appendix B. Finally, we have suggested the following relationships which results in a trend consistent with experimental observations of Lenormand and Zarcone (1984): pci ¼ 2rnw ji ; ji 8  a  3:5  a > simb sw simb simb > 1 i i i > þ < r1ij  R1i si dr Ri 1simb sdr i i i ¼  a > sw > 1 1 i > : R sdr ; a ¼ 2:98sdr 3:85 i i i imb sw i P si imb smin < sw i i < si ð23Þ where simb is the pore body saturation at which one pore throat is i still filled with the nonwetting fluid. It should be noted that this saturation is not predetermined in the simulations. Its value for each pore body will be determined during the simulation. 2.3.2. Entry capillary pressure for a pore throat We assume that a pore throat will be invaded by the nonwetting fluid when the capillary pressure in a neighboring pore body becomes larger than the entry capillary pressure of the pore throat. For a pore throat with square cross-section, the entry capillary pressure can be calculated as follows (due to (Mayer and Stowe, 1965; Princen, 1969a,b, 1970; Ma et al., 1996): pce;ij ¼ rnw r ij h þ cos2 h  p=4  sin h cos h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos h  p=4  h þ sin h cos h ! ð24Þ fined to be the minimum entry capillary pressure of all pore throats neighboring the vicinity of the nonwetting fluid but not invaded by it yet. Thus, using the pci —sw i relationship given by Eq. (21), the local minimum wetting fluid saturation in a pore body may be determined as follows: dr sw i;min ¼ si  Ri 2rnw n o2:98sdri 3:85 min Pcglobal ; Pceblock ð22Þ A capillary pressure should be also assigned to a pore throat once it is invaded and both fluids are present. We assume that the capillary pressure of a pore throat pcij is equal to the capillary pressure of the upstream pore body. Fig. 3. Pore filling mechanisms during imbibition referred to I-mechanism. Wetting fluid is shown as hashed areas. With decrease of number of pore throats filled with the nonwetting fluid, radius of the interface may decrease. Consequently capillary pressure assigned to the interface can increase (due to Lenormand and Zarcone, 1984). Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx where h is the contact angle. 2.3.3. Conductivities of pore throats Conductivities of pore throats are determined based on their size and fluid occupancy. One of the following two cases may occur. (a) A pore throat is occupied by the wetting fluid only. For this case, the following equation was obtained by Azzam and Dullien (1977): Kw ij ¼ p  eff 4 rij 8l w ij l ð25Þ K nij ¼ 0 where lw is the viscosity of the wetting fluid, lij is the length of pore throat, and reff ij rffiffiffiffi 4 r ij ¼ p ð26Þ (b) A pore throat is occupied by both fluids. Then, following Ransohoff and Radke (1988), we can write: 4  p  c 4 r ij blw ij l p  eff 4 K nij ¼ r 8ln lij ij Kw ij ¼ ð27Þ ð28Þ reff ij r nw pcij 0sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 c2 1 @ rij  ð4  pÞr ij ¼ þ r ij A 2 p ever, we have added five buffer layers at each boundary to reduce the boundary effect. Thus the network has a length of 45 pore bodies along the main flow direction. The buffer layers are not included in the averaging window. 3.2. Boundary conditions For our simulations, we assume that the network is connected to a nonwetting fluid reservoir on one side and a wetting fluid reservoir on the other side. Fluid pressures are specified at these boundaries. Side boundary conditions are assumed to be periodic. For drainage and imbibition simulations the following procedure is followed: (a) Drainage: Pressure at the nonwetting fluid reservoir is fixed to Pntop and pressure at the wetting fluid reservoir is fixed to zero. The difference between the two boundary pressures during drainage is referred to as ‘‘global pressure difference” P cglobal ¼ Pntop . After the nonwetting fluid breaks into the wetting fluid reservoir a boundary condition for its pressure is needed. One may choose to set it equal to the wetting fluid boundary pressure. But, this will cause an unrealistic relaxation of the nonwetting fluid pressure field throughout the network. Instead, as an assumption we have chosen to set the gradient of capillary pressure within an invaded pore throat (connected @pc  to the wetting fluid boundary) to be equal to zero @lijij ¼ 0 . (b) Imbibition: Similar to the drainage process, pressure at the nonwetting fluid reservoir is fixed to P ntop and pressure at where rcij ¼ 7 ð29Þ the wetting fluid reservoir is fixed to zero. However, P ntop should be so small that imbibition process can occur continuously. During imbibition similar to the drainage, it is ð30Þ assumed that In Eq. (27), b is a resistance factor that depends on geometry, surface roughness, crevice roundness and other specifications of the cross section. Detailed explanation about b can be found in Zhou et al. (1997). As mentioned earlier, the pore throat capillary pressure pcij is set equal to the capillary pressure of the upstream pore body. @pcij @lij ¼ 0 as long as a pore throat at the nonwett- ing fluid boundary is filled with both fluids. 3.3. Drainage simulations 3. Simulations and analysis For primary drainage, the network is initially fully saturated with the wetting fluid. Simulation of drainage starts with raising the pressure of the nonwetting fluid reservoir, and establishing a global pressure difference, P cglobal , across the network. When the pressure difference is larger than the entry pressure of the largest pore throat connected to the nonwetting fluid reservoir, drainage starts. In quasi-static simulations, the nonwetting fluid reservoir pressure is increased in incremental steps so that the network is invaded in steps. At the end of each step, when there is no flow (static conditions), the overall saturation of the network and average specific interfacial area are determined. Then, global pressure difference is increased again. In dynamic simulations, the imposed P cglobal is so large that the whole network (or a large part of it) could be flooded in one step. The simulations are continued till the change of average saturation in a selected averaging window is negligible. For main drainage simulations, the initial saturation occupancy is based on the last snapshot of the quasi-static main imbibition experiment. The simulation procedure is otherwise similar to the primary drainage simulation. 3.1. Network size 3.4. Imbibition simulations To analyze Darcy-scale equations using pore-network models, the size of the pore-network should be at least one REV. REV size was determined by performing quasi-static drainage simulations in networks with different sizes but with the same statistical parameters. Our simulations show that the REV size for these statistical parameters is a cube with length of 35 pore bodies. How- Consider a pore network filled by the nonwetting fluid at the end of a (main or primary) drainage experiment. The wetting fluid is still present along edges of pore bodies and pore throats. Starting from an equilibrium condition, all pore bodies have the same capillary pressure. The global capillary pressure is decreased by reducing the nonwetting fluid reservoir pressure or increasing 2.3.4. Snap-off During imbibition, if the local capillary pressure in a pore throat becomes smaller than a critical value (defined below), the corner interfaces become unstable and snap-off will occur. The criterion for snap-off in a square cross-sectional pore throat has been defined as follows (Vidales et al., 1998): pcij 6 rnw r ij ðcos h  sin hÞ ð31Þ We assume that as soon as snap-off occurs, the nonwetting fluid retreats instantaneously into the two neighboring pore bodies, and the pore throat is filled up with the wetting fluid. Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 8 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx Our simulations result in local-scale variables such as pressure, saturation, and fluxes. These have to be averaged over an averaging window to obtain macroscopic variables. Average saturation is simply defined as the ratio of the wetting fluid to the total pore volume of the network: Pnpb w Vw i¼1 si V i n ¼ Pnpb V þV i¼1 V i ð32Þ w Sn ¼ 1  S w in which npb is the total number of pore bodies within the averaging window. The total flux across any given surface is equal to the sum of fluxes of all pore throats intersecting that surface. The averaging of pressure is, however, less straightforward. Commonly, average pressure is obtained using an intrinsic phase average operator (see e.g. Whitaker, 1977). However, recently it has been shown that the intrinsic phase average pressure introduces numerical artifacts when both pressure and saturation are spatially variable (see Nordbotten et al., 2007, 2008; Korteland et al., 2010; Bottero, 2009). Instead, a centroid-corrected averaging operator has been suggested by Nordbotten et al. (2008) to alleviate problems associated with intrinsic phase averaging. However, it has been shown numerically that this approach also suffers from artifacts at large wetting saturations (Korteland et al., 2010). Thus, we use the intrinsic phase average since it is still the most commonly used operator: Pa ¼ 1 dV a Z Pnpb pai sai V i Pa dV ¼ Pi¼1 ; npb a dV a i¼1 si V i a ¼ n; w 4.1. Non-equilibrium capillarity effects As mentioned in the Introduction, in standard models of twophase flow in porous media, the difference in average fluid pressures, Pn–Pw, is assumed to be an algebraic function of saturation. In this section, we show that the curves of Pn–Pw versus Sw depend strongly on boundary conditions and dynamics of the system, as shown by Yang et al. (2009) for bundle-of-tube model. We performed a large number of dynamic primary drainage, main drainage, and main imbibition simulations for various viscosity ratios and under a range of boundary pressure values. For every simulation, we calculated the variations of average saturation and the ‘‘fluids pressure difference” with time. From those results, plots of Pn–Pw versus saturation were made. For primary and main drainage simulations, three different global pressure differences were imposed over the domain: 45, 60 and 90 kPa. The entry capillary pressure is 15 kPa. For main imbibition, four different global pressure differences were imposed over the domain, equal to +5, 0, 5, 10 kPa. Also, all simulations were performed for three different viscosity ratios, M = 0.1, 1, and 10 (ratio of viscosity of the nonwetting fluid to that of the wetting fluid). Computational time increased with decrease of viscosity ratio M=0.1 30 25 PD 20 MD MI 15 10 0 0.2 0.4 0.6 0.8 1 Saturation (Sw) (b) 35 M=1.0 30 25 PD 20 MD MI 15 10 5 ð33Þ 4. Results and discussion 35 5 0 0.2 0.4 0.6 0.8 1 Saturation (Sw) (c) Average Pressure Difference (kPa) Sw ¼ (a) Average Pressure Difference (kPa) 3.5. Averaging procedure and decrease of global pressure difference. A simulation took from 24 to 72 h on Intel(R) CPU 6600, 2.4 GHz with 2GB RAM. Plots of Pn–Pw as well as Pc versus average saturation, Sw, are shown in Fig. 4. It is clear that Pn–Pw curves strongly depend on boundary conditions and significantly deviate from the Pc curves. Their behavior is in agreement with Eq. (3): at any given saturation, the fluids pressure difference is larger than the capillary pressure during drainage and smaller during imbibition. The differences are larger for larger nonwetting fluid boundary pressure during drainage and for large wetting fluid boundary pressure Average Pressure Difference (kPa) the wetting fluid reservoir pressure. A decrease in the global capillary pressure causes the interfaces to relax gradually and main imbibition experiment will start. The imbibition simulation will stop when all pore throats at the outflow boundary (nonwetting fluid reservoir) are fully filled with the wetting fluid. For the quasi-static imbibition simulations, an approach similar to the drainage has been employed. 40 M=10 35 30 25 PD 20 MD 15 MI 10 5 0 0.2 0.4 0.6 0.8 1 w Saturation (S ) PD:+45 kPa PD:+60 kPa MD:+45 kPa MD:+60 kPa MI:-100.2 kPa MI:+5kPa MI:-5 0.6 kPa w Saturation Pc-S(S Equil.) 0.4 PD:+90 kPa MD:+90 kPa 0.8 1 MI:0 kPa Fig. 4. Average fluids pressure difference curves, Pn–Pw for three different global pressure differences during primary and main drainage (PD, MD), and four different global pressure differences during main imbibition (MI), for viscosity ratios (a) M = 0.1, (b) M = 1 and (c) M = 10 and for ln = 0.001 Pa s. These curves should be compared to the equilibrium Pc–Sw curves resulted from quasi-static simulations. Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 9 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx Fluid entrapment at pore scale, and consequently fluid distribution at macroscale is controlled by pore-scale invasion mechanisms, such as piston-like movement and snap-off. Snap-off refers to the mechanisms in which residual fluid can cause a local rupture in the bulk fluid. The importance of these mechanisms varies depending on the process (drainage or imbibition). For instance, snap-off is more important during imbibition than drainage, since the nonwetting fluid (bulk fluid) is the receding one and it is more vulnerable to rupture. As a result, at the end of imbibition process, a significant amount of the nonwetting fluid remains behind, known as the residual saturation. Due to the importance of imbibition in reservoir engineering, the dependence of residual saturation on capillary number, viscosity ratio, contact angle and pores 0.4 0.5 According to Eq. (3), the deviation of non-equilibrium fluids pressure difference from the macroscopic capillary pressure is 0.3 4.2. Effect of viscosity ratio on fluid distribution 4.3. Dependence of s on fluid viscosities and saturation 0.2 It must be noted that the effective viscosity (or effective resistance) of a system is related to fluid viscosities as well as spatial distribution of fluids. The definition given in Eq. (34) is only suitable for a stable invading front where a piston-like movement is dominant, as the viscosities are linearly weighted by saturation. This equation is analogous to Kirchhoff’s law for electrical resistors. If the resistors are in series, the effective resistance is the summation of all resistors (similar to piston-like movement). However, unstable flow can be analyzed as a circuit with mixed combination of resistors in parallel and series. Every single capillary blockage increases the resistance of the system against flow. Consequently, resistance due to capillary blockage is more pronounced under unstable flow compared to stable one. For unstable fronts, e.g. in a viscous fingering regime, some researchers, such as Koval (1963), have suggested other empirical relationships for effective viscosity (see e.g. Fayers et al., 1990), which are not linearly weighted by saturation. During both drainage and imbibition, the effective viscosity changes considerably if fluid viscosities are not equal. During primary drainage, effective viscosity is initially equal to lw, which was chosen to be different in simulations with different viscosity ratios (i.e. nonwetting fluid saturation was kept constant and lw is varied to get various M values). As the drainage proceeds, effective viscosity decreases for M = 0.1 and increases for M = 10. At the end of drainage, the effective viscosity is almost equal to ln. According to Fig. 4 the effect of viscosity ratio on the transient behavior of fluids pressure difference is significant. During drainage, initially (i.e. at higher saturations), Pn–Pw is larger for small M and smaller for larger M (for equal ln). This is due to the fact that the pressure drop over the whole domain for M = 10 is smaller than the other cases. The decline in the fluids pressure difference during drainage as residual saturation is approached, as an artefact of the intrinsic phase averaging. As the nonwetting fluid invades the network, the wetting fluid pressure behind the invasion front increases along with the nonwetting fluid pressure. As a result, after breakthrough, the average fluid pressure difference may even decrease. Consequently, the average pressure difference approaches the Pc–Sw curve. These results illustrate that under dynamic conditions, phase pressures difference is not equal to capillary pressure–saturation curve obtained under equilibrium conditions. In the other words, the classical two-phase flow simulators, which do not include dynamic effects in their capillary pressure–saturation curve, may not hold under dynamic conditions. 0.1 ð34Þ Residual Nonwetting Phase (Sor) leff ¼ ln Sn þ lw Sw aspect ratio has been studied significantly (see e.g. Dias and Payatakes, 1986a,b; Vizika et al., 1994; Mogensen and Stenby, 1998; Hughes and Blunt, 2000). The typical curve relating residual (trapped) nonwetting fluid saturation to the capillary number is shown in Fig. 5. This curve was obtained from imbibition n simulations at M ¼ llw ¼ 1:0 under constant flux boundary conditions. The residual nonwetting fluid saturation was determined as a function of capillary number at various flow rates. This figure shows that with the increase of capillary number, the residual nonwetting fluid saturation decreases. The largest decrease occurs for capillary numbers between 106 and 104. At higher capillary numbers, large flow rate cause a piston-like movement and suppresses the snap-off mechanism. As mentioned before, all previous studies were performed under a constant capillary number. In this section, we show that under Dirichlet boundary condition, a variable spatial distribution of the trapped phase can be resulted, since the flow rate and capillary number change continuously in time and space. Our imbibition results (not presented here) show that, for M = 0.1, with the increase of wetting fluid saturation, effective viscosity increases, and consequently flow rate decreases. But, for M = 1 and M = 10 with the increase of wetting fluid saturation, flow rate increases. This increase is very pronounced for M = 10, which changes about two order of magnitude. As shown in Fig. 5, change of flow rate can influence the trapped nonwetting fluid saturation. This effect can be observed when examining the average saturation profile along the network, shown in Fig. 6a and b. Here, the saturation is averaged over a cross section of the network located at position x and then plotted against x/l at different times, and for two different viscosity ratios. An interesting result here is the non-monotonic behavior of saturation for M = 10. Moreover, the saturation front for the case of M = 10 is steeper than for M = 1.0. With the significant increase of flow rate for M = 10, snap-off is suppressed. Suppression of snap-off helps increasing the flooding efficiency. To show the effect of viscosity ratio on saturation profile during drainage, saturation profiles for M = 1.0 and 0.1 for DP = 45 kPa have been shown in Fig. 7. As it can be seen, for M = 0.1, the invading front is unstable and the slope of the front is much smaller than for M = 1.0. The comparison between invasion front for these two cases shows how significant the effect of local heterogeneities on saturation profile is. Consequently, the flooding efficiency is smaller for M = 0.1 compared with M = 10. M=1 0 during imbibition. These curves are used later to calculate nonequilibrium capillarity coefficient, s. Resulting Pn–Pw curves also depend strongly on viscosity ratio and effective viscosity. In order to explain this, let us define the effective viscosity, leff, as follows: 1.E-08 1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 Capillary number Fig. 5. Effect of capillary number on residual nonwetting saturation for n M ¼ llw ¼ 1:0. Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 10 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx (a) 1 (a) M=10 M=0.1 0.8 0.6 Saturation (Sn) 0.8 Saturation (Sw) 1 t=0.018s 0.4 0.6 0.4 0.2 0.2 t=0.31 s 0 0 0 0 0.2 0.4 0.6 0.8 0.2 1 (b) 0.6 0.8 1 Length (X/L) Length (X/L) (b) 1 0.4 1 M=1 M=1.0 0.8 0.6 Saturation (Sn) Saturation (Sw) 0.8 t=0.048 s 0.4 0.6 0.4 0.2 t=0.11 s 0.2 0 0 0 0 0.2 0.4 0.6 0.8 1 Length (X/L) Fig. 6. Wetting fluid saturation profiles during imbibition for DP = 10 kPa for (a) M = 10 and (b) M = 1. Dark-colored curves show the wetting fluid saturation distribution at the time of breakthrough of the wetting fluid. Light-colored curves show the saturation profile at earlier times. In all cases ln = 0.001 Pa s. related to the time rate of change of saturation. Hassanizadeh and Gray (1993a) have suggested that the non-equilibrium capillarity coefficient, s, is a function of saturation. It is also believed to depend on other medium and/or fluid properties. To compute s at a given saturation, a plot of Pn–Pw versus @Sw/@t at that saturation must be made. This is done by starting from one of the graphs in Fig. 4 for a given viscosity ratio. At a given saturation, different values of Pn–Pw–Pc can be obtained pertaining to different global pressure differences. Corresponding values of @Sw/@t are known for that saturation. Note that such plots are made for drainage and imbibition separately, whereby the corresponding (drainage or imbibition) Pc–Sw curve is used to calculate Pn–Pw– Pc. By plotting Pn–Pw–Pc versus @Sw/@t at various values of Sw and fitting a linear relationship to the data points, s is found as the slope of the fitted line. Results for primary and main drainage as well as for main imbibition are plotted in Fig. 8a and b, respectively. Under drainage, the fitted line passed through the origin for different M values, as it should. However, under imbibition, there was an intercept. This means that for the curve to go through the origin a nonlinear relationship should be employed. However, with the increase of the pressure gradient, this nonlinearity decreases, as the capillary forces are suppressed by the viscous forces. Two different aspects of these data are important; (a) order of magnitude of the non-equilibrium capillarity coefficient s and (b) its variation with saturation for different viscosity ratios as well as for drainage and imbibition. 0.2 0.4 0.6 0.8 1 Length (X/L) Fig. 7. Effect of viscosity ratio on invasion process during drainage. Nonwetting fluid saturation profiles during drainage for DP = 45 kPa for (a) M = 0.1 and (b) M = 1. For M = 0.1, the invasion front is unstable and the role of porous medium geometry in viscous fingering is obvious. In all cases ln = 0.001 Pa s. Values of s obtained here (ranging from 100 to 1000 Pa s for M = 1.0) are in agreement with results obtained in other pore-network modeling studies as shown in Table 3. The dimensions of our network correspond to a porous medium sample size of 1.9  1.9  1.9 mm3 with an intrinsic permeability of 1.43  1012 m2. While the permeability corresponds to a fine sand, the sample size is very small. The general understanding is that the magnitude of s increases with the size of the observation or averaging window and is inversely correlated with permeability. Dahle et al. (2005), using a bundle-of-tube model, concluded that s value can be proportional to L2, where L is the length of averaging window. A similar result was found by Manthey et al. (2005) based on simulations at continuum scale. In laboratory experiments by Hassanizadeh et al. (2004), the value of s for a fine sand sample of 3 cm in height was found to be 5  105 Pa s. The pressure measurements were actually done by transducers with a diameter of around one centimeter. Using similar transducers in experiments with the same fine sand, Bottero and Hassanizadeh (2006) found a s value of around 105 Pa s. But, when Bottero (2009) upscaled the results to the column scale (18 cm), the average s value was found to be around 106 Pa s. Our results shown in Fig. 8a and b illustrate that for the same nonwetting fluids, smaller viscosity ratio is, larger non-equilibrium capillarity coefficient will be. Since effective viscosity is defined as the sum of linear saturation-weighted viscosities, for the same saturation smaller viscosity ratio results in larger effective viscosity. Thus, we conclude that larger effective viscosity results in larger non-equilibrium capillarity coefficient (s(M = 0.1) > s(M = 1) > s(M = 10), ln = constant). This trend is obvious for drainage Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 11 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx (a) 1.E+04 (b) MI: M=0.1 MI: M=1.0 MI: M=10 1.E+03 (Pa.s) 1.E+03 (Pa.s) 1.E+04 1.E+02 PD:M=0.1 PD:M=1.0 PD:M=10 MD:M=0.1 MD:M=1.0 MD:M=10 1.E+01 1.E+00 0 0.2 0.4 0.6 0.8 1.E+02 1.E+01 1.E+00 0 1 0.2 w (c) (d) 1.E+07 1.E+06 1.E+05 1.E+05 PD:M=0.1 PD:M=1.0 PD:M=10 MD:M=0.1 MD:M=1.0 MD:M=10 1.E+04 1.E+03 0.2 0.4 0.6 0.8 1.E+07 1.E+06 0 0.4 1 Saturation (Sw) Saturation (S ) 0.6 0.8 MI: M=0.1 MI: M=1.0 MI: M=10 1.E+04 1.E+03 1 0 0.2 Saturation (Sw) 0.4 0.6 0.8 1 Saturation (Sw) Fig. 8. Variation of non-normalized (a and b) and normalized (c and d) non-equilibrium capillarity coefficient as a function of saturation for different viscosity ratios M = 0.1,1,10 for (a and c) primary and main drainage, (b and d) main imbibition simulations. s is normalized by effective viscosity defined as leff = (1  Sw)ln + Swlw. The vertical axis is in logarithmic scale and in all cases ln = 0.001 Pa s. Table 3 Values of s in literature for computational and experimental works. Reference Exp. type Process Fluids M Pd (kPa) K (m2)  1012 s (Pa s) @s @Sw Domain dimensions (cm) Dahle et al. (2005) Gielen et al. (2005) Gielen (2007) Joekar-Niasar et al. (2010a) Joekar-Niasar et al. (2010a) Das et al. (2007) Manthey et al. (2004) BTM PNM Drain Drain – Oil–water 1 10 0.8 5 4.7 – 274 1.2  105 0 <0 0.1 (L) 0.3  0.3  1 PNM PNM CM Lab Drain Drain Drain Drain – – Oil–PCE–water PCE–water 1,10 0.1 >0.6 0.9 4 4 1.2 5.58 150 150 5 3 900,2750 350 105–107 2–7  104 0 0 0 <0 >0 0 0 0 0.5  0.5  0.5 0.5  0.5  0.5 O’Carroll et al. (2005) O’Carroll et al. (2005) Bottero and Hassanizadeh (2006) Bottero (2009) Camps-Roach et al. (2010) Camps-Roach et al. (2010) 7 Lab Lab Lab Drain Drain Drain PCE–water PCE–water PCE–water 0.8 0.8 0.9 2.1 2.4 6 15.8 12.6 – 5.64  10 1.99  107 105–107 Lab Lab Drain Drain Air–water Air–water 0.02 0.02 2.2 4.6 53 14.7 2  105–106 105–8  105 0 0 – 5.07(D), 9.62(L) 5.07(D), 9.62(L) 9.8 (D), 19(L) 10 (D), 20(L) 10 (D), 20(L) BTM: Bundle of tubes, PNM: Pore-network model, CM: Continuum model, Lab: Laboratory experiment D: Diameter, L: Length. (Fig. 8a) and imbibition (Fig. 8b). This agrees with the explanation given by Barenblatt et al. (2003) and Entov (1980), who stated that dynamic effects in capillary pressure are related to the finite time required for the fluids in the pore structure to rearrange themselves. Indeed, for larger effective viscosity values, more time is required for fluids to reach the equilibrium condition, which corresponds to a larger s value. As it is known, dimension of s is the same as viscosity. Since resistance of the system is related to viscosities of both fluids and spatial distribution of fluids, effective viscosity can be used to normalize the non-equilibrium capillarity coefficient. The normalized s is expected to be less dependent on saturation and discrepancies in Fig. 8a and b can be reduced. But, results given in Fig. 8c and d show that the dimensionless s is still a strong function of saturation specially for main drainage and main imbibition curves. As it can be seen in Fig. 8c, there is not any significant difference among the curves obtained under main drainage. But the primary drainage curves are not identical, although the order of magnitude of them are the same. The difference among the curves is even larger under main imbibition (Fig. 8d). This difference is probably due to the fact that saturation alone is not sufficient to characterize the spatial distribution of fluids in the pores. As mentioned earlier, capillary processes and related coefficients (and thus s coefficient) may not only depend on saturation. As pointed out by Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 12 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx O’Carroll et al. (2010), the magnitude of the observed dynamic effect may be a function of fluid/fluid/solid contact line friction in addition to fluid viscosities. It is known that contact lines are the most detailed description of two-phase flow invasion and its fluid–fluid spatial distribution that include all pore-scale fluid– fluid properties such as equilibrium contact angle, interfacial tension, displacement distance, and front topology. These pore-scale issues are too difficult to be measured and simulated. As a macroscopic property, specific interfacial area may be an effective parameter in addition to effective viscosity and saturation. This dependence, however, is not being investigated in this paper, as it requires major computational effort. That is, in addition to main drainage and imbibition simulations carried out here, one has to perform a large number of (primary and secondary) scanning drainage and imbibition simulations to create capillary pressure– saturation–interfacial area relationship. Such a computational work is beyond the scope of the present work. Furthermore, as explained in Eq. (34), effective viscosity defined as linearly saturation-weighted viscosity does not account for all nonlinearities related to the fluid distribution as well as macroscopic interface front; the fluids can be distributed as disconnected blobs and macroscopic interface represent fingering or stable front. It is interesting to look at the variation of s versus two system parameters, namely the invading phase saturation (Sinv) and dimensionless effective viscosity (l), which is defined as the effective viscosity divided by the viscosity of invading fluid. Fig. 9 shows the contour plot of variation of log(s/leff) versus saturation and l. The contour plot is obtained by fitting a polynomial function to the data points of Fig. 8 resulted from drainage and imbibition simulations. The variation ranges for fitted parameters have been calculated for 95% confidence level based on t-distribution. The trajectories of the primary drainage and main imbibition simulations for different viscosity ratios have been also shown. Since the effective viscosity is divided by invading fluid’s viscosity, starting points (S) for M = 0.1 and M = 10 are different. But endpoints (E) are the same, almost equal to unity. As it is known, changing the invasion process from drainage to imbibition, may change the stability. This figure can be interpreted in analogy to the bifurcation diagrams. The upper branch is the unstable one, which is the tra- jectory  of drainage for M = 0.1 (S0.1 ? E) and imbibition for  M = 10 S010 ! E . The lower branch is the stable one, which is the trajectory for M = 10 (S10 ? E) and imbibition for  of drainage  M = 0.1 S00:1 ! E . The stable and unstable branches are separated from each other by the line associated with the trajectory of M = 1. Fig. 8a shows that during stable primary drainage (M P 1), for a wide range of saturation, non-equilibrium capillarity coefficient (s) increases with the decrease of wetting fluid saturation, which is similar to most findings as reported in Table 3. For M < 1, the reverse trend has been observed; namely s decreases with the decrease of wetting fluid saturation. An empirical formula relating s to medium and fluid properties was proposed by Stauffer (1978). In that formula, however, no dependence on saturation was included. We propose to modify Stauffer formula as follows: s¼  aleff Pcd kk qg 2 ð35Þ where l is replaced by leff, which is a function of saturation. In this equation, a is a constant,  is porosity, k and Pcd are the coefficients in Brooks-Corey formula (Brooks and Corey, 1964), k is the saturated permeability, q is the water mass density, and g is the gravity. Eq. (35) suggests that the change of s with saturation is proportional to the change of leff with saturation. Recalling the simple linear algebraic equation for leff, we can write: @ leff @Sw ¼ lw  ln ð36Þ   @l which means that @S@ sw / @Seff is negative for M > 1 and positive for w M < 1: @s / @Sw lw  ln 6 0; M P 1 lw  ln > 0; M < 1 This trend is in full agreement with the observations reported in the literature, as shown in Table 3. It was found in Mirzaei and Das (2007) in column-scale drainage simulations, in Joekar-Niasar et al. (2010a) in pore-network model simulations and in Bottero and Hassanizadeh (2006) in PCE–water experiments (M  0.9) that s increased with the decrease of wetting fluid saturation. For 10 4.5 4.75 5 5.25 5.5 5.75 6 6.25 above S0.1 , S' 10 unstable branch 8 6 4 2 S1.0 , S' 1.0 E S10 ,S' 0.1 0 0.0 stable branch 0.2 0.4 0.6 0.8 1.0 Saturation (invading) Fig. 9. Contour plot of log(s/leff) versus l and saturation of the invading phase obtained from dynamic simulations. l is defined as the ratio of leff to the viscosity of the invading phase. The arrows show the trajectory of drainage simulations from starting point, S, to end point, E, for different viscosity ratios. The corresponding starting points for imbibition are denoted by S0 . A simple polynomial surface fitted to the drainage and imbibition data points has the following form: 2 s=leff ¼ 6:945  1:526E5 þ 1:623  4:592E5Sinv  5:500  1:906E5l þ 3:591  1:615E5Sinv l  3:494  1:177E5Sinv þ 5:117  1:804E4l2 . Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 13 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx M = 0.1, the inverse trend of s as a function of saturation has been reported in pore-network modeling results of Joekar-Niasar et al. (2010a) and in air-water experiments by Camps-Roach et al. (2010). Main drainage and imbibition curves shown in Fig. 8a and b show a strong non-monotonic behavior in variation of s versus saturation. This might be due to contribution of spatial distribution of disconnected phase over the domain to the average phase pressures. While the disconnected phase does not contribute to the change of saturation, it does contribute to the averaged fluid pressure. The contribution of the trapped phase is obviously more important under imbibition than drainage. This can be more visibly seen in main imbibition curves shown in Fig. 8b. As it can be seen for all viscosity ratios, there is a sharp increase of s with increase of wetting phase saturation and afterwards s value decreases with the increase of wetting phase saturation and there is not a dominant trend over the whole range of saturation. 5. Summary and conclusion We have developed a DYnamic POre-network SImulator for Two-phase flow (DYPOSIT) to investigate non-equilibrium effects in fluid pressure fields under dynamic imbibition and drainage conditions. The pore-network model consists of truncated octahedron pore bodies and parallelepiped pore throats. The angularity in cross-sections allows for simultaneous flow of both nonwetting and wetting fluids. Two pressure fields, one for each fluid, are computed separately. Thus, counter-current flow within pore throats can be simulated. This means that capillary diffusion in the network is properly taken into account. To improve numerical stability of the model under capillary-dominated flow, a semi-implicit algorithm is employed. This allows us to simulate flow dynamics for different flow regimes and viscosity ratios for drainage as well as imbibition. We find that in spite of the fact that we have not included localscale dynamics of the interface (such as dynamic contact angle), the average fluids pressure difference (under large pressure gradients) is significantly distinct from the capillary pressure obtained under equilibrium conditions. Our analysis shows that the invasion regime has a significant effect on deviation of fluids pressure difference from capillary pressure. This deviation can be much larger for stable front regime compared to the viscous fingering regime. It means that in capillary fingering regime or viscous fingering regime less dynamic effect is resulted compared with stable front regime. Several numerical simulations for primary drainage, main imbibition, and main drainage are implemented. We have calculated non-equilibrium capillarity coefficient (s) as a function of saturation, for three different viscosity ratios under Dirichlet boundary conditions. We found that there is hysteresis in s–Sw relationship; i.e. although the order of magnitude of s during drainage and imbibition for a given viscosity ratio does not change significantly, the curves for primary drainage and main imbibition are different. This is probably due to the fact that saturation alone is not sufficient to characterize the distribution of fluids in the pores. As mentioned earlier, capillary processes and related coefficients (and thus s coefficient) may not only depend on saturation but also on specific interfacial area. This dependence, however, is not investigated in this paper, as it requires major computational effort. That is, in addition to main drainage and imbibition simulations carried out here, one has to perform a large number of (primary and secondary) scanning drainage and imbibition simulations. Such a computational work is beyond the scope of the present work. Our analysis also shows a strong dependence of s on viscosity ratio as well as effective viscosity. Variation of s with saturation is strongly dependent on the viscosity ratio. This is in agreement with laboratory experiments reported in Table 3. Given the definition of viscosity ratio as the ratio of nonwetting fluid viscosity to the wetting fluid’s, if it is smaller than one under drainage or larger than one under imbibition can create unstable invading fronts. Since different viscosity ratios lead to different distribution of phases, interfacial area between the two phases will be different. Thus, the dependence of s on fluid viscosities will be probably also modeled by the inclusion of specific interfacial area. Moreover, we have shown that snap-off is highly related to the viscosity ratio and capillary number. A viscous fingering during imbibition under Dirichlet boundary conditions creates non-monotonic distribution of trapped nonwetting fluid. With the invasion of wetting fluid during imbibition, flooding efficiency increases and less nonwetting fluid remains in the domain. Acknowledgments Comments by one of the anonymous reviewers were most valuable and helped to improve the manuscript. The authors are members of the International Research Training Group NUPUS, financed by the German Research Foundation (DFG) and The Netherlands Organization for Scientific Research (NWO). This research is financed by Utrecht Centre of Geosciences. Appendix A. Local pci —sw i relationship for drainage Consider the octahedron shown in Fig. 1 to be associated with pore body i. As mentioned earlier, the size of a pore body is specified by the radius of inscribed sphere of the p octahedron, Ri. So ffiffiffi the length of a typical edge O1A is equal to ai ¼pffiffi 6Ri . Total pffiffiffi volume of an an octahedron is then given by V oct:hedi ¼ 32 a3i ¼ 4 3R3i . Next, consider the corners to be truncated by pore throats connecting to the pore body. A pore throat ij results in a pyramid with square base and edge length 2rij to be cut off the octahedron (rij is the radius of inscribed circle p offfiffiffithe pore throat). The height of the pyramid will be equal to 2rij . Thus, the volume of the truncated pffiffi section is equal to 4 3 2 r3ij . So, the volume of the truncated octahedron forming pore body i is given by: pffiffiffi pffiffiffi 4 2X V i ¼ 4 3R3i  r3 j2Ni ij 3 ðA:1Þ During primary drainage, the wetting fluid saturation of a pore body specified below. An important varies from 1 to residual value Smin i intermediate saturation corresponds to the situation that the nonwetting fluid fills up the inscribed sphere of the pore body. Denoted by sdr i , the corresponding wetting fluid saturation is given by: sdr i ¼ 1 4pR3i p R3 ¼ 1  pffiffiffi 3 piffiffiffiP 3V i 3 3Ri  2 j2N r 3ij i ðA:2Þ For the derivation of local capillary pressure–saturation curves, we identify two different ranges as described below. w A.1. Capillary pressure for the range sdr i 6 si 6 1 (Zone 1) We assume that the entrance of nonwetting fluid into a pore body and the position of fluid–fluid interface may be idealized by sequences shown in Fig. A.1a. Position 1 shows the moment of invasion of nonwetting fluid into the pore body. The interface has the same radius of curvature as when it was in the pore throat ij and the local capillary pressure is the same as the (entry) capillary pressure of the pore throat, denoted by pce;ij (the starting point w of pci —sw i curve in Fig. A.2a at si ¼ 1). We assume that as the interface moves into the pore body, its radius remains unchanged until position 2 is reached. At this point, the saturation of nonwetting Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 14 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx r 3ij 2 3 r3 p Viji . Thus, for the range 2 1 P sw i P 1  3 p V i , the local capillary pressure will be constant equal to the entry capillary pressure of the pore throat; pci ¼ pce;ij (see the plateau in Fig. A.2a). As the interface moves into the pore body, it will expand and its capillary pressure decreases (interface 3 in Fig. A.1a). To simplify the geometry of the interface, we assume that from this saturation onwards, the nonwetting fluid fills up a sphere of increasingly larger radius until it fills up the inscribed sphere of the pore body, at which point the wetting fluid (a) ðA:3Þ r ci In this range, as increases, the local capillary pressure, given by pci ¼ 2rnw =r ci decreases until sdr i is reached (see Fig. A.2). The irregular shape of the resulting curves is approximated by the following formula for the range 1 to sdr i for sake of computational efficiency. sw i ¼ sdr i  þ 1 sdr i  1  sdr i  1  sdr i  4prci 3 3V i ðA:4Þ 4pr 3ij 3V i 2 1.5 1 Si dr 0 0 0.2 0.4 0.6 0.8 1 Saturation (Sw) (b) Imbibition_Analytical Imbibition_Fitted 3 2.5 Mean Curvature ( ) 4pr ci 3 3V i Drainage_Analytical Drainage_Fitted 0.5 saturation is sdr i given by Eq. (A.2). If the radius of such an intermediate sphere is r ci , then the wetting fluid saturation in this range is given by the following relation: sw i ¼ 1 3 2.5 Mean Curvature ( ) fluid in the pore body is equal to 2 1.5 1 0.5 A.2. Capillary pressure for the range 0 < sw i 6 sdr i (Zone 2) Once the nonwetting fluid fills the inscribed sphere, and after that, the wetting fluid will be present only along the edges of the pore body. So, interfaces will be formed along pore body’s 12 edges as well as in the opening of pore throats not invaded by the nonwetting fluid yet. The latter interfaces are formed in the vertices of the truncated octahedron (where non-invaded pore throats are connected). Edge and vertex interfaces are shown in Fig. A.1b and d, respectively. Edge interfaces form part of a cylinder (with only one finite radius of curvature), and vertex interfaces form part of a spherical surface. When all pore throats have been invaded by the nonwetting fluid, there will be no vertex interface in a pore body. Si imb 0 0 0.2 0.4 0.6 0.8 1 w Saturation (S ) Fig. A.2. Local pci —sw i curves during drainage (a) and imbibition (b). Thick curves show the analytical calculations and thin curves show the fitted curves. Let the radius of edge interfaces be denoted by rci . We assume that the radius of curvature of interfaces along the edges is negligible compared r ci . As shown in Fig. A.1b, we denote the half corner angle of all edges by b and half angle of interfaces by a. We can show that NMO = b = 0.5cos1(1/3) = 0.9553, and 0 0 NOM = a = (0.61548  h). For given b and rci , area of the MNM N (a) (b) (c) (d) Fig. A.1. Schematic presentation of interfaces in a pore body (a) expansion of an interface into a pore body before filling the inscribed sphere (b and c) interface geometry along the edges (d) interface geometry in a vertex of pore body when the corresponding pore throat is filled with wetting fluid. Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx is equal to 2  ðAMMNO  A follows: AMMNO ¼ 0:5r ci A ONM0 2 ONM0 Þ, which are themselves defined as B.1. All pore throats are (partially) filled by the nonwetting fluid, I = Z (Zone 1) ðA:5Þ In this range all pore throats are (partially) filled by the nonwetting fluid and pci —sw relationship is the same as Zone 2 i   dr for drainage process. Thus, pci —sw 0 < sw i 6 si i curve follows the same curve shown in Fig. A.2(a) till one of the pore throats is fully filled by the wetting fluid. sin a cos h sin b 2 ¼ 0:5r ci a ðA:6Þ Consequently, the cross-sectional area of residual wetting fluid may be written as: Awet ¼ r ci i 2 15   sin a cos h a sin b ðA:7Þ Total length of edges filled with the wetting fluid is given by: X pffiffiffi Li ¼ 12 6Ri  8 r ij j2Ni ðA:8Þ With the total pore body volume given by Eq. (A.1), the saturation of wetting fluid in the pore body is determined as:   pffiffiffi  pffiffiffi P a cos h  a 3 2rci 2 3 6Ri  2 j2Ni r ij sinsin b pffiffiffi sw P i ¼ 3 6R3i  2 j2Ni r 3ij ðA:9Þ This relation is valid only when all pore throats connected to the pore body i are invaded by the nonwetting fluid, so that there is no vertex interface in that pore body. This may not be the case when the nonwetting   fluid has filleddronly the inscribed sphere of the pore body r ci ¼ Ri for which sw i ¼ si given by Eq. (A.2). In fact, substituting r ci ¼ Ri in Eq. (A.9), results in a saturation different from sdr i . We denote this saturation by si . Forcing saturation to vary between sdr i and smin (the latter one was defined by Eq. (22)), results in the foli lowing equation: min sdr i  si si  smin i   2 pffiffiffi pffiffiffi 3 P a cos h  a 3 2 3 6Ri  2 j2Ni r ij sinsin b 5  pffiffiffi   smin 4 P i 4j2i 3 6R3i  2 j2Ni r 3ij min sw þ i ¼ si B.2. More than one of the pore throats and not all of them are (partially) filled by the nonwetting fluid, 1 < I < Z (Zone 2) As long as all the pore throats are (partially) filled by the nonwetting fluid, the wetting fluid will be present only along the 12 edges of the pore body. When nonwetting fluid recedes completely from one of the pore throats, an interface will be created in the opening of that pore throat. This vertex interfaces is shown in Fig. A.1d. For simplicity, we assume that all wetting fluid’s volume existing in that pore body will remain in that vertex. If nonwetting fluid recedes from two or more pore throats, the wetting fluid will be distributed among them. Knowing the volume of the wetting fluid and the number of pore throats nonwetting fluid has receded from, radius of the curvature of the interface can be estimated. Re0 0 0 0 fer to Fig. A.1d, let’s denote the length of edge of pyramid O1A B C D 0 0 base by a . We can calculate a based on the geometry of pore body, pore throats, and the number of pore throats partially filled with nonwetting fluid (I), as follows: a0 ¼ " pffiffiffi #1=3 6I 3 2 w 8 X r3ij si V i þ 6I 6  I j¼1 0 in which Vi can be calculated from Eq. (A.1). Using the value of a , mean radius of curvature of the interface is given by 0  a0 pffi . As it can be seen with decrease of I, a rci ¼ will 2 cos hþsin1 ðA:10Þ This relationship is shown as a thick curve in Fig. A.2a. A.3. Local pci —sw i relationship for the full range of saturation Finally, to reduce the computational efforts, we have chosen to fit one single curve to pci for the full range of saturation. Eq. (21) shows the fitted pci —sw i relationship for the full range of saturation in a pore body. This relationship is shown by a thin curve in Fig. A.2a. Appendix B. Local pci —sw i relationship for imbibition Local variation of capillary pressure is much more complex under imbibition than drainage. As mentioned in Section 2.3.4, interface topology, and consequently capillary pressure in a pore body are controlled by the number of pore throats which are still (partially) filled by the nonwetting fluid. Introducing Z as the total number of pore throats connected a pore body, and I as the number of pore throats (partially filled by the nonwetting fluid), the following zones are introduced to defining local pci —sw i relationship:  All pore throats are (partially) filled by the nonwetting fluid, I = Z (Zone 1).  More than one of the pore throats and not all of them are (partially) filled by the nonwetting fluid, 1 < I < Z (Zone 2).  Only one pore throat is (partially) filled by the nonwetting fluid, I = 1 (Zone 3). ðB:1Þ 3 3 decrease and consequently capillary pressure will increase. This trend is shown schematically in Fig. A.2b as a step-wise curve. The jumps in the curve occur when the nonwetting fluid recedes from one of the pore throats (I decreases), and refer to Eq. (B.1), radius of curvature decreases. We follow this trend till only one of the pore throats remains filled with the nonwetting fluid. There is a qualitative consistency between the trend of this local pci —sw i curve with the experimental observations of Lenormand and Zarcone (1984), who observed with decrease of number of pore throats (partially) filled with the nonwetting fluid, local capillary pressure may increase. B.3. Only one pore throat is (partially) filled by the nonwetting fluid, I = 1 (Zone 3) When only one of the pore throats is (partially) filled with the nonwetting fluid, the behavior of the interface follows the explanations given in Zone I during drainage. B.4. Local pci —sw i relationship for the full range of saturation As mentioned before, the variation of capillary pressure with saturation is very irregular and it depends on the fluid occupancy of pore throats. The thick curve in Fig. A.2b shows an example of this curve for a pore body. But, handling such a discontinuous curve in the pore-network model is nontrivial. Thus, to reduce the computational efforts, we have chosen to fit one single curve to pci for the full range of saturation. Eq. (23) shows the fitted pci —sw i relationship for the full range of saturation in a pore body. This relationship is shown by a thin curve in Fig. A.2b. The fitted curve results decrease of capillary pressure with increase of Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 16 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx saturation to a minimum capillary pressure. This minimum capillary pressure occurs at a saturation corresponding to the start of Zone 3 during imbibition, where only one pore throat is (partially) filled by the nonwetting fluid. We denote this saturation by simb as i shown in Fig. A.2b. References Aker, E., Maloy, J., Hansen, A., Batrouni, G.G., 1998a. A two-dimensional network simulator for two-phase flow in porous media. Transport Porous Media 32, 163–186. Aker, E.K., Maloy, K.J., Hansen, A., 1998b. Simulating temporal evolution of pressure in two-phase flow in porous media. Phys. Rev. E 58, 2217–2226. Al-Gharbi, M.S., Blunt, M.J., 2005. Dynamic network modeling of two-phase drainage in porous media. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 71, 016308. Avraam, D.G., Payatakes, A.C., 1995a. Generalized relative permeability coefficients during steady-state, two-phase flow in porous media and correlation with the flow mechanisms. Transport Porous Media 20, 135–168. Azzam, M.I.S., Dullien, F.A.L., 1977. Flow in tubes with periodic step changes in diameter: a numerical solution. Chem. Eng. Sci. 32, 1445–1455. Bakke, S., Øren, P.-E., 1997. 3-d pore-scale modelling of sandstones and flow simulations in the pore networks. SPE J. 2, 136–149. Barenblatt, G., Patzek, T., Silin, D., 2003. The mathematical model of nonequilibrium effects in water–oil displacement. SPE J. 8, 409–416. Berentsen, C.W.J., Hassanizadeh, S.M., 2006. On the equivalence of measured and averaged pressures in porous media. In: 16th International Conference on Computational Methods in Water Resources. Bottero, S., 2009. Advances in Theories of Capillarity in Porous Media. Ph.D. Thesis, Utrecht University. Bottero, S., Hassanizadeh, S.M., 2006. Experimental and numerical study to investigate dynamic capillary pressure effect in two-phase flow in porous media. In: 16th International Conference on Computational Methods in Water Resources. Brooks, R.H., Corey, A.T., 1964. Hydraulic Properties of Porous Media. Tech. Rep. Hydrol. Pap. No. 3, Colo. State Univ. Camps-Roach, G., O’Carroll, D.M., Newson, T.A., Sakaki, T., Illangasekare, T.H., 2010. Experimental investigation of dynamic effects in capillary pressure: grain size dependency and upscaling. Water Resour. Res. 46, W08544. Chen, D.Q., Pyrak-Nolte, L.J., Griffin, J., Giordano, N.J., 2007. Measurement of interfacial area per volume for drainage and imbibition. Water Resour. Res. 43, W12504. Cheng, J.T., Pyrak-Nolte, L.J., Nolte, D.D., 2004. Linking pressure and saturation through interfacial area in porous media. Geophys. Res. Lett. 31, L08502. Constantinides, G.N., Payatakes, A.C., 1996. Network simulation of steady-state twophase flow in consolidated porous media. AIChE J. 42, 369–382. Dahle, H.K., Celia, M.A., 1999. A dynamic network model for two-phase immiscible flow. Comput. Geosci. 3, 1–22. Dahle, H.K., Celia, M.A., Hassanizadeh, S.M., 2005. Bundle-of-tubes model for calculating dynamic effects in the capillary–pressure–saturation relationship. Transport Porous Media 58, 5–22. Das, D., Gauldie, R., Mirzaei, M., 2007. Dynamic effects for two-phase flow in porous media: fluid property effects. AIChE J. 53, 2505–2520. Das, D.B., Mirzaei, M., Widdows, N., 2006. Non-uniqueness in capillary pressure– saturation–relative permeability relationships for two-phase flow in porous media: interplay between intensity and distribution of random microheterogeneities. Chem. Eng. Sci. 61, 6786–6803. Dias, M.M., Payatakes, A.C., 1986a. Network models for two-phase flow in porous media part 1. Immiscible microdisplacement of non-wetting fluids. J. Fluid Mech. 164, 305–336. Dias, M.M., Payatakes, A.C., 1986b. Network models for two-phase flow in porous media part 2. Motion of oil ganglia. J. Fluid Mech. 164, 337–358. Entov, V.M., 1980. Theory of nonequilibrium effects associated with the flow of nonuniform fluids in porous media. Fluid Dynam. 15, 365–369. Fayers, F.J., Blunt, M.J., Christie, M.A., 1990. Accurate calibration of empirical viscous fingering models. In: 2nd European Conference on the Mathematics of Oil Recovery, 75015 Paris, pp. 45–55. Fong, K.W., Jefferson, T.H., Suyehiro, T., Walton, L., 1993. Guide to the slatec common mathematical library. <http://www.netlib.org/slatec>. Gielen, T., Hassanizadeh, S.M., Leijnse, A., Nordhaug, H.F., 2005. Dynamic effects in multiphase flow: a pore-scale network approach. In: Das, D.B.M.H.S. (Ed.), Upscaling Multiphase Flow in Porous Media. Springer, The Netherlands, pp. 217–236. Gielen, T.W.J., 2007. Upscaling Multiphase Transport Processes in Porous Media. Ph.D. Thesis, Delft University of Technology. Hassanizadeh, S., Celia, M., Dahle, H., 2002. Dynamic effects in the capillary pressure–saturation relationship and their impacts on unsaturated flow. Vadose Zone J. 1, 38–57. Hassanizadeh, S.M., Gray, W.G., 1990. Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries. Adv. Water Resour. 13, 169–186. Hassanizadeh, S.M., Gray, W.G., 1993a. Thermodynamic basis of capillary pressure in porous media. Water Resour. Res. 29, 3389–3405. Hassanizadeh, S.M., Gray, W.G., 1993b. Toward an improved description of the physics of two-phase flow. Adv. Water Resour. 16, 53–67. Hassanizadeh, S.M., Oung, O., Manthey, S., 2004. Laboratory experiments and simulations on the significance of non-equilibrium effect in the capillary pressure–saturation relationship. In: Schanz, T. (Ed.), Unsaturated Soils: Experimental Studies. Proceedings of the International Conference: From Experimental Evidence towards Numerical Modeling of Unsaturated Soils, vol. 1. Springer Verlag, Berlin, pp. 3–14. Held, R.J., Celia, M.A., 2001. Modeling support of functional relationships between capillary pressure, saturation, interfacial area and common lines. Adv. Water Resour. 24, 325–343. Hughes, R.G., Blunt, M.J., 2000. Pore scale modeling of rate effect in imbibition. Transport Porous Media 40, 295–322. Joekar-Niasar, V., Hassanizadeh, S., Leijnse, A., 2008. Insights into the relationships among capillary pressure, saturation, interfacial area and relative permeability using pore-network modeling. Transport Porous Media 74, 201–219. Joekar-Niasar, V., Hassanizadeh, S.M., Dahle, H.K., 2010a. Non-equilibrium effects in capillarity and interfacial area in two-phase flow: dynamic pore-network modelling. J. Fluid Mech. 655, 38–71. Joekar-Niasar, V., Hassanizadeh, S.M., Pyrak-Nolte, L.J., Berentsen, C., 2009. Simulating drainage and imbibition experiments in a high-porosity micromodel using an unstructured pore network model. Water Resour. Res. 45, W02430. Joekar-Niasar, V., Prodanović, M., Wildenschild, D., Hassanizadeh, S.M., 2010b. Network model investigation of interfacial area, capillary pressure and saturation relationships in granular porous media. Water Resour. Res. 46, W06526. Kalaydjian, F., Marle, C.M., 1987. Thermodynamic aspects of multiphase flow in porous media. Collection colloques et séminaires – Institut français du pétrole 45, 513–531. Koplik, J., Lasseter, T.J., 1985. Two-phase flow in random network models of porous media. Soc. Petrol. Eng. J. 25, 89–110. Korteland, S., Bottero, S., Hassanizadeh, S.M., Berentsen, C.W.J., 2010. What is the correct definition of macroscale pressure? Transport Porous Media 84, 153– 175. Koval, E.J., 1963. A method for predicting the performance of unstable miscible displacements in heterogenous media. Trans. AIME 219, 145–150. Lenormand, R., Zarcone, C., 1983. Mechanism of the displacement of one fluid by another in a network of capillary ducts. J. Fluid Mech. 135, 337–353. Lenormand, R., Zarcone, C., 1984. Role of roughness and edges during imbibition in square capillaries. In: 59th Annual Technical Conference of the SPE. Houston, Richardson, TX, No. 13264. Lewalle, J., Singh, K., Bambacht, J., 1994. Analysis of the wet pressing of paper pulp. Int. J. Multiphase Flow 20, 415–437. Ma, S., Mason, G., Morrow, N.R., 1996. Effect of contact angle on drainage and imbibition in regular polygonal tubes. Colloids Surfaces A 117, 273– 291. Manthey, S., Hassanizadeh, S.M., Helmig, R., 2005. Macro-scale dynamic effects in homogeneous and heterogeneous porous media. Transport Porous Media 58, 121–145. Manthey, S., Hassanizadeh, S.M., Oung, O., Helmig, R., 2004. Dynamic capillary pressure effects in two-phase flow through heterogeneous porous media. In: 15th International Conference on Computational Methods in Water Resources. Mason, G., Morrow, N.R., 1987. Meniscus configurations and curvatures in nonaxisymmetric pores of open and closed uniform cross section. Proc. Roy. Soc. Lond. Series A, Math. Phys. Sci. 414, 111–133. http://www.jstor.org/stable/ 2398111. Mayer, R.P., Stowe, R.A., 1965. Mercury porosimetry-breakthrough pressure for penetration between packed spheres. J. Colloid Sci. 20, 891–911. Mirzaei, M., Das, D.B., 2007. Dynamic effects in capillary pressure–saturations relationships for two-phase flow in 3d porous media: implications of microheterogeneities. Chem. Eng. Sci. 62, 1927–1947. Mogensen, K., Stenby, E.H., 1998. A dynamic two-phase pore-scale model for imbibition. Transport Porous Media 32, 299–327. Nguyen, V.H., Sheppard, A.P., Knackstedt, M.A., Pinczewski, W., 2006. The effect of displacement rate on imbibition relative permeability and residual saturation. J. Petrol. Sci. Eng. 52, 54–70. Nordbotten, J., Celia, M., Dahle, H., Hassanizadeh, S., 2007. Interpretation of macroscale variables in Darcy’s law. Water Resour. Res. 43, W08430. Nordbotten, J., Celia, M., Dahle, H., Hassanizadeh, S., 2008. On the definition of macroscale pressure for multiphase flow in porous media. Water Resour. Res. 44, W06S02. Nordhaug, H.F., Celia, M., Dahle, H.K., 2003. A pore network model for calculation of interfacial velocities. Adv. Water Resour. 26, 1061–1074. O’Carroll, D.M., Phelan, T.J., Abriola, L.M., 2005. Exploring dynamic effects in capillary pressure in multistep outflow experiments. Water Resour. Res. 41, W11419. O’Carroll, D.M., Mumford, K.G., Abriola, L.M., Gerhard, J.I., 2010. Influence of wettability variations on dynamic effects in capillary pressure. Water Resour. Res. 46, W08505. Øren, P.E., Bakke, S., Arntzen, O.J., 1998. Extending predictive capabilities to network models. SPE J. 3, 324–336. Oung, O., Hassanizadeh, S.M., Bezuijen, A., 2005. Two-phase flow experiments in a geocentrifuge and the significance of dynamic capillary pressure effect. J. Porous Media 8, 247–257. Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007 V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx Pereira, G.G., Pinczewski, W.V., Chan, D.Y.C., Paterson, L., Øren, P.E., 1996. Pore-scale network model for drainage-dominated three-phase flow in porous media. Transport Porous Media 24, 167–201. Porter, M.L., Schaap, M.G., Wildenschild, D., 2009. Lattice-boltzmann simulations of the capillary pressure–saturation–interfacial area relationship for porous media. Adv. Water Resour. 32, 1632–1640. Prat, M., 2002. Recent advances in pore-scale models for drying of porous media. Chem. Eng. J. 86, 153–164. Princen, H.M.J., 1969a. Capillary phenomena in assemblies of parallel cylinders i. capillary rise between two cylinders. J. Colloid Interface Sci. 30, 69–75. Princen, H.M.J., 1969b. Capillary phenomena in assemblies of parallel cylinders ii. capillary rise in systems with more than two cylinders. J. Colloid Interface Sci. 30, 359–371. Princen, H.M.J., 1970. Capillary phenomena in assemblies of parallel cylinders iii. liquid columns between horizontal parallel cylinders. J. Colloid Interface Sci. 34, 171–184. Pyrak-Nolte, L.J., 2007. Url of Website with Original Data. <http://www.physics. purdue.edu/rockphys/dataimages/>. <http://www.physics.purdue.edu/rockphys/ DataImages/>. Ransohoff, T.C., Radke, C.J., 1988. Laminar flow of a wetting liquid along the corners of a predominantly gas-occupied noncircular pore. J. Colloid Interface Sci. 121, 392–401. Raoof, A., Hassanizadeh, S., 2010. A new method for generating pore-network models of porous media. Transport Porous Media 81, 391–407. 17 Reeves, P.C., Celia, M.A., 1996. A functional relationship between capillary pressure, saturation, and interfacial area as revealed by a pore-scale network model. Water Resour. Res. 32, 2345–2358. Singh, M., Mohanty, K.K., 2003. Dynamic modelling of drainage through threedimensional porous materials. Chem. Eng. Sci. 58, 1–18. Stauffer, F., 1978. Time dependence of the relations between capillary pressure, water content and conductivity during drainage of porous media. IAHR. Thessaloniki, pp. 3.35–3.52. Thompson, K.E., 2002. Pore-scale modeling of fluid transport in disordered fibrous materials. AIChE J. 48, 1369–1389. Van der Marck, S.C., Matsuura, T., Glas, J., 1997. Viscous and capillary pressures during drainage: network simulations and experiments. Phys. Rev. E 56, 5675–5687. Vidales, A.M., Riccardo, J.L., Zgrabli, G., 1998. Pore-level modelling of wetting on correlated porous media. J. Phys. D – Appl. Phys. 31, 2861–2868. Vizika, O., Avraam, D.G., Payatakes, A.C., 1994. On the role of viscosity ratio during low-capillary-number forced imbibition in porous media. J. Colloid Interface Sci. 165, 386–401. Whitaker, S., 1977. Simultaneous heat, mass, and momentum transfer in porous media: a theory of drying. Adv. Heat Transfer 13, 119–200. Yang, D., Currier, R.P., Zhang, D.Z., 2009. Ensemble phase averaged equations for multiphase flows in porous media. Part 1: the bundle-of-tubes model. Int. J. Multiphase Flow 35, 628–639. Zhou, D., Blunt, M.J., Orr, F.M., 1997. Hydrocarbon drainage along corners of noncircular capillaries. J. Colloid Interface Sci. 187, 11–21. Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007