Uniqueness of Specific Interfacial Area–Capillary Pressure–Saturation Relationship Under Non-Equilibrium Conditions in Two-Phase Porous Media Flow

The capillary pressure–saturation (Pc–Sw) relationship is one of the central constitutive relationships used in two-phase flow simulations. There are two major concerns regarding this relation. These concerns are partially studied in a hypothetical porous medium using a dynamic pore-network model called DYPOSIT, which has been employed and extended for this study: (a) Pc–Sw relationship is measured empirically under equilibrium conditions. It is then used in Darcy-based simulations for all dynamic conditions. This is only valid if there is a guarantee that this relationship is unique for a given flow process (drainage or imbibition) independent of dynamic conditions; (b) It is also known that Pc–Sw relationship is flow process dependent. Depending on drainage and imbibition, different curves can be achieved, which are referred to as “hysteresis”. A thermodynamically derived theory (Hassanizadeh and Gray, Water Resour Res 29: 3389–3904, 1993a) suggests that, by introducing a new state variable, called the specific interfacial area (anw, defined as the ratio of fluid–fluid interfacial area to the total volume of the domain), it is possible to define a unique relation between capillary pressure, saturation, and interfacial area. This study investigates these two aspects of capillary pressure–saturation relationship using a dynamic pore-network model. The simulation results imply that Pc–Sw relation not only depends on flow process (drainage and imbibition) but also on dynamic conditions for a given flow process. Moreover, this study attempts to obtain the first preliminary insights into the global functionality of capillary pressure–saturation–interfacial area relationship under equilibrium and non-equilibrium conditions and the uniqueness of Pc–Sw–anw relationship.


List of symbols
Half-corner angle in a pore throat (rad) Γ A geometrical coefficient for calculating the corner conductance κ i Curvature of the interface in the pore body i (m −1 ) σ nw Interfacial tension (N/m) π 3.1415 θ Contact angle (rad) Resistance factor μ α Viscosity of phase α (Pa s) 1 Introduction

Specific Interfacial Area and its Role in Two-Phase Flow
Physically consistent modeling of two-phase flow in porous media remains a challenge, as the governing equations that are currently used do not fully account for the dynamics of the system. In particular, the capillary pressure-saturation function, which relates the pressures of the two fluids to a given saturation, is empirically measured under equilibrium conditions but applied to all dynamic conditions. It does not provide any information about the dynamics of the moving interfaces between the fluids. It is hysteretic and non-unique. A given system of fluids and porous medium may have many scanning drainage and imbibition curves, in addition to primary and main drainage and imbibition curves, depending on the initial condition and flow direction (drainage or imbibition). Using rational thermodynamics, Gray (1990, 1993b) proposed that the hysteresis in the P c -S w relationship can be modeled by introducing the specific interfacial area a nw , defined as the ratio of the total capillary interfacial area to the porous media volume. Based on their theory, the following relation has been suggested: which is conjectured to be valid under equilibrium and non-equilibrium conditions. The dynamics of fluid-fluid interfaces is not only important for a proper modeling of two-phase flow but also in all cases where mass or heat transfer among phases occurs, including evaporation in porous media, enhanced oil recovery methods, biodegradation of non-aqueous phases, and CO 2 dissolution in brine. For all these applications, it is important not only to estimate the capillary pressure-saturation relation but also the fate of interfaces under equilibrium and non-equilibrium conditions. In addition to a new capillary pressure-saturation relationship, Hassanizadeh and Gray (1993b) proposed an extended Darcy's law for two-phase flow, which includes gradient of specific interfacial area and gradient of saturation as well as new material coefficients. To close the system for these state variables, they propose a new equation that describes the balance of specific interfacial area. A numerical analysis of these parameters and their behavior in a simple drainage case has been provided in Joekar-Niasar and Hassanizadeh (2011b). However, non-equilibrium effects on P c -S w -a nw relation have yet to be studied. A number of computational and experimental works have shown that for a wide range of drainage and imbibition histories under quasi-static conditions, P c -S w -a nw surfaces more or less coincide. The authors' previous studies investigated the validity of this theory under equilibrium conditions. In the study of Joekar-Niasar et al. (2009), micromodel two-phase flow experiments were simulated using a medial-axis based pore-network model. The pore-network model results match well with the experimental data not only for P c -S w but also for a nw -S w curves. The simulation and experimental results show that the P c -S w -a nw surface obtained for drainage and imbibition match well under equilibrium conditions. The average error is around 0.17. In another study that compares experimental and pore-network model results for glass beads experiments, a similar conclusion was reached (Joekar-Niasar et al. 2010b). These results indicated that inclusion of a nw leads to the removal or significant reduction of hysteresis in P c -S w relationship (see e.g., Reeves and Celia 1996;Held and Celia 2001;Cheng et al. 2004;Chen et al. 2007;Joekar-Niasar et al. 2008;Porter et al. 2009). All these studies have a common characteristics: the analyses were conducted under equilibrium conditions in which the capillary interfaces do not move in mostly granular media. A major question is whether a relationship determined under equilibrium conditions is also applicable under non-equilibrium conditions. For example, even if a unique capillary pressure-saturation-interfacial area relationship exists under quasi-static conditions, does such a surface also exist under dynamic conditions? Answering such a question requires the measurement of capillary pressure, saturation, and specific interfacial area under transient flow conditions. While such variables can be measured under quasi-static or steady-state flow conditions, it is almost impossible to measure them under transient conditions. In 3D imaging (for example, X-ray tomography), the exact position of moving interfaces cannot be captured because of the time required for imaging versus the physical time scale of the interface displacements. Also, in 2D imaging techniques, it is only possible to visualize interfaces that are normal to the solid boundary (referred to as main terminal interfaces). Measuring local or average capillary pressure under transient flow conditions is even more difficult. Commonly, capillary pressure is measured as the difference between the non-wetting and wetting phase pressures, which is a correct way of measuring capillary pressure under static conditions, but not under dynamic conditions because of non-equilibrium capillarity effects (Hassanizadeh et al. 2002). The local capillary pressure of a meniscus could be estimated under flow conditions from its curvature, but, as explained above, this is difficult to determine experimentally. In the absence of experimental measurement methods, the pore-scale modeling may be the most feasible option for investigating the dynamics of capillary interfaces.

Pore-Network Modeling
Among the pore-scale and subpore-scale models (for example, pore-network model, Lattice Boltzmann, smoothed particle hydrodynamics, level set, direct simulation, density functional method, etc.), the pore-network models have the coarser computational discretization. In this approach, a porous medium is delineated by idealized interconnected geometries representing the pore space. Then, each geometry is usually considered as a computational node, while in other methods it can be discretized to hundreds or thousands of nodes. Because this study aimed to investigate a continuum-scale theory, it was necessary to simulate a sufficiently large porous medium. In addition to this advantage, pore-network modeling was applied successfully for many two-phase flow problems. Specifically, this study found good agreement between modeling and experimental results for the relationship of capillary pressure, saturation, and interfacial area under drainage and imbibition under equilibrium conditions (Joekar-Niasar et al. 2009, 2010b. Finally, given the idealization of the geometries of porous media and interfaces, their description is based on an infinite resolutions, which prevents errors introduced by the discretization of the geometries. Pore-network models can be divided into quasi-static and dynamic ones. In the first category, invasion percolation simulation using Young-Laplace equation is performed and only equilibrium status corresponding to the boundary pressures can be presented. (For an extensive review of its applications, see Blunt et al. 2002.) However, in dynamic pore-network models, the two-phase flow is simulated with time, which is much more computationally demanding. An extensive review can be found in Joekar-Niasar and Hassanizadeh (2010).
Here, this study employed the DYPOSIT model (DYnamic POre-network SImulator for Two-phase flow), developed and discussed in Joekar-Niasar et al. (2010a), Joekar-Niasar and Hassanizadeh (2011a). The model was improved to include mobilization of a disconnected single blob (isolated in a pore body), which was absent in the previous version. DYPOSIT model is a complete dynamic pore-network model, which can simulate drainage and imbibition. Either constant-pressure or constant-flux boundary conditions may be chosen. Major pore-scale invasion mechanisms, such as break-up, snap-off, pinch-off, coalescence, and mobilization of the discontinuous phase, have been included. This study aimed to investigate whether the capillary pressure-saturation-interfacial area relationship determined under quasi-static conditions is equally applicable under transient flow conditions. So, first the capillary pressure-saturation-interfacial area surface was obtained, based on results of a large number of quasi-static drainage and imbibition simulations, including many scanning curves. Next, non-equilibrium drainage and imbibition processes in a wide range of capillary numbers (10 −7 -10 −5 ) were simulated. Simulation results are presented for three different viscosity ratios (ratio of the viscosity of the invading fluid to the receding one) namely 1, 0.1, and 10. Capillary number (C a ) is traditionally defined as the ratio of viscous forces of the invading phase to capillary forces ( μ inv q inv σ nw ). Based on these results, a non-equilibrium capillary pressure-saturation-interfacial area surface was constructed. Then, equilibrium and nonequilibrium surfaces were compared to determine whether they were significantly different.

Model Description
Details about the structure and geometry of the model, as well as its governing equations are described in Joekar-Niasar and Hassanizadeh (2011a). For brevity, only a short description is provided here.

Structure and Geometry
The pore network is a regular three-dimensional lattice with a fixed coordination number of six. Pore bodies and pore throats are presented by "truncated octahedron" and "parallelepiped," respectively. This allows simultaneous existence of both phases 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. 1. Truncated sections of the octahedron have the shape of square pyramids with base width of 2r i j (which is equal to the size of the pore throat i j) as shown in Fig. 1. The size distribution of pore bodies (radius R i of the inscribed spheres) is specified by a truncated log-normal distribution, with no spatial correlation. The algorithm for generation of pore throats based on the inscribed radius of pore bodies is explained in detail in Joekar-Niasar and Hassanizadeh (2011a). The size of the pore network is chosen to be larger than one REV, which has been determined based on a sensitivity analysis of the capillary pressure-saturation and relative permeability curves for Fig. 1 Schematic presentation of a pore body and its connected pore throats. Truncated parts of the pore body are square pyramids with the base width of 2r i j . r i j is the inscribed radius in the cross-section of pore throat i j different sizes of networks with given statistical properties (Joekar-Niasar and Hassanizadeh 2011b). Given the size of REV equal to a cubic 35 pore bodies and considering five pore bodies for reducing the reservoir boundary effects, a network with 35 × 35 × 45 pore bodies along the flow direction is considered. Statistical properties of the distributions of pore bodies and pore throats inscribed spheres are given in Table 1.

General Pore-Scale Equations for Two-phase Flow
Given a pore channel i j, the volumetric flux of phase α through the channel in a laminar Newtonian flow is given by Washburn formula (Washburn 1921): where K α i j is the conductivity of the pore channel as a function of geometry and fluid occupancy. The positive flow direction is from pore body i to pore body j. Note that assigning different pressures and conductivities to each phase in a pore element, as opposed to one single pressure to a pore body and one effective conductivity to a pore throat (see e.g., Al-Gharbi and Blunt 2005; Mogensen and Stenby 1998) has major advantages. For example, mechanisms related to the local capillary pressure (such as snap-off, counter-current flow) can be included in simulations. Assuming incompressibility of the phases, the volumetric balance equation for the pore body i reads: where N i is the set of all pore throats connected to the pore body i. Equation 3 cannot be solved individually because p n i and p w i are both unknown. For a given fluid configuration, Eq. 3 is augmented by an equation for the pore-scale capillary pressure p c i , which is known to be a function of the local curvature of the interface For computational ease, the curvature is related to the local volume fraction of the phases (see Sect. 2.3). Thus, (4) can be rewritten as follows: where s α i denotes the volume fraction of phase α in the pore body i. The equation of volume balance for phase α within a pore body may be written as: where V i is the volume of the pore body i. Given that s w i + s n i = 1, summation of Eq. 6 results directly in Eq. 3.

Pressure Field Solver
In order to reduce the computational load, the governing equations are reformulated in terms of a total pressurep i in a pore body, defined as the volume fraction-weighted average of fluid pressures:p Equation 3 can be rewritten forp i : In this equation, the right-hand side, as well as the coefficients of the left-hand side, depend on local volume fraction only. This linear system of equations was solved forp i by the diagonally scaled biconjugate gradient method using the SLATEC mathematical library (Fong et al. 1993).

Saturation Update
After calculatingp i and based on s α i values, pressures of phases can be back-calculated. Commonly in the literature, Eq. 6 has been used to calculate the new s w i values explicitly. However, this procedure will result in numerical problems especially for capillary-dominated and viscous-fingering flow regimes (see for example, Koplik and Lasseter 1985). This study used a procedure analogous to fractional flow formulation to control the nonlinearities under such flow conditions. Details are described in Joekar-Niasar et al. (2010a). This results in a semi-implicit discretization of Eq. 6 as follows: ⎛ where superscript k denotes time step level, and is calculated from the upstream pore body. Note that, since Q tot i j and K α i j are calculated from time step k, this scheme is not fully implicit. Again, the diagonally scaled biconjugate gradient method from SLATEC mathematical library (Fong et al. 1993) was used to solve Eq. 9. The stability of the model for a wide range of flow conditions (from capillary-dominated regime to viscous-fingering regime) and a wide range of viscosity ratios was demonstrated in Joekar-Niasar et al. (2010a).

Time Stepping
The time step, denoted by t i , was determined based on the time of filling of pore bodies by the non-wetting phase or wetting phase, for drainage and imbibition simulations, respectively. The wetting phase volume fraction in a pore body varied between 1 and s w i,min as the assumption was that a pore body may be drained down to a minimum volume fraction (Joekar-Niasar and Hassanizadeh 2011a). On the other hand, if local imbibition occurs, the wetting phase volume fraction in a pore body can go back to 1. So, t i was calculated for all pore bodies, depending on the local process, from the following formulas: where, the accumulation rate of the non-wetting phase is defined as q n i = j∈N i Q n i j . Then, the global time step is chosen to be the minimum t i .
It should be noted that we have imposed a truncation criterion of 10 −6 for s w i − s w i,min and 1 − s w i to reduce the computational cost for solving the equations.

Local Rules
The Darcy scale behavior of two-phase flow in porous media is strongly controlled by porescale mechanisms. Pore invasion, snap-off, trapping, and mobilization are some of these pore-scale mechanisms that should be included in pore-network modeling. For brevity, the details of piston-like movement, entry pressure, trapping, and snap-off are not mentioned here and can be found in Joekar-Niasar and Hassanizadeh (2011a).

Conductivities of Pore Throats and Pore Bodies
Conductivities of pore throats are determined based on their size and fluid occupancy. One of the following three cases may occur. (a) A pore throat and the connecting pore bodies are filled by the wetting phase only. For this case, the following equation was obtained by (Azzam and Dullien 1977): where μ w is the viscosity of the wetting phase, l i j is the length of pore throat. (b) A pore throat is occupied by both phases. Following Ransohoff and Radke (1988), it is possible to write: In Eq. 13, is a resistance factor that depends on geometry, surface roughness, crevice roundness, and other specifications of the cross-section. (Detailed explanation about can be found in Zhou et al. 1997). As mentioned above, the pore throat capillary pressure p c i j is set as equal to the capillary pressure of the upstream pore body. (c) A pore throat is filled by the wetting phase, but the neighboring pore bodies are in the non-wetting phase. The wetting phase conductivity in the pore body i next to the pore throat i j is denoted by K w i,i j , but is considered only if both phases are present in a pore body. The conductivity of the non-wetting phase in a pore body is usually much larger than that of the pore throat. Therefore, only pore throat conductivity for the non-wetting phase is considered. However, relations given by Patzek (2001) are employed to calculate the wetting phase conductivity along the edges in a pore body with half-corner angle of β = 11 36 π.
, and the fitting parameters a 1 = −18.206, a 2 = 5.883, a 3 = −0.352. Where there is a disconnected blob of the non-wetting phase in a pore body, it is assumed that the blob will cover the opening of the pore throat that has the maximum outwards wetting phase gradient. Thus, the effective length for resistance should be corrected. Then, in Eq. 14, √ 6 2 R i − 2r i j will be replaced by the length of the blob. The length of a blob, l b , with a given volume V b constricted within the octahedron pore bodies can be geometrically approximated by solving the following equation for l b : Finally, the total wetting phase conductivity of a pore throat, bordered by two pore bodies containing the non-wetting phase, can be approximated as the 1 where K w i j is given by Eq. 13.

Local Capillary Pressure Curves
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. A corresponding capillary pressure and wetting phase volume fraction can be determined for a given fluid-fluid interface position within a pore body; therefore, a unique relationship between capillary pressure, p c i and local wetting phase volume fraction, s w i is obtained. Resulting local p c i -s w i curves for drainage and imbibition are discussed in Joekar-Niasar and Hassanizadeh (2011a).

Local Interfacial Area Curves
Drainage To calculate the interfacial area within the pore bodies, three different ranges of local s w i (and consequently local capillary pressure) are distinguished: (i) s w i is larger than the s dr i , which is the volume fraction obtained by an inscribed sphere within the pore body.
In this range, the non-wetting phase is assumed to fill up a sphere with radius r c i . The interfacial area is then equal to 4πr c i 2 . (ii) when s w i is smaller than s * i . s * i is the wetting phase volume fraction at which, all pore throats connected to the pore body i can be filled by both phases. To obtain this condition, a minimum capillary pressure is required that is denoted by p i c e and the corresponding interfacial area is denoted by A nw dr . Therefore, the interface along the edges in a pore body will have a shape as shown in Fig. 2b. This type of interface, which is sandwiched against the vertices is referred to as the "arc"or "corner" interface (Mason and Morrow 1987). The total length of the edges is a pore body is given by: Given L i , the total interfacial area reads Then, A nw i = A nw dr , when p c i = p i c e . (iii) When s w i is between s dr i and s * i . In this case, pore throats may be filled by the wetting and non-wetting phases. The interface geometry may be complicated and not necessarily along the flow channel, even bridging over several pores (cooperative filling). To resolve the exact geometrical configuration of an interface within a pore element, an interface simulator for given geometries is required (e.g., see the Surface Evolver software Brakke 1992). In this way, exact equilibrium pore-filling mechanisms (Lenormand and Zarcone 1983) at the pore level may be calculated. However, this approach is computationally expensive. Since the geometry of pore space is well defined and an approximate quantification of interfacial area is required at the Darcy scale, the idealization of the interface geometry may result in an estimation of specific interfacial area to coincide with the experiments (for example, see Joekar-Niasar et al. 2010b). Thus, we assume that the interfacial area in this range is obtained from a simple interpolation of the previous two ranges. It means that, if s where A nw dr is obtained from Eq. 18. The power of the 2/3 is assumed based on the contribution of area to the volume. Now, we can summarize the analysis in the following three ranges of s w i (related to the mean curvature).
Imbibition Similar to the drainage case, a relation for interfacial area versus s w i for a pore body, under imbibition can be written for three ranges of local fluid volume fraction. The only difference is that one must account for cooperative filling (the so-called I-pore-filling mechanisms in Lenormand and Zarcone 1983). As with the p c i − s w i relation (Joekar-Niasar and Hassanizadeh 2011a), this is done by considering the s w i at which only a single pore throat is filled by both phases. The corresponding wetting phase volume fraction is denoted by s imb i . Obviously, since snap-off capillary pressure is smaller than entry capillary pressure, the value of s w i at which snap-off happens will be larger than that of entry capillary pressure. Therefore, s imb i > s dr i . Thus, the ratio , defined above, can be larger than one. This results in larger local interfacial area under drainage in comparison to the drainage, which is in agreement with the pore-scale observations of Lenormand and Zarcone (1983). Finally, the A nw i − s w i relation under imbibition is as follows: Range I only has main terminal interfaces, while range II only has corner interfaces. However, both types of interfaces exist in range III and their ratio needs to be determined. This decomposition is done by considering pore throats' non-wetting phase occupancy. Thus, if all pore throats of a given pore body are filled by the non-wetting phase, the interfacial area is considered to be all corner interfacial area. However, if only one pore throat is filled by the non-wetting phase and the s w i is in range III, one-sixth of the total interfacial area is assumed to be the corner interfacial area and the rest would be main terminal interfacial area. Assuming a constant non-wetting phase pressure inside the blob, the capillary pressure at the tip of the blob can be estimated from the local wetting phase pressure gradient (see Eq. 23)

Mobilization of a Non-Wetting Blob
A discontinuous phase blob will be mobilized if the local capillary pressure, in the downstream side of the blob, can overcome the local entry capillary pressure. At the moment of an invasion, due to the continuity of the pressure field, an interface at upstream will have to recede. The local capillary pressure, which is related to the curvature of the interface is controlled by the local non-wetting and wetting phases pressure difference. In order to include the mobilization, the pressure gradient in both (continuous and discontinuous) phases must be calculated. The pore-scale two-pressure algorithm of DYPOSIT model provides the distribution of local capillary pressures at the interfaces of (the discontinuous) phases. This approach is much more efficient than the other pore-scale approaches for mobilization such as the search algorithm proposed by Ng and Payatakes (1980) or the definition of single pressure as proposed by Dias and Payatakes (1986). As shown schematically in Fig. 3, breakup, mobilization, stranding, and coalescence are all included in the model. However, coalescence is assumed to occur spontaneously. Because the discretization level is the pore body, there is only one pressure value for each phase in the pore body. But, for determining whether a disconnected blob in a single pore body (Fig. 4) will be mobilized, sub-scale information is required. The wetting phase pressure drop along a single non-wetting blob with length l b , is estimated as follows: where, p w i−i j is the wetting phase pressure drop along the blob. The conductivities are defined in Eqs. 13-14. The non-wetting phase pressure inside the blob is assumed to be constant, and assign the values of p w i and p n i from global calculations to the center of l b . Then, the capillary pressure at the tip of the blob, p c i tip , is estimated as follows: If the p c i tip is larger than the entry capillary pressure, the blob will be mobilized.

Boundary and Initial Conditions
Two types of simulations-quasi-static and dynamic-were performed in this study. Details of the simulation procedure under equilibrium and non-equilibrium conditions can be found in Joekar-Niasar and Hassanizadeh (2011a). The quasi-static simulations cover primary drainage, main imbibition, and several scanning drainage and imbibition curves. For quasi-static simulations, constant phase pressures at inlet and outlet were applied. In the non-equilibrium simulations, a constant flux at the inlet and a constant pressure at the outlet were applied. In all simulations, periodic boundary conditions were imposed along the network sides. For non-equilibrium simulations, varieties of initial conditions can be assumed: (a) an equilibrium initial conditions, where there is no movement of interfaces; (b) a non-equilibrium initial conditions, which non-equilibrium effects are still present in the spatial distribution of s w i and interfacial area. Table 2 shows specifications of dynamic simulations.

Averaging Procedure
Our simulations result in local-scale variables such as phase pressures, capillary pressures, phase volume fractions, and fluxes at each time step. These must be averaged over an averaging window to obtain macroscopic variables. The average saturation is simply defined as the ratio of the wetting phase volume to the total pore volume of the network: S n = 1 − S w in which n pb is the total number of pore bodies. Specific interfacial area is calculated using the following equation: Also, an average capillary pressure is calculated as follows: This capillary pressure is not equal to the phase pressure difference (P n −P w ) under non-equilibrium conditions. In the literature, the macroscopic phase pressure difference (P n − P w ) under non-equilibrium conditions is referred to as capillary pressure. However, these two quantities are separate entities and are equal only under equilibrium conditions. These issues are beyond the scope of this study and have been discussed in detail in Joekar-Niasar and Hassanizadeh (2011a) and reviewed in Joekar-Niasar and Hassanizadeh (2010).

On the Uniqueness of Equilibrium P c -S w -a nw Surfaces
This section analyzes the uniqueness of P c -S w -a nw surfaces under equilibrium conditions. This analysis corroborates previous experimental and computational studies (e.g., Reeves and Celia 1996;Held and Celia 2001;Cheng et al. 2004;Chen et al. 2007;Joekar-Niasar et al. 2008, 2010bPorter et al. 2009). It is presented here in order to establish equilibrium P c -S w -a nw surface needed for comparison to the corresponding non-equilibrium surface. The simulation results relating interfacial area to saturation are shown in Fig. 5. Similar to the capillary pressure-saturation curves, specific interfacial area-saturation curves are also non-unique. Under equilibrium conditions, the dynamic effects are absent; therefore, the relationship between capillary pressure, saturation, and interfacial area is controlled only by intrinsic properties of the porous medium such as pore topology and geometry. The influence of the geometry on the magnitude and trend of specific interfacial area is shown here. As previously mentioned, capillary interfaces can be classified into main terminal interfaces and arc menisci (or corner interfaces). The existence of arc menisci is strongly related to the angularity of the grains. For instance, in glass bead packings, where the solid surface is smooth, the arc menisci are almost absent and the main terminal interface has the major contribution to the specific interfacial area. In a very angular cross-section, however, the arc menisci comprise most of the specific interfacial area and their maximum value approaches the specific surface of the solid phase. The main terminal specific interfacial area increases with invasion as long as the newly generated interfacial area compensates for the loss of interfaces as a result of coalescence. To visualize this behavior, the variations of main terminal interfacial area with saturation for drainage and imbibition have been shown in Fig. 5b, d, respectively. The magnitude of the main terminal interfacial area is small compared with the corner interfacial area in this specific porous medium (because of the aspect ratio and the angularity of pore bodies). Furthermore, a peak of main terminal interfacial area at the saturation of approximately 0.2 clearly occurs, while the peak for total specific interfacial area occurs at a smaller saturation. This result is similar to the calculations of Raeesi and Piri (2009) for Berea sandstones, who found a peak saturation as low as 0.1-0.2. The variation of total specific interfacial area is almost monotonic under drainage, with fewer kinks compared with the main terminal interfaces. Finally, full set of data points are used to generate contour lines of equal interfacial area in the plain of capillary pressure and saturation. Figure 6a, c shows the mean equilibrium contour lines of total specific interfacial area and main terminal specific interfacial area, respectively. It should be noted that a linear interpolation has been used to generate the contour plots.
To show the uniqueness of the contour plots, the relative error contour plots have been calculated and presented in Fig. 6. The relative error in interfacial area was calculated as (X − < X >)/ < X >, in which X is the objective plot, and < X > is the average of drainage and imbibition contour plots. The relative error contour for drainage is the same as imbibition, but with the opposite sign. The deviation from the mean surface is minimum at Fig. 6 Mean equilibrium contour plots of specific interfacial area a nw versus capillary pressure P c and saturation S w for a all interfaces and c main terminal interfaces. The plots have been generated using a linear interpolation. The relative difference between the drainage contour plots (not shown) and the mean contour plots (a, c) are shown in (b, d), respectively the intermediate saturation and capillary pressure values. Comparison between the relative difference contour plots show that uniqueness for total specific interfacial area is much better than for the main terminal interfacial area. The average absolute relative errors are ∓8.5% and ∓17.5% for total specific interfacial area and main terminal interfacial area, respectively. In previous studies, a simple quadratic polynomial equation has been fitted to such a surface (Joekar-Niasar et al. 2009;Porter et al. 2009). However, the polynomial equations are not restricted within the physical range of variation. For instance, at saturation S w equal to 0 and 1, the specific interfacial area can be non-zero. To have a physical-based range of variation, the following positive-valued relation has been fitted to the data: Please note that the units of P c and a nw are kPa and 1/m, respectively. The fitting coefficients are α 1 = 6.462, α 2 = 3.057 × 10 −12 , ł pha 3 = 1.244, α 4 = −0.963. In this specific porous medium, the term S wα2 can be set to 1 as the value of a 2 is very small. Correlation between the fitted surface and the mean surface (interpolated from the data) and the mean relative error are 0.99 and 0.1, respectively. Note that in contrary to polynomial equations, this equation will result in zero interfacial value for S w = 0 and S w = 1. As can be seen, the uniqueness of the P c -S w -a nw surfaces under drainage and imbibition in this study is significantly better than the reported values in the literature (Reeves and Celia 1996;Held and Celia 2001;Joekar-Niasar et al. 2008, 2010bPorter et al. 2009). The reason is that in the previous modeling studies, the corner interfaces were almost absent or not calculated. In visualization studies, too, either the corner interfaces were absent (e.g., in glass beads) or they could not be visualized because of technical or resolution limitations (e.g., Cheng et al. 2004). In natural porous media, such as crushed sands or rocks, roughness is a major property of the solid phase, and the probability of having significant corner flow and arc menisci is high. In other words, in these natural porous media, the uniqueness of P c -S w -a nw surface may be statistically better than that in granular smooth porous media such as glass beads.
4.2 Non-Equilibrium P c -S w and a nw − S w Curves There are major differences between equilibrium and non-equilibrium situations. While in equilibrium simulations, fluid and capillary pressures are constant throughout the domain (in the absence of gravity), they are variable in non-equilibrium simulations. The non-uniform distribution of capillary pressure at a given time and fluid pressures depend on boundary conditions and fluid viscosities, which also leads to many different configurations of fluid phase and interfaces. Because of local invasion phenomena, there might be a distribution of disconnected non-wetting phases with their own capillary pressures especially under imbibition conditions. Note that all non-wetting phases (connected or not), capillary pressures, and their interfacial areas have been included in calculating of P c -S w -a nw surfaces. The disconnected phases can be mobilized, become fragmented, or be united to other blobs; all of these will lead to the change of capillary pressure and interfacial area. Figure 7 shows two snapshots of an imbibition simulation. The initial condition of the simulation is the end point of the main imbibition curve at S w = 0.49 obtained by the quasi-static simulator ( Fig. 7a. Thus, the non-wetting phase has lost its connection to the outlet completely and saturation increase will no longer be possible under quasi-static simulations. However, with a high capillary number of 10 −5 , the non-wetting phase can be mobilized to move out of the domain so that S w = 0.86 (Fig. 7b). Obviously, during this non-equilibrium process, many new interfaces will be generated or destroyed. The macroscopic values of capillary pressure, and specific interfacial area (for both total and main terminal) are plotted versus saturation in Fig. 8 for equilibrium and non-equilibrium conditions. For comparison between the nonequilibrium results (where the interfaces are unstable) and equilibrium conditions (where the interfaces are stable), equilibrium curves have been also plotted. As shown in Fig. 8, the non-equilibrium and equilibrium curves are not identical. Figure 8a shows that, at a given saturation, the transient capillary pressure under drainage condition can be even larger than the equilibrium capillary pressure. Under constant flux conditions at the inlet boundary, the internal capillary pressure at the pores will be controlled by the geometrical constrictions.
The capillary pressure has to build up locally so that interface can invade new pores and the flux boundary conditions can be satisfied. If the interface encounters another constriction, capillary pressure built-up will recur. Because of this issue, the difference between equilibrium and non-equilibrium capillary pressure-saturation curves has to be finite. Moreover, this difference is more pronounced under imbibition conditions than drainage. Under equilibrium conditions during imbibition, if the non-wetting phase is connected to the outlet boundary, the interface will be relaxed (small capillary pressure). Under imbibition, the viscous forces in the system do not allow all interfaces to be relaxed. The interfaces connected to the invading wetting fluid's front can be relaxed (or even disappear completely by the pore filling), while the interfaces far from the invading front will keep their capillary pressures. That is why, at low saturations under a large capillary number, the macroscopic imbibition capillary pressure-saturation curve follows the drainage curve. With the reduction of viscous forces in the system, the non-equilibrium capillary pressure-saturation curve will be close to the Fig. 7 Mobilization of the non-wetting phase under a forced imbibition process: a initial condition obtained at the end of the main quasi-static imbibition curve, S w = 0.49. b End of the non-equilibrium simulation. The major part of the trapped non-wetting phase has been mobilized due to a high capillary number of 10 −5 ; S w = 0.86 equilibrium curve. This trend can also be seen for total specific interfacial area versus saturation curve (shown in Fig. 8b). In this case, the variation of total specific interfacial area under imbibition is more sensitive to the non-equilibrium conditions than under drainage. This difference can be explained by the larger number of pores that can be imbibed by the wetting phase at a given saturation compared with those that can be drained under drainage. Greater possibility for invasions provides greater possibility for change of interfacial area. The sensitivity of the variation of main terminal and total specific interfacial area to non-equilibrium conditions is illustrated by Fig. 8c. From the order of magnitude of the interfacial area associated with main terminal interfaces and their variation with saturation, we can conclude that non-equilibrium conditions strongly influence the main terminal interfaces. In other words, the non-equilibrium effects on the corner interfacial area are smaller.
4.3 On the Uniqueness of Non-Equilibrium P c -S w -a nw Surfaces Figure 9a, b shows the contour plots of mean total specific interfacial area versus saturation and capillary pressure for capillary numbers 10 −7 and 10 −6 , respectively. The difference between the two non-equilibrium cases shown here (and more simulations were done which are not shown here) is small with a maximum relative difference of 0.09. Next, contour maps of interfacial area from equilibrium and non-equilibrium cases were compared. The contour plots of absolute values of relative difference between non-equilibrium and equilibrium contour plots are shown in Fig. 9c, d. With the increase of capillary number (increase of viscous forces), the deviation from the equilibrium surface increases. The analyses show that the magnitude of relative difference between the mean non-equilibrium and mean equilibrium P c -S w -a nw surfaces increases from 0.03 to 0.13 with the increase of capillary number from Ca = 10 −7 to Ca = 10 −6 . A closer look at the relative (a) (b) (c) (d) Fig. 9 Mean contour plots of specific interfacial area versus capillary pressure and saturation for capillary numbers of a 10 −7 and b 10 −6 . The contour plots of absolute values of relative difference between non-equilibrium and equilibrium contour plots have been shown in (c, d) difference contour plots in Fig. 9c, d shows that the major discrepancies occur at high saturations (out of the main loop), where the non-wetting phase at the end of imbibition is completely discontinuous. In other words, uniqueness of P c -S w -a nw surfaces decreases with changing from a continuous flow regime to a discontinuous flow regime. Note that in the original theory developed by Hassanizadeh and Gray (1993b), both phases were assumed to be continuous. Moreover, in their theory, solid-fluid interfacial areas may be as important as the fluid-fluid interfacial areas in describing the dynamics of the flow.

Final Remarks
This study employed a dynamic pore-network model, called DYPOSIT, that can simulate a wide range of viscosity ratios and capillary numbers. The DYPOSIT model includes major pore-scale mechanisms such as piston-like movement, snap-off, partial pore filling, change of capillary interfaces and their effect on local pore conductivity, mobilization of discontinuous blobs, capillary pinning, and relaxation of the interface due to local phase pressure variations. As visualization techniques in measuring non-equilibrium effects on capillary interfaces have limitations, the DYPOSIT model can serve as an investigative tool for the study of these effects as well as the interplay between specific interfacial area, macroscopic capillary pressure, and saturation. The macroscopic capillary pressure is totally different from macroscopic phase pressure difference, which is commonly (and falsely) referred to as capillary pressure. The macroscopic capillary pressure is an intrinsic property of the capillary interfaces, and basically depends on the configuration of the interfaces. Contrary to the macroscopic phase pressure difference, the capillary pressure does not scale with size of the averaging domain, as the specific interfacial area and local capillary pressure are not scale dependent. Furthermore, this equilibrium analysis illustrates that in angular media the uniqueness of P c -S w -a nw surface is more pronounced than that of smooth granular medium because of the contribution of main terminal and corner interfacial area to the total specific interfacial area. The total interfacial area is a monotonic function of saturation at saturation values larger than the residual saturation. This finding may be interpreted to indicate that uniqueness of the surface improves in angular/crushed media as the arc(corner) menisci can suppress the role of main terminal menisci in change of interfacial area versus saturation. This study's non-equilibrium simulations show that capillary pressure-saturation curves are not unique and depend on capillary number. Similarly, specific interfacial area-saturation curves also depend on capillary number. However, non-equilibrium drainage and imbibition P c -S w -a nw surfaces collapse on each other. But, the difference between mean equilibrium and non-equilibrium surfaces increases as the viscous forces increase, and mostly results from the change of pore-scale invasion mechanisms under imbibition conditions, such as suppression of snap-off and the fragmentation of the non-wetting phase with the increase of capillary number. These specific pore-scale mechanisms will result in discontinuity of the non-wetting phase that has not been considered in the original theory of Hassanizadeh and Gray (1993b). Nevertheless, the equilibrium P c -S w -a nw surface can provide an acceptable estimate of the non-equilibrium interfacial area within the range of practical uncertainties. For instance, for a capillary number of 10 −6 , the maximum discrepancy is about 30%, which is mostly observed at large saturations.