Transp Porous Med DOI 10.1007/s11242-016-0814-8
Pore-Scale Network Modeling of Three-Phase Flow Based on Thermodynamically Consistent Threshold Capillary Pressures. I. Cusp Formation and Collapse Arsalan Zolfaghari1 · Mohammad Piri1
Received: 8 February 2016 / Accepted: 20 December 2016 © Springer Science+Business Media Dordrecht 2017
Abstract We present a pore-scale network model of two- and three-phase flow in disordered porous media. The model reads three-dimensional pore networks representing the pore space in different porous materials. It simulates wide range of two- and three-phase pore-scale displacements in porous media with mixed-wet wettability. The networks are composed of pores and throats with circular and angular cross sections. The model allows the presence of multiple phases in each angular pore. It uses Helmholtz free energy balance and Mayer–Stowe–Princen (MSP) method to compute threshold capillary pressures for two- and three-phase displacements (fluid configuration changes) based on pore wettability, pore geometry, interfacial tension, and initial pore fluid occupancy. In particular, it generates thermodynamically consistent threshold capillary pressures for wetting and spreading fluid layers resulting from different displacement events. Threshold capillary pressure equations are presented for various possible fluid configuration changes. By solving the equations for the most favorable displacements, we show how threshold capillary pressures and final fluid configurations may vary with wettability, shape factor, and the maximum capillary pressure reached during preceding displacement processes. A new cusp pore fluid configuration is introduced to handle the connectivity of the intermediate wetting phase at low saturations and to improve model’s predictive capabilities. Based on energy balance and geometric equations, we show that, for instance, a gas-to-oil piston-like displacement in an angular pore can result in a pore fluid configuration with no oil, with oil layers, or with oil cusps. Oil layers can then collapse to form cusps. Cusps can shrink and disappear leaving no oil behind. Different displacement mechanisms for layer and cusp formation and collapse based on the MSP analysis are implemented in the model. We introduce four different layer collapse rules. A selected collapse rule may generate different corner configuration depending on fluid occupancies of the neighboring elements and capillary pressures. A new methodology based on the MSP method is introduced to handle newly created gas/water interfaces that eliminates inconsistencies in relation between capillary pressures and pore fluid occupancies. Minimization of Helmholtz free energy for each relevant displacement enables the model to accurately
B 1
Arsalan Zolfaghari
[email protected];
[email protected] Dept. 3295, 1000 E. University Avenue, Laramie, WY 82071, USA
123
A. Zolfaghari, M. Piri
determine the most favorable displacement, and hence, improve its predictive capabilities for relative permeabilities, capillary pressures, and residual saturations. The results indicate that absence of oil cusps and the previously used geometric criterion for the collapse of oil layers could yield lower residual oil saturations than the experimentally measured values in two- and three-phase systems. Keywords Wetting and spreading oil layers · Relative permeability · Spreading and non-spreading three-phase systems · Threshold capillary pressures · Two- and three-phase flow in porous media
Acronyms DI LC LF MSP MTM NAPL P3P PL
Displacement index Layer collapse Layer formation Mayer–Stowe–Princen Main terminal meniscus Non-aqueous phase liquid Point of three-phase Piston-like
1 Introduction A wide range of applications in different areas of science and engineering are concerned with the flow of three fluid phases in porous media. This includes the flow of oil, brine, and gas in geologic formations as well as flow of non-aqueous phase liquid (NAPL) in partially saturated soils. Partial differential equations that govern the flow at the continuum scale are often solved numerically to analyze the performance of these flow processes. The equations use as input macroscopic multiphase flow properties, e.g., relative permeabilities and capillary pressures, that have to be either measured or approximated using some other tools, such as empirical correlations or pore-scale models. It is particularly challenging to obtain these parameters experimentally for three-phase flow processes. They can be time consuming, expensive, and inaccurate (Oak et al. 1990; Dria et al. 1993; Honarpour et al. 1986). Number of parameters affecting the final results are numerous. For example, during three-phase flow processes (e.g., gas injection, wateralternating-gas drive, carbonated water injection, and flow of NAPL in aquifers) oil can flow at very low saturations establishing higher oil recoveries (Maloney and Brinkmeyer 1992; Maloney et al. 1997; Oak et al. 1990; Oak 1990; DiCarlo et al. 2000; Saraf et al. 1982; Alizadeh et al. 2014). At the pore scale, the flow of oil/NAPL at low saturations is controlled by spreading and wetting layers sandwiched between brine and gas in the crevices. This means that relative permeability of oil at low saturations depends on stability of the oil layers. An experimental observation of spreading oil layers and their ability to maintain oil connectivity at the pore scale is published by Alizadeh et al. (2014) where threephase intrapore fluid configurations are studied using high resolution micro-CT images. The stability, and therefore connectivity, of these layers are dependent upon parameters such as saturation history, spreading coefficient, wettability, capillary pressure, and pore geometry. Furthermore, we do not a priori know all the parameters/conditions, e.g., saturation history,
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
that the flow properties depend on, which in turn makes it impractical to measure them in the laboratory for all possible scenarios and displacement paths. A recent review by Alizadeh and Piri (2014) provides a detailed account of experimental studies of three-phase relative permeability and their controlling factors. As a consequence, empirical correlations, such as Stone (1970, 1973) and Baker (1988), with little or no physical basis have been used to predict three-phase relative permeabilities using two-phase properties as input. Therefore, properties predicted by these correlations may introduce significant uncertainties into the predicted flow behavior (Saraf et al. 1982; Grader and O’Meara 1988; Oliveira and Demond 2003; Delshad and Pope 1989). These models usually deviate further from the experimental data as the oil saturation decreases. One possible explanation for this behavior is that at low saturations of the intermediate wetting phase, three-phase flow physics is significantly different than that of two-phase flow. The uncertainty in prediction of three-phase relative permeability, particularly at low oil saturations, can drastically affect estimates of final oil recovery. This in turn can translate into loss of millions of dollars in field-scale applications. A physically based alternative is to use pore-scale modeling techniques to obtain more accurate estimates of these properties. They are developed and used to model two- and threephase flow processes and associated properties for a wide range of displacement scenarios and wettabilities. They are divided into two main categories: (1) direct numerical models and (2) pore network models. Direct numerical models, e.g., Lattice-Boltzmann (Gunstensen et al. 1991; Grunau et al. 1993; Ferreol and Rothman 1995; Chen and Doolen 1998; Kats and Egberts 1999), smoothed particle hydrodynamic (Monaghan 1988; Morris et al. 1997; Zhu et al. 1999; Tartakovsky and Meakin 2005; Liu and Liu 2003; Tartakovsky and Meakin 2006; Hu and Adams 2006), and moving particle semi-implicit (Ovaysi and Piri 2010, 2011) are computationally expensive and often exhibit poor numerical stability. Furthermore, these models have not been able to quantitatively predict the flow of multiple fluid phases in porous systems with varying wettabilities. The pore network models (Fatt 1956a, b, c; Blunt and Scher 1995; McDougall and Sorbie 1995; Blunt 1997, 1998; Øren et al. 1998; Xu et al. 1999; Dixit et al. 1999; Tsakiroglou and Payatakes 2000; Blunt 2001; Patzek 2001; Blunt et al. 2002; Suicmez et al. 2007, 2008), however, are not only computationally much less expensive but also have been successfully used to predict two- and three-phase flow properties (Øren et al. 1998; Patzek 2001; Valvatne et al. 2005; Valvatne and Blunt 2004; Piri and Blunt 2005a, b). This is mainly because they use as input pore networks that faithfully represent the pore space topology and also incorporate relevant pore-scale displacement mechanisms (Celia et al. 1995; Blunt 2001; Blunt et al. 2002, 2013). Three-phase pore-scale network models have been used over the last three decades to model displacement processes relevant to various three-phase flow applications. The models have become more sophisticated over the years by incorporating random pore networks, pores with angular cross sections, contact angle hysteresis, more realistic displacement mechanisms, spreading and wetting layers, and thermodynamically consistent threshold capillary pressures. Incorporation of each of these elements has consequently led to improvements in predictive capabilities of the models. Piri and Blunt (2005a) provide a detailed review of pore-scale network models of three-phase flow published prior to year 2005 (see Refs. Heiba et al. 1984; Soll and Celia 1993; Øren et al. 1994; Pereira et al. 1996; Paterson et al. 1997; Fenwick and Blunt 1998a, b; Mani and Mohanty 1997, 1998; Laroche et al. 1999; Laroche 1999; Hui and Blunt 2000; Lerdahl et al. 2000; Larsen et al. 2000; Van Dijke et al. 2001a, b, 2004, 2000a; Van Dijke and Sorbie 2002b). Except for the work of Øren et al. (1998), the models developed to that date had used regular pore networks.
123
A. Zolfaghari, M. Piri
In order to review the models developed since year 2005, we first discuss a brief review of threshold capillary pressures. Authors have used variety of techniques to estimate threshold capillary pressures of different displacements. They were often empirically based or derived from geometric criteria, see, for instance, snap-off and oil layer collapse (LC) threshold capillary pressure equations in Hui and Blunt (2000). For piston-like (PL) displacements, however, thermodynamically consistent Mayer–Stowe–Princen (MSP) method (Mayer and Stowe 1965; Princen 1969a, b, 1970) was used to find the relevant threshold capillary pressures. Lago and Araujo (2001, 2003) applied the MSP method in order to calculate threshold capillary pressures of the two-phase PL displacements in capillary tubes with polygonal cross sections and those with the curved sides. The authors studied only two-phase drainage process that led to quadratic threshold capillary pressure equations. Van Dijke and Sorbie (2003) and Van Dijke et al. (2004) extended the MSP method to three-phase (e.g., oil, water, and gas) systems for capillaries with uniform wettability and regular cross sections. They showed that capillary entry pressure for two-phase PL displacements is dependent on the third phase pressure, if it exists in the pore. Only phases with stronger wetting, compared to that of the bulk fluid, were allowed to form wetting layers. The authors ignored wettability alteration and consequent contact angle hysteresis. Extension of the MSP theory into pores with mixed-wet (non-uniform wettability) was presented by Piri and Blunt (2004). They used MSP method to study two- and three-phase PL displacements in capillary tubes with irregular triangular cross sections. They adopted Kovscek wettability alteration scenario (Kovscek et al. 1993) to form mixed-wet pores. Van Dijke et al. (2007) later extended this work to other three-phase PL displacements in mixedwet pores. They, however, used pores with regular cross sections. Van Dijke and Sorbie (2006) also studied the existence of wetting layers in capillary tubes with star cross section and mixed wettability. The authors used MSP method to analyze the formation and collapse of wetting oil layers in two-phase systems of varying wettability. Here, we present a brief review of the network models developed since year 2005 and the advances they have achieved. Piri and Blunt (2005a, b) developed a quasi-static random network model incorporating numerous fluid configurations for two- and three-phase flow. Wetting and spreading layers, wettability alteration, contact angle hysteresis, and double displacements were all incorporated in the model. Different displacement mechanisms such as piston-like, snap-off, pore-body filling, layer formation, and layer collapse were considered. A new methodology for saturation path tracking was used to compare predicted and measured relative permeabilities for the exact three-phase saturation paths followed by the experiments. They used a Berea sandstone network constructed by Øren et al. (1998) to match two-phase oil/water and gas/oil relative permeabilities measured in water-wet Berea sandstone core samples (Oak 1990). Using Bartell and Osterhof (1927) constraint, they calculated the range of gas/water contact angles. Using the calculated values, there were able to predict two-phase gas/water relative permeabilities. Keeping all values of contact angles fixed, they compared their predictions of three-phase relative permeabilities of various saturation histories with their experimental counterparts reported by Oak (1990). They studied the impact of initial oil saturation, wettability, and saturation history on the relative permeability curves. Oil relative permeabilities were overpredicted at low oil saturations, and as a consequence, residual oil saturation at the end of gas injection were underestimated. These were believed to be due to the use of geometric criteria for formation and collapse of spreading oil layers. Suicmez et al. (2006, 2007) incorporated imbibition–imbibition (water–oil–gas) and imbibition–drainage (water–gas–oil) double displacements into the model developed by Piri and Blunt (2005a) to study water-alternating-gas flooding. For water-wet systems, water relative permeability was found to be smaller when gas was present. They attributed this
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
observation to the fact that gas pushes trapped oil into smaller pores and throats, consequently increasing oil/water capillary pressure. For the gas injection into waterflood residual oil, they found a larger hysteresis in the water-wet media due to double drainage mechanism. They showed that relative permeabilities were independent of saturation path when plotted as a function of flowing (non-trapped) saturation. By anchoring a relatively simple regular pore network model to the measured two-phase relative permeability and capillary pressure data, Svirsky et al. (2007) predicted three-phase relative permeabilities and capillary pressures for water-wet Berea sandstone. Spreading oil layers were not included in the model and hence, three-phase oil relative permeability was underpredicted when compared to the experimental values reported by Oak (1990). The model could not capture saturation paths observed in the experiments. Al-Dhahli et al. (2013) presented a three-phase pore-scale network model with variable wettability and thermodynamic criteria for oil layer formation and collapse. All pores and throats had regular shapes. Consequently, layers formed or collapsed simultaneously in all corners of a given pore or throat. Authors compared the computed relative permeabilities against the experimental values reported by Oak (1990) for Berea sandstone. Compared to the results reported by Piri and Blunt (2005a, b), the model was able to reduce the discrepancy with the experimental data at low oil saturations, but the values were still overpredicted. The model did not follow the experimental saturation path and was not able to capture the experimental residual oil saturation at the end of gas injection. Three-phase gas relative permeabilities were not reported for the same saturation histories. In this work, we improve the predictive capabilities of the model developed by Piri and Blunt (2005a, b, 2004). We incorporate thermodynamically consistent threshold capillary pressures for all possible two- and three-phase displacement scenarios relevant to gas and water injections using MSP analysis. Similar to the work of Piri and Blunt (2004), Kovscek wettability alteration scenario is used to establish pores with non-uniform wettability (mixedwet). We perform stability analysis of wetting and spreading layers in pores with irregular angular cross sections. We focus on PL displacements, with and without layer formation (LF), and LC events in systems with varying wettabilities. A new method to handle gas/water interfaces which are formed due to collapse of oil layers in three-phase systems is presented. We also implement, for the first time, oil cusp analysis in our three-phase network model. We show that gas injection into oil can directly form oil cusps as well as spreading oil layers. Furthermore, these oil layers collapse to cusps. Both of these in turn impact the behavior of three-phase oil relative permeability at low oil saturations significantly. We rigorously validate our predictions against three groups of experimentally measured two- and threephase relative permeabilities in Zolfaghari and Piri (2016). In this paper, we first briefly describe our three-phase network model. This is followed by discussions of pore fluid configurations and displacement mechanisms, layer formation and collapse, cusp formation and shrinkage, and threshold capillary pressures.
2 Three-Phase Network Model Pore-scale network models combine various displacement mechanisms to mimic multiphase flow in naturally occurring porous media (Celia et al. 1995; Blunt 2001; Blunt et al. 2002, 2013). Different displacements such as piston-like, snap-off, pore-body filling, and layer formation and collapse are performed in a network of pores connected by throats to simulate fluid flow. Displacements are sequenced according to the state of wettability, saturation history,
123
A. Zolfaghari, M. Piri
pore geometry, and displacing and defending fluids (Lenormand et al. 1983; Lenormand and Zarcone 1984; Øren et al. 1998; Piri and Blunt 2005a). These displacements are modeled in random pore networks representing three-dimensional pore space in a given porous medium. The pores may have different cross-sectional shapes, e.g., circular, triangular, square, starshape, and lenticular (Van Dijke and Sorbie 2002b; Fenwick and Blunt 1998a, b; Larsen et al. 2000; Goode and Ramakrishnan 1993; Man and Jing 1999; Pereira et al. 1996), depending on, for instance, the geometry and size of the corresponding pores in the actual medium. Most of the elements in network models have angular cross sections, which allow them to hold more than one fluid phase at a time (Ehrlich and Davies 1989; Mason and Morrow 1991; Øren et al. 1998; Patzek 2001). A displacement involves a fluid configuration change in a pore during which a defending phase is replaced by an invading phase. The order by which these displacements take place depends, in capillary-dominated flow, on their threshold capillary pressures and other constraints such as the displacement site being accessible to the invading phase. This means that the threshold capillary pressure directly affects the evolution of fluid occupancy during a flow process and consequently the predicted flow properties. Three-phase pore-scale network model presented here is the result of new developments and substantial improvements implemented on the platform developed by Piri and Blunt (2005a, b, 2004). The model reads as input random pore networks and handles mixed-wet pores and contact angle hysteresis. New features and capabilities incorporated in the model include: thermodynamically consistent threshold capillary pressure calculations for various two- and three-phase displacements, including spreading and wetting layer formation and collapse, in irregular triangular elements, cusp formation and shrinkage analysis, different LC rules, and new continuous–continuous double displacement mechanisms. In this work, the pores are considered to be full of water and strongly water-wet initially r ) for oil primary drainage during before oil injection. We use a receding contact angle (θow which oil displaces water in the center of the pores leaving water as wetting layers in the corners. The displacement mechanism taking place during this process is PL displacement for which we calculate threshold capillary pressure using the MSP method. We then allow the wettability of the solid surfaces contacted by oil to alter. This is consistent with the wettability alteration scenario presented by Kovscek et al. (1993). The amount of solid surface exposed to oil depends on the maximum oil/water capillary pressure reached during oil drainage process and also the receding contact angle. For each pore or throat, the altered wettability boundary is defined as the location where two surfaces with different wettabilities meet. This boundary is stored for each element that has been exposed to oil at least once during the course of simulation. If multiple oil injection scenarios occur in a pore or throat, the scenario, in which oil/water capillary pressure reaches the maximum value for that element, is used to calculate the altered wettability boundary. After oil drainage, water or gas may be injected into the network depending on the selected injection sequences as the model allows for a variety of relevant displacement mechanisms (see Sect. 3 for more details). The wetting oil layers may form due to PL and snap-off displacements in weakly and strongly oil-wet systems. The formation of spreading oil layers, however, may take place in systems with varying wettabilities. They are strongly dependent on the spreading coefficient. In this work, we perform a detailed analysis of PL and LC events using the MSP method for both wetting and spreading oil layers. The threshold capillary pressures of those events are affected by the receding and advancing oil/water contact angles, receding gas/oil contact angle, geometry of the pore, and the maximum oil/water capillary pressure reached during primary drainage. For gas injection into oil, the relevant contact angle is the receding gas/oil value on the alt ). Oil/water capillary pressure and consequently oil/water altered wettability surface (θrgo
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
interfaces may stay unchanged throughout the network during a gas injection process. In other words, we do not allow any changes in capillary pressure between non-invading phases under any displacements as long as both non-invading phases are continuous. Therefore, for calculation of threshold capillary pressures, oil/water contact angle is kept fixed in any corner with a spreading oil layer. This is an important postulation in deriving all appropriate equations for cusp formation, cusp collapse, PL displacements with and without spreading oil layers as well as different collapse displacements presented in this paper.
3 Pore Fluid Configurations and Displacement Mechanisms 3.1 An Overview of Key Features A change in pore fluid occupancy constitutes a displacement in network modeling. Every pore has a fluid configuration at any point during a simulation by which the location of the fluid(s) residing in the element are defined. Pores may experience different fluid occupancies (or configurations) throughout a flow process depending on pore geometry, wettability, saturation history, and fluid/fluid interfacial tensions. Figure 1 presents the generic fluid configurations considered in this work for pores with irregular triangular cross sections. We have introduced a new pore fluid configuration (see configuration E in Fig. 1) that holds oil cusps at the corners. This configuration is different than those studied in literature (Van Dijke and Sorbie 2006; Kats et al. 2001) where two bulks of fluids are separated by a cusp interface, and capillary entry pressure is independent of cusp phase pressure and its geometry. In this work, we allow cusp interfaces to form at corners of angular pores by computing three geometrically related AM interfaces. That impacts capillary entry pressures and oil connectivity. As we
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 1 Two- and three-phase fluid configurations used in this study for capillary tubes with irregular triangular cross sections. Three of the most important displacement sequences implemented in the model include: A → F, A → D → E → F, and A → E → F
123
A. Zolfaghari, M. Piri
will explain later in this paper, the new configuration plays an important role in simulating some of the three-phase flow processes discussed. The number of fluid configurations is directly related to the number of corners with different corner half angles. For example, since irregular triangles have three different corner half angles, the number of fluid configurations, and hence the number of possible displacement scenarios, is much greater than those of equilateral triangles. In the interest of brevity, all of the possible fluid configurations are not shown here. We handle a numerous pore occupancy combinations with different corners having different configurations. For instance, a gas invasion into an irregular triangular pore that contains oil and water can form configurations with (a) oil layers in one, two, or all corners (configuration changes A → D), (b) corners with oil cusps (configuration changes A → E), and (c) corners with no oil (configuration changes A → F). We investigate all of the possibilities and identify the most favorable configuration change during each step of a flow simulation. We use thermodynamic analysis to obtain threshold capillary pressures for various valid displacements, which are then used to decide the next displacement during a flow simulation. We use the concept of corner-based configuration change, even if a displacement mechanism involves configuration change of multiple corners (e.g., displacements with multi-layer formation/collapse). Any part of a pore cross section that holds a fluid surrounded by solid walls and fluid/fluid interfaces is called phase location. Each phase location has a continuity flag which indicates whether the phase location is trapped or continuous. For each displacement, continuity flags are assigned to the newly created phase locations, and are updated throughout simulation of a flow process. There might be different displacement scenarios for a certain configuration change in an element. Each scenario is added to an appropriate list as a separate displacement. In the next few paragraphs, we briefly explain some of the new features of the model. The model allows the oil cusps (configuration E) to form due to gas injection into oil. Later in this document we will show that these cusps can exist under non-spreading conditions. This makes the displacement scenarios even more complicated for three-phase systems. These cusps are then allowed to shrink with the changes in the relevant capillary pressures. An additional complexity emerges when one handles the collapsing of the spreading layers, since they can collapse into cusp or directly into no-cusp configurations in each corner. In other words, the spreading oil layer(s) can collapse into gas/water interface through: (1) individual layer collapse into cusp configuration followed by a cusp shrinkage, (2) singular layer collapse, and (3) simultaneous multi-layer collapse (2LC or 3LC). For the case of spreading oil layers or cusp configurations, either reduction of oil/water capillary pressure (due to water injection) or increase in gas/oil capillary pressure (due to gas injection) can cause a collapse displacement. Oil layers should be connected to the outlet for the LC displacements to take place. This is because the displaced volume needs to be produced at the outlet. If oil is trapped, the double drainage mechanism could occur which involves two different types of displacements in a sequence (i.e., gas → oil followed by oil → water). More details are given by Piri and Blunt (2005a, b). We allow layers to collapse based on energy balance or geometric criteria depending on the type of invading phase locations present and collapse rule of the layers (see Sect. 3.4 for details on various LC rules implemented in the model). For the MSP layer collapse, an eligible Main Terminal Meniscus (MTM) in the neighboring element is required. A geometric layer collapse happens by increasing the pressure of fluids on either side of the layer to a threshold value, when an MTM in the neighborhood of the element is not present. Stability of spreading oil layers and cusp configurations control the connectivity of oil clusters during any gas injection process, and hence, profoundly affect the residual oil saturation and oil relative permeability at low oil saturations.
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
Different layer collapse events have different threshold capillary pressures calculated from the energy balance equations. These equations are derived for each layer collapse displacement. The layer collapse displacement associated with the lowest invading phase pressure occurs as the most favorable collapse event. Therefore, in theory, the collapse of spreading or wetting oil layers in irregular capillary tubes may happen individually or simultaneously. In regular capillary tubes, however, the collapse of oil layers will always be simultaneous. The model adds simultaneous multi-layer collapse events to the appropriate lists of displacements, if the following two conditions are met: (1) The oil layer is spreading, and (2) MSP layer collapse is allowed based on the selected collapse rule and configuration of neighboring elements (see Sect. 3.4 for the descriptions of various LC rules). Our results show that wetting oil layers collapse individually without any exception, see Sect. 4.1 for more details. Two different phase locations are involved in any displacement, i.e., invading and defending. The appropriate set of contact angles is determined based on the invading and defending phases as well as the surface (altered or original) on which the new interface is being created. However, the situation becomes unclear when a contact angle with the third phase (neither displacing nor displaced phase) is desired. For instance, consider collapse of a spreading oil layer in a corner of a pore by an eligible gas MTM located in a neighboring throat. Invading and defending phase locations are gas center in the neighboring throat and oil layer in the target pore, respectively. Since gas displaces oil on the altered surface, the appropriate gas/oil alt ). However, one does not, a contact angle is receding on the altered wettability surface (θrgo priori, know the appropriate gas/water contact angle for this scenario since water is neither a displacing nor a displaced phase for the collapse displacement. In the final configuration, the gas/water interface can be on the original or altered surface, or even hinging at the altered wettability boundary. Therefore, gas/water contact angle in the final configuration can have any value between receding on the original and advancing on the altered wettability surfaces org alt ). In this work, the gas/water contact angle in each corner is treated (i.e., θrgw ≤ θgw ≤ θagw as an unknown parameter in the system of nonlinear equations formed to obtain threshold capillary pressures. Additional water corner area constraints are added to this system of equations with the newly created gas/water interface. For instance, water area in each individual corner is not allowed to increase due to a gas injection into a pore, but it can remain constant or shrink. Threshold capillary pressure is calculated based on the minimum invading phase pressure, e.g., gas pressure for a gas invasion process, found from the system of nonlinear equations written for a given displacement. In the case of gas injection, the model investigates the consistency of calculated gas/water contact angle(s) with the appropriate surface that the interface is located on. One should note that, if any gas/water interface forms on the original surface during any iterations, code uses both surface parameters in the MSP method to form appropriate system of nonlinear equations accordingly. This situation can occur for any displacement creating a new interface involving the third phase. This is different from the work of Piri and Blunt (2005a) where a newly created interface was assigned the furthest contact angle of the previous interface just before the displacement. In the next two subsections, we provide details on the formation and shrinkage of oil cusps, as two of the critical novelties incorporated in the model to significantly enhance its predictive capabilities.
3.2 Cusp Formation In three-phase systems oil layers could collapse into a no-layer phase configuration or a cusp configuration. Collapse of a layer into a cusp happens for systems with negative spreading
123
A. Zolfaghari, M. Piri
coefficient if the cusp collapse displacement is favored over other collapse displacements. Oil cusp may form due to PL displacement by gas or water. Cusp connectivity is similar to the layer connectivity. Cusp phase is assumed to be connected to cusps, layers, and center phase locations of the neighboring elements if they hold the same fluid phase. Cusp configurations can also form, in theory, through geometric collapses during gas or water injection processes. For non-spreading systems, a droplet of oil placed on a water surface in the presence of gas forms an oil lens with a point of three-phase (P3P) contact. According to Neumann’s triangle (Rowlinson and Widom 1982), triple point’s force balance requires Ces < 0 at equilibrium. The net force at the P3P contact vanishes at equilibrium. Writing the relevant force balances at the P3P contact leads to (Rowlinson and Widom 1982): 2 − σ2 − σ2 σgw ow go ψo = ar ccos 2σow σgo σow + σgo − σgw σgw + σgo + σow = ar ccos 1 − (1) 2σow σgo where ψo is the phase angle for oil measured at the P3P contact (see Fig. 2). Similar equations could be written for other phases, i.e., gas and water. We use these phase angles in the cusp formation analysis. Three important facts are inferred from the Neumann’s triangle for any three-phase system at thermodynamic equilibrium: (1) Summation of any two interfacial tensions must be greater than the third one (triangle inequalities). This means the spreading coefficient cannot be positive for any system at equilibrium; (2) for spreading systems, the Neumann’s triangle becomes a slit-like triangle when one of the above-mentioned inequalities becomes an equality. This occurs when the summation of the two smaller interfacial tensions, i.e., oil/water and gas/oil, becomes equal to the largest tension, i.e., gas/water (ψo → 0). For a given system, there is a range of thermodynamic parameters, such as temperature, pressure, composition, and chemical potential, over which system behaves completely spreading (Rowlinson and Widom 1982); and (3) phase angles are dependent on the relative magnitude of equilibrium interfacial tensions rather than the difference in their magnitude. Mathematically, this means that the same phase angles could exist for different spreading coefficient values. Figure 2 illustrates a cusp interface in half of a corner. The other cusp in the corner is a mirror image of this configuration with respect to the middle line. The model assumes that there is no changes out of plane. Point A is the corner of the pore with the x- and ycoordinates of zero (i.e., x A = 0 and y A = 0). The positive direction of X and Y axes are toward left and down, respectively. Geometrically, two points and a radius are required to define a circle on a given plane. If bgo is the length of apex to the intersection of gas/oil interface with the pore wall (see Fig. 2), its coordinates can be uniquely defined by: xbgo r go cos θrgo + α = (2) ybgo r go cos θrgo + α cot (α) where α, θrgo , and r go are the corner half angle, receding gas/oil contact angle on the altered wettability surface, and curvature of the gas/oil interface, respectively. Since cusp formation is a result of gas-to-oil displacement, the appropriate contact angle would be receding value on the altered wettability surface for the gas/oil interface. Geometrically, the coordinates of P3P contact can be formulated as: x P3P bcontc sin (α) − H cos (α) = (3) y P3P bcontc cos (α) + H sin (α)
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
Fig. 2 Various parameters in a cusp configuration illustrated for half of a corner
where bcontc is the distance from apex to the intersection between pore wall and the perpendicular line passing through P3P contact, and H is the cusp height as indicated in Fig. 2. Gas/oil interface should pass through P3P contact and bgo where it meets the pore wall with contact angle of θrgo . Substituting Eqs. (2) and (3) into the general form of a circle equation, and considering positive gas/oil interface curvature, the x-coordinate for the center of gas/oil interface, i.e., xcgo , is calculated. The y-coordinate is calculated from the general circle equation passing through bgo coordinates with the radius of r go . Having xcgo , ycgo , and r go , the location of gas/oil interface is determined at the corner. The general solution for a gas/oil interface described above is not as simple as a quadratic equation which may arise for finding a circle with a given radius passing through two fixed points. In our case, both coordinates are dependent on the radius and hence capillary pressure. Here we are interested in finding an interface with the radius of r go passing through two points, i.e., P3P contact and bgo , whose coordinates are a function of gas/oil capillary pressure. This gives a system of nonlinear equations that must be solved numerically. Cusp analysis is an iterative procedure in which intersection of two interfaces are often needed to calculate P3P contact. Below, we present appropriate equations for such calculations. Details of iterative steps are given later in this section. Writing two circles’ equations with center coordinates of (xc1 , yc1 ) and (xc2 , yc2 ), the general equation for x coordinates of the intersection points can be formulated as: √ (yc2 − yc1 ) c1 + c2 ⎜ c3 x =⎜ ⎝ (yc1 − yc2 ) √c1 + c2 c3 ⎛
⎞ ⎟ ⎟ ⎠
(4)
123
A. Zolfaghari, M. Piri
where three constants of c1 , c2 , and c3 are defined as: c1 = (r2 + r1 )2 − (xc2 − xc1 )2 − (yc2 − yc1 )2 × (xc2 − xc1 )2 + (yc2 − yc1 )2 − (r2 − r1 )2 2 2 c2 = (xc2 + xc1 ) (yc2 − yc1 )2 + (xc1 − xc2 ) r22 − r12 + xc1 − xc2 2 2 2 2 2 2 c3 = 2 xc2 − xc1 + yc2 − yc1
(5) (6) (7)
where subscript c stands for center, and numbers 1 and 2 represent different circles. Two circles have no intersection if c1 < 0. For the case of c1 > 0, two circles have two distinct intersections even if yc2 = yc1 [see Eq. (4)]. Two circles have a tangent point, if c1 = 0. The model uses a two-circle intersection subroutine to determine the coordinates of P3P contact for a given pair of interfaces, i.e., gas/oil and oil/water, at any iteration of cusp calculations. During iterations, the model needs to calculate the intersection of a gas/oil interface with the pore wall to update the value of bgo . Considering point A as the origin in Fig. 2, the general form of a line equation that passes through the origin and makes an angle of size α with the positive direction of Y-axis is: π x y= −α x (8) = tan tan (α) 2 Combining a general interface equation with Eq. (8) gives a quadratic equation, which could be solved for the x-coordinate(s) of the intersection point(s): 2 π π 2 − α + 1 − 2x tan − α yc + xc + xc2 + yc2 − r 2 = 0 x tan (9) 2 2 The y-coordinate(s) for the intersection point(s) are obtained using Eq. (8) or the appropriate interface equation. For an intersection of a pore line and a given interface, contact angle toward the phase in the apex, i.e., θ apex , is derived as: ⎛ 2 2 ⎞ − x + y − y x int,2 int,1 int,2 int,1 π ⎠ (10) θ apex = − arccos ⎝ 2 2 |r | where xint,2 and xint,1 are the x-coordinates of intersections between a given interface and pore line, and yint,2 and yint,1 are the corresponding y-coordinates. |r | is the absolute value of the interface curvature. Here we explain the sequence of the calculations. For given values of θow and bow , the model first fixes the position of an oil/water interface depending on the oil/water capillary pressure. bow is the length of apex to the intersection of oil/water interface with the pore wall. In the next step, for a given gas/oil capillary pressure and gas/oil contact angle, we determine a locale of coordinates for gas/oil interface center with which the interface has an appropriate intersection with oil/water interface. Using this locale, an iterative method is developed to search for a gas/oil center location, which forms a cusp configuration with the minimum phase angle error. The concept of phase angle error is explained in details below. Note that here iterations are done for a fixed set of capillary pressures, bow , and θow values. P3P contact is then calculated as both interfaces for oil/water and gas/oil are fixed. If the interfaces have two distinct legitimate intersections, the one with a lower phase error is selected as the P3P contact. Having P3P’s coordinates and gas/water interface curvature
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
(both oil/water and gas/oil capillary pressures are known at this stage), gas/water interface is uniquely determined in the pore space with the criterion of x cgw = 0. By writing a force balance at the P3P contact, different phase angles are calculated for a given three-phase system (see, for instance, Eq. (1) for the oil phase angle). The model compares calculated phase angles from the force balance equations and their counterparts for a given cusp geometry at each iteration. Phase angle error is calculated for every geometrically feasible cusp configuration using Eq. (11): cusp
E phase angle =
|ψw
cusp
− ψw | + |ψo 2
− ψo |
(11)
Phase angles calculated from the force balance equations (i.e., ψw and ψo ) have fixed values for a given set of interfacial tensions [see Eq. (1)]. Their counterparts, on the other cusp cusp hand, calculated from geometrically constructed cusp interfaces (i.e., ψw and ψo ) vary at each iteration. That causes variations in phase angle errors at different iterations. Upon convergence, we report phase angles for the calculated cusp geometry. Here, we present equations used in the phase angle calculations for a given cusp configuration. All of the parameters are illustrated in Fig. 2. For any cusp configuration with given interfaces equations, water phase angle can be calculated from: ψwcusp = π + βow − ϕ
(12)
where ϕ is the angle of gas/water interfacial tension vector with the horizontal line at P3P contact, i.e., y = y P3P (see Fig. 2), which is given by: xcgw − x P3P (13) ϕ = arctan − ycgw − y P3P and βow is the angle of oil/water interfacial tension vector at P3P contact with the negative direction of x axis, which is expressed by: xcow − x P3P (14) βow = arctan − ycow − y P3P Similarly, gas phase angle from cusp configuration is written as: cusp
ψg
= π − βgo + ϕ
(15)
where βgo is the angle of gas/oil interfacial tension vector at P3P contact with the negative direction of x axis (see Fig. 2), and is calculated from: xcgo − x P3P βgo = arctan − (16) ycgo − y P3P Equations (11)–(16) are used to determine whether a cusp configuration is geometrically compatible with the given three-phase interfacial tensions. For every geometrically feasible cusp, MSP analysis and energy balance equation are invoked. This is to minimize the free energy of the displacement. A system of nonlinear equations written for a given cusp formation displacement includes two equations for each corner, i.e., corner water area [see Eq. (49)] and energy balance equations. Corner water area equation is written based on the similar concept introduced for creation of the new gas/water interfaces. Energy balance equation, on the other hand, is written based on the MSP analysis for a cusp formation displacement. Similar to the optimization procedure used for the newly created gas/water interfaces explained earlier, the
123
A. Zolfaghari, M. Piri
change in water area is found through an optimization algorithm for each corner to determine the minimum gas invasion pressure for a cusp formation displacement. One should note that the contact angle of oil/water interface is one of the unknown parameters in the system of nonlinear equations written for a cusp formation displacement. On the other hand, the location of oil/water interface is determined through the optimization process mentioned earlier. If final location of the oil/water interface on the pore wall is not compatible with the calculated oil/water contact angle, the interface location is adjusted accordingly. That changes water area in the corner. Since the wetting phase is well connected throughout the medium, for a given pore, any changes of the wetting phase volume is admissible. This is similar to the creation of new gas/water interfaces where gas injection changes water corner areas. As an example, imagine the optimized solution of the system of nonlinear equations suggests that the cusp configuration should form entirely on the altered wettability surface. If the calculated oil/water contact angle is greater than the advancing oil/water contact angle on the altered wettability surface, and water in the corner is a continuous phase, the oil/water interface should therefore move away from the apex to reach advancing oil/water contact angle. Then, it searches for the optimized location of gas/oil interface to minimize the phase angle error. The curvature and contact angle of the gas/oil interface are constant during the search. All possible geometric considerations are checked in order to detect any interface conflict for a cusp configuration. They are explained in detail here. The model checks for a physically meaningful P3P contact and intersections of the pore wall with the oil/water and gas/oil interfaces. It also ensures the oil/solid contact length in a cusp configuration has a positive value, and the location of P3P contact on oil/water and gas/oil interfaces is in agreement with each interface’s curvature. To explain the latter, P3P contact, for instance, cannot be on the lower part of oil/water interface (i.e., lower hemisphere in the corresponding circle equation) with a positive curvature (Pcow > 0) when oil/water interface crosses the pore’s middle line at y = ycow . Similar situations are considered for a negative oil/water curvature and a gas/oil interface, which can only have positive curvatures. The model ensures that there are no intersections of the gas/oil and oil/water interfaces between contact points on the pore wall and P3P contact. It also checks that gas/water interface is geometrically feasible with respect to oil/water and gas/oil interfaces. For instance, the model searches for any possible conflict between gas/water and oil/water interfaces when P3P contact is on the lower part of the oil/water interface with a positive curvature. Furthermore, we ensure that the entire oil/water interface is within the pore space for a given oil/water capillary pressure and bow value. As an example, consider a cusp configuration with positive curvature of oil/water interface. If contact point (calculated from bow value) falls beneath the center location, i.e., bow cos α ≥ ycow , the model checks for any possible conflict between oil/water interface and pore wall at x = xcow + row and y = ycow . This could be extended for the negative oil/water curvatures with contact point being above the interface’s center. Note that such a test for the gas/oil interface is not necessary. This is due to the fact that the center of gas/oil interface at each iteration is determined from its calculated locale. All of the interfaces from that locale create contact angle of θrgo with the pore wall. That prevents any conflicts with the pore wall for all gas/oil interfaces created in various iterations. For each geometrically feasible cusp configuration, all fluid–fluid and fluid–solid areas as well as fluid volumes are needed in order to calculate threshold pressure of invading phase for the cusp formation displacement. Appendix A lists all equations for different geometrical properties of a cusp configuration. The conductance of an oil cusp is estimated using the expression proposed by Hui and Blunt (2000) for the conductance of spreading layers. The cross-sectional area of the cusp is used
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
as an input to compute its conductance. Similar to the work of Piri and Blunt (2005a), no-flow and free boundary conditions are used for the oil/water and gas/oil interfaces, respectively.
3.3 Cusp Shrinkage Once a cusp is formed, it can shrink to no-cusp configuration and create a complete gas/water interface stretching from one side of the pore wall to another. Shrinkage of a cusp configuration can occur due to the pressure increase of either of the fluid phases surrounding it, i.e., gas or water. Similar to a layer shrinkage, no MTM is required for the shrinkage of a cusp. Shrinkage by increase in gas pressure is handled by fixing the oil/water interface at its original position determined at the last stage of the cusp analysis. The corresponding error during shrinkage is related to the phase angle error introduced earlier. The new curvature for the gas/oil interface is calculated for a given gas/oil capillary pressure. Next, the locale of the gas/oil center is determined for a fixed oil/water interface. The locale is calculated to have a constant gas/oil contact angle with the pore wall, i.e., receding value on the altered wettability surface, to create a valid P3P contact (valid intersection with the fixed oil/water interface), and to fulfill the bgo ≥ bow inequality. Final position of the gas/oil interface is found from the locale by minimizing the phase angle error for the newly calculated gas/oil curvature. Collapse of a cusp occurs at a threshold gas/oil curvature at which the cusp phase area, i.e., oil area, reduces to zero. Later in this document, we compare two important features of cusps and layers, i.e., oil areas and collapse threshold capillary pressures (see Sect. 4.2 for more details). The numerical procedure for the cusp shrinkage due to water injection is similar to that of the gas injection shrinkage explained earlier. That is because the corresponding error is related to the phase angle error. For the shrinkage of a cusp due to water injection, the curvature for the gas/oil interface remains unchanged. New curvatures are calculated for the current values of oil/water and gas/water capillary pressures due to water injection. Depending on the contact angle of the oil/water interface at the time of formation, the interface may hinge at its original position with the newly calculated curvature. The hinging continues until the oil/water interface reaches the advancing oil/water contact angle and move toward the center of the pore. Collapse threshold capillary pressure corresponds to an oil/water curvature at which oil area reduces to zero for a fixed gas/oil interface curvature.
3.4 Layer Collapse Rules During different flow processes, fluid phases may form spreading or wetting layers. Based on the connectivity criteria presented by Piri and Blunt (2005a), the layers in two neighboring elements are considered connected if they hold the same fluid. We use the same connectivity concept for layer–center, layer–cusp, and cusp–cusp connections between two neighbors. As mentioned earlier, an MSP LC displacement is allowed only if there is at least one legitimate MTM in the neighboring elements. Two fundamental questions arise from the above-mentioned concept: (1) What is the definition of a legitimate MTM?, and (2) Which of the followings has a priority, MTM availability or layer connectivity? In order to answer these questions and cover all possible scenarios, we have introduced four different rules for layer collapse displacements. They determine whether a layer is allowed to collapse based on the geometric criteria, energy balance or both. Some rules are strongly dependent on the configuration of the neighboring elements. The collapse rules are: 1. Geometrical under this rule, all layers holding the phase of interest (oil or gas) are collapsed only based on geometric criteria.
123
A. Zolfaghari, M. Piri
(a)
(b)
Fig. 3 Illustration of the minimum requirements for a geometric or MSP layer collapse displacement under different LC rules: a layer collapse would be geometric under MSP-Geometrical-a, and b layer collapse would be MSP under MSP-Geometrical-b. A triangular pore is shown with four different connecting throats (indicated by arrows). Here blue, red, and green represent water, oil, and gas, respectively. The red bar next to the arrows indicates the existence of a spreading or wetting oil layer in the connecting throat.
2. MSP-Geometrical-a this rule favors geometric criteria over the energy balance counterparts as long as there are at least two neighbors with the displaced phase in their layer or center phase locations. If this requirement is not met, MSP LC occurs. Figure 3a illustrates the minimum requirement for this rule under which geometric layer collapse is used for the layer in the target element. 3. MSP-Geometrical-b this rule favors the MSP LC over the geometric displacements as long as at least one legitimate MTM is found in one of the neighboring elements. A legitimate MTM exists in a neighboring element if the element has displacing phase in either its layer or center without any phase location holding the displaced phase. Figure 3b illustrates the condition under which at least one legitimate MTM is available in the neighboring elements to perform the MSP LC for layers in the target pore. 4. MSP All layer collapse events occur based on the MSP criteria regardless of the availability of legitimate MTMs in the neighboring elements and the type of their fluid configurations.
3.5 Double Displacements and Saturation Path Tracking In order to predict three-phase relative permeabilities, it is essential to accurately follow the experimental saturation path for which the properties are predicted. The model incorporates both direct and indirect (i.e., double displacements) displacements during gas injection to improve its ability to track experimentally measured saturation histories. The direct displacement of, for instance, water by gas during gas injection is implemented in the model as one of the main categories of single displacements. However, this is not sufficient to track the experimentally measured saturation paths. This is partly because of often significantly different interfacial tensions between the invading phase and the defending phases. For instance, in gas injections, gas/water interfacial tension values are usually much higher than their gas/oil counterparts. This results in displacement of most of the water at the end of gas injection in quasi-static pore-scale network models since direct displacements of gas into oil are far
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
more favorable than gas to water displacements. Gas injection experiments performed under capillary-dominated displacement regimes with relatively high initial oil and water saturations show saturation paths corresponding to significant simultaneous displacement of water and oil even though gas/oil interfacial tensions are significantly lower than their gas/water counterparts (Oak et al. 1990; Oak 1990; Alizadeh and Piri 2014). This could be due to the double drainage mechanisms at the pore scale causing an increase in oil/water capillary pressure during gas injections (Alizadeh and Piri 2014). Local pressure drops caused by the transport of the defending phase (i.e., oil) volume can trigger double drainage displacements. This could particularly occur when displacement site for the first drainage (i.e., gas-to-oil) is far away from the outlet or the spreading layers restrict oil conductivity. This is because pressure drops in the oil phase are higher if volume of the displaced oil is too high with respect to the oil conductivity between the displacement site and the outlet. Wetting phase is always connected throughout the medium and can be produced more easily. The increase in oil pressure should be high enough to allow invasion into water-filled elements accessible by the invading phase in the second drainage (i.e., oil-to-water). Experimental observations of the double drainage mechanisms in micro-models are reported by Øren and Pinczewski (1995). In these pore-scale displacement mechanisms, gas injection causes gas/oil and oil/water interface movements in a sequence. Similar evidences are also reported during gas injections in a naturally occurring sandstone with continuous oil and water phases using micro-CT and miniature core-flooding technologies (Khishvand et al. 2016). Convectional continuous-continuous single displacements in quasi-static network models always displace most of the connected oil before displacing any water. This leads to very different saturation path than those observed experimentally (for the same initial conditions). In order to reproduce the experimental saturation path, the model must allow oil/water capillary pressure to increase during gas injection. This could be the result of double displacement mechanisms as explained above. The model presented by Piri and Blunt (2005a) allowed double displacements in which a continuous phase displaces a trapped phase that in turn displaces another continuous phase. In this work, we also allow double displacements in which all three phases are continuous (“continuous–continuous” double displacements). To do this, oil pressure must increase during gas injection so that gas → oil → water displacement can take place. At each stage of gas injection, we tune the oil phase pressure to closely reproduce the experimental saturation path. Higher oil pressure increases the possibility of oil invasion into water as the second part of the above-mentioned double-displacement mechanism. The increase in gas pressure compensates for the threshold capillary pressures associated with both displacements. The model allows quasi-static displacement sequences for both gas-to-oil and oil-to-water displacements. That means that all simulated saturations and their calculated relative permeabilities are obtained under capillary equilibrium conditions where no more displacements can occur with the current values of the capillary pressures. Oil pressure is increased to carry out a set of most favorable oil-to-water displacements. Gas pressure is consequently increased to maintain constant gas/oil capillary pressure. This prevents any oil-to-gas displacements to happen. We also make sure that all gas/oil interfaces remain unchanged. During this process, we only allow currently available oil-to-water displacements to take place. In other words, oil is not injected at the inlet of the network. Oil-to-water displacements continue until system reaches capillary equilibrium, where no more oil displacement could occur with the assigned oil/water capillary pressure. Model decides whether to continue with the oil injection (i.e., another increase in the oil/water capillary pressure) or to resume gas injection by comparing gas and water saturations to those of the experiment. This approach was inspired by the saturation path tracking technique first introduced by Piri and Blunt (2005a, b), which was used to
123
A. Zolfaghari, M. Piri
follow an exact saturation path of a three-phase experiment and predict relative permeability curves. The difference is that during each flow process, we do not add new displacements to the sorted lists of available displacements. Generally, there are no measured values of oil/water capillary pressure available for gas injection experiments to quantitatively guide its evolution during saturation path tracking algorithms. We therefore use the experimental saturation history to adjust the oil/water capillary pressure at each stage of gas injection process to closely reproduce the measured path.
4 Threshold Capillary Pressures Each displacement mechanism is characterized by a threshold capillary pressure, and one can use MSP method to calculate them. For a displacement in an angular pore under capillarydominated displacement regime and in the absence of gravitational forces, the pressure difference between invading and defending phases across an MTM is equal to that of the AMs’ bounding similar phases at the corners of the same capillary tube. Assuming a thermodynamically reversible PL displacement, minimization of Helmholtz free energy (F) of the fluid configuration change associated with a displacement gives the threshold capillary pressure. This means that appropriate threshold capillary pressure is found such that the derivative of the Helmholtz free energy becomes zero for that configuration change. The application of this method to calculate threshold capillary pressure of various displacements have been presented by numerous authors (Lago and Araujo 2001, 2003; Van Dijke and Sorbie 2003, 2006; Van Dijke et al. 2004, 2007; Piri and Blunt 2004). Therefore, for brevity, we do not include details of the method here and only present the final equations for some of the critically important two- and three-phase displacements (PL and LC) in capillary tubes with irregular triangular cross sections. For each displacement we include the associated configuration change as depicted in Fig. 1. 1. Two-phase displacements – PL displacement without formation of any wetting oil layers (configuration changes A → C): −π a 2 h h At − row L t cos(θow ) + row + θow,2, + θow,1, 1 1 2 ξ1 + ξ2 + ξ3 h =0 (17) + θow,3, + 1 sin(α1 ) sin(α2 ) sin(α3 ) where j is the jth interface in a corner numbered from apex downward and ξ1 -ξ3 are given by: h a h ξ1 = sin(α2 ) sin(α3 ) cos(θow,1, + α1 ) 2 cos(θow ) − cos(θow,1, ) (18) 1 1 h a h ξ2 = sin(α1 ) sin(α3 ) cos(θow,2, + α2 ) 2 cos(θow ) − cos(θow,2, ) (19) 1 1 h a h + α3 ) 2 cos(θow ) − cos(θow,3, ) (20) ξ3 = sin(α1 ) sin(α2 ) cos(θow,3, 1 1 h where θow,i, is the oil/water contact angle hinging at the point of altered wettability 1 surface in the i th corner. Three auxiliary equations are written for oil/water interfaces
at every corner (i.e., bow,i = row r,orig
dr , θ and 3 ). row ow
123
h +αi ) cos(θow,i, 1 sin(αi )
r,orig
dr cos(θow +αi ) , where i=1, 2, = row sin(αi )
, At , and L t are the minimum oil/water curvature reached during
Pore-Scale Network Modeling of Three-Phase Flow Based on…
drainage, receding oil/water contact angle on the original surface, total pore area, and total pore perimeter, respectively. – PL displacement with the formation of three wetting oil layers (configuration changes A → B): a 2 a At − row L t cos(θow ) + row 2π − 3θow ξ1 + ξ2 + ξ3 a + cos(θow ) =0 (21) sin(α1 ) sin(α2 ) sin(α3 ) where a − α1 ) ξ1 = sin(α2 ) sin(α3 ) cos(θow
ξ2 = ξ3 =
a sin(α1 ) sin(α3 ) cos(θow a sin(α1 ) sin(α2 ) cos(θow
(22)
− α2 )
(23)
− α3 )
(24)
– PL displacement with the single collapse of a wetting oil layer in corner i: a h a a −π + θow + θow,i, − cos(θow ) sin(θow ) 1 h h a h + cos(θow,i, ) sin(θow,i, ) − 2 cos(θow ) sin(θow,i, ) 1 1 1 2 a ) − cos(θ h cos(αi ) cos(θow ow,i,1 ) − =0 sin(αi )
(25)
To solve for the threshold capillary pressure associated with each displacement, first we solve for their corresponding curvature (row ), which is explicitly involved in Eqs. h (17) and (21). In Eq. (25), however, the only unknown parameter is θow,i, , which 1 relates to the oil/water curvature through water/solid contact length equation (bow = cos(θ h
+αi )
ow,i,1 ). Knowing the curvature, threshold capillary pressure is calculated row sin(αi ) for each displacement based on the Young–Laplace equation written for capillary tubes ow with straight pore walls, i.e., Pcow = σrow . 2. Three-phase displacements
– PL displacement without formation of any spreading oil layers (configuration changes A → F): σow C3 2 r gw C2 + r gw C1 + − σgw C3 = 0 (26) row where C1 –C3 are given by:
σow
Aw,1 + Aw,2 + Aw,3 row + 2 σow row θow,1 + θow,2 + θow,3 − π cos θow,2 + α2 cos θow,3 + α3 cos θow,1 + α1 r + + L t − 2 row + σgo cos θgo sin (α1 ) sin (α2 ) sin (α3 ) ⎤ ⎡ orig r,orig/alt orig r,orig/alt orig r,orig/alt L + L s,2 cos θgw,2 + L s,3 cos θgw,3 cos θgw,1 ⎥ ⎢ s,1 ⎢ ⎫⎥ ⎧ ⎥ ⎢ cos θ + α + α + α cos θ cos θ ow,1 1 ow,2 2 ow,3 3 ⎪ + σgw ⎢ ⎬⎥ ⎨ 2 row ⎪ + + ⎥ ⎢ r,alt sin (α1 ) sin (α2 ) sin (α3 ) ⎥ ⎢ cos θgw ⎦ ⎣ ⎪ ⎪ ⎭ ⎩ orig orig orig − L s,1 − L s,2 − L s,3
C1 =
(27)
123
A. Zolfaghari, M. Piri ⎤ cos θgw,1 + α1 r,orig/alt cos θgw,1 ⎢ π − θgw,1 + θgw,2 + θgw,3 − ⎥ sin (α1 ) ⎢ ⎥ C2 = 2 σgw ⎢ ⎥ ⎣ ⎦ cos θgw,3 + α3 cos θgw,2 + α2 r,orig/alt r,orig/alt cos θgw,2 cos θgw,3 − − sin (α2 ) sin (α3 ) 2 θow,1 + θow,2 + θow,3 − π C3 = At − Aw,1 − Aw,2 − Aw,3 − row ⎡
cos θow,1 + α1 cos θow,2 + α2 + cos θow,2 + cos θow,1 sin (α1 ) sin (α2 ) cos θow,3 + α3 + cos θow,3 sin (α3 )
(28)
(29)
where θgw,1 , θgw,2 , θgw,3 , and r gw are four unknowns in this equation. Equation (26) and three independent corner area equations (written for the water phase in each corner) form the system of nonlinear equations that is solved for the unknown parameters. Equation (26) is derived in a general form and therefore each corner can have an independent water reduction value ( Aw,i = f × Aw,i , where f is the fraction of original water area in the i th corner that satisfies 0 < f < 1). Furthermore, r each interface can be on any of the surfaces with different wettabilities (θgw,i on original or altered surface, where i=1, 2, and 3). We determine these parameters such that gas invasion pressure for the displacement is minimized. The minimization procedure is done in a generic way that allows each Aw,i to be treated as an independent unknown parameter. For a faster convergence, one could assume that for all involved corners, the percentage of water area reduction is similar. A similar concept is used for other displacements with at least one newly created gas/water interface. – PL displacement with the formation of three spreading oil layers (configuration change A → D): ⎧ ⎨ r 2 r At − r go L t cos θgo + r go 3 θgo −π ⎩ ⎞⎫ ⎛ r +α r +α r +α ⎬ cos θgo cos θgo cos θgo 1 2 3 r ⎝ ⎠ =0 + + + cos θgo ⎭ sin (α1 ) sin (α2 ) sin (α3 ) (30) – PL displacement with the single collapse of a spreading oil layer in corner i due to gas injection: 3 2 r gw C1 + r gw C2 + r gw C3 + σgw C4 = 0
where C1 –C4 are given by:
(31)
cos(θgw,i + αi ) π θgw,i + αi − + cos θgw,i (32) 2 sin(αi ) 2 σow σgw cos θgw,i + αi r,orig/alt cos θgw,i − cos θgw,i C2 = (33) row sin (αi ) ( ) A cos θow,i + αi π w,i r + 2 σgo cos θgo C3 = σow σow 2 − θow,i − αi − 2 2 row sin (αi ) cos θow,i + αi σow σgw orig r,orig/alt orig r,alt + cos θgw L s,i cos θgw,i − L s,i − 2 row row sin (αi )
C1 =
123
σow row
2
Pore-Scale Network Modeling of Three-Phase Flow Based on…
r +α cos θgo i π 2 r r θgo + σgo + αi − + cos θgo 2 sin (αi ) cos θgw,i + αi π r,orig/alt 2 − cos θgw,i 2 cos θgw,i + σgw − θgw,i − αi − 2 sin (αi ) cos θow,i + αi
Aw,i r − 2 σgo row cos θgo C4 = σow + 2 row θow,i + αi − row 2 sin (αi ) cos θ + α ow,i i orig r,orig/alt orig r,alt − L s,i + cos θgw + σgw L s,i cos θgw,i 2 row sin (αi )
π
(34)
(35) – PL displacement with a single cusp formation in corner i due to gas injection: cusp cusp r go C1 + C2 r go , θow,i , f bow (36) − σgo C3 r go , θow,i , f bow = 0 where C1 is a constant value, and C2 and C3 are functions of the unknown parameters for a cusp analysis. They are given by:
cos θow,i + αi ⎤ π ⎥ ⎢ σow θow,i + αi − 2 − cos θow,i sin(αi ) ⎥ ⎢ ⎥ ⎢ cos θow,i + αi ⎥ ⎢ r C1 = row × ⎢ − 2 σgo cos θgo (37) ⎥ ⎥ ⎢ sin(α ) i ⎥ ⎢ ⎦ ⎣ π r + 2 σgw cos θgw − θow,i − αi (1 − f bow ) 2 ⎡ cusp ⎤ cusp Aw r go , θow,i , f bow cusp cusp cusp C2 r go , θow,i , f bow = σow ⎣ + L ow r go , θow,i , f bow ⎦ row ⎤ ⎡ π cusp cusp r , θ , f − α − θ r − 2 r L go bow go i go go ow,i ⎥ ⎢ 2 ⎥ ⎢ ⎥ ⎢ ⎫ ⎥ ⎧ bcusp r , θ cusp , f ⎢ go ow,i bow go ⎥ ⎢ ⎪ ⎪ ⎪ ⎥ ⎢ ⎪ ⎪ ⎪ ⎪ ⎪ + σgo ⎢ ⎬ ⎥ ⎨ − bcusp r , θ cusp , f ⎥ ⎢ go bow ow r ow,i ⎥ ⎢ − 2 cos θgo ⎥ ⎢ ⎪ ⎪ ⎪ ⎪ ⎥ ⎢ r ⎪ ⎪ ⎪ ⎦ ⎪ cos θgo + αi ⎣ ⎭ ⎩ − r go sin(αi ) cusp cusp + σgw × L gw r go , θow,i , f bow (38) r +α cos θgo i π cusp 2 r r θgo C3 r go , θow,i , f bow = r go + αi − + cos θgo 2 sin(αi ) cusp cusp cusp − Acusp (39) r go , θow,i , f bow − Ao r go , θow,i , f bow w ⎡
Once the geometry of a cusp has been fixed (through minimization of the phase angle error, see Eq. 11), energy balance equation (Eq. 36) is invoked to investigate the convergency of r go by solving the system of nonlinear equations written for cusp formation. The equations for all the other possible two- and three-phase displacements in capillary tubes with irregular triangular cross sections are listed in the Appendices B and C.
123
A. Zolfaghari, M. Piri
4.1 Two-Phase Threshold Capillary Pressures In this section, we present the thermodynamically consistent two-phase threshold capillary pressures of displacements in capillary tubes with irregular triangular cross sections. We study the impacts of parameters such as advancing and receding contact angles, capillary pressure ratio, and shape factor. In the interest of brevity, in this section, we include only two sets of results with a summary of our findings about the effect of different parameters on the stability of wetting layers. To generate the results presented here, we used 48 mN/m and 13.6 μm for oil/water interfacial tension and inscribed radius, respectively. Inscribed radius is the average inscribed radii of the pores and throats in the Berea sandstone network studied by Lopez et al. (2003), Valvatne and Blunt (2004), and Piri and Blunt (2005a, b). The interfacial tension is that of distilled water and n-hexane reported by Firincioglu et al. (1999). Here, we investigate which of the displacement routes are the most favorable for water invasion into oil under different wetting conditions and capillary pressure histories. We have performed an extensive study on several thousand randomly generated irregular triangles (Zolfaghari 2014), but here only one of the most interesting cases is presented to illustrate the MSP analysis for the multi-layer formation and collapse displacements. The results presented in this section are divided into two main categories. In the first category, we analyze PL displacement of oil by water, while in the second, we study the collapse of wetting oil layers that may form due to PL displacements discussed in the first category. The initial condition for the first group is the fluid configuration at the end of primary drainage (see configuration A in Fig. 1). The initial condition for the second group, however, is one of the fluid configurations with oil layer(s) sandwiched between water in the corner and water in the center. We not only present the threshold capillary pressures at which the layers collapse, but also the sequence by which they occur. The analysis in both categories is based on the MSP theory. Furthermore, we, for comparison, include threshold capillary pressures of oil LC obtained from the geometric criterion and snap-off (see Ref. Zolfaghari (2014) for pertinent equations). During waterflooding, the displacement with the highest threshold oil/water capillary pressure is favored, if it is relevant. Below we present the results for each category.
Piston-like displacements: Here we focus on weakly and strongly oil-wet systems and oil layer formation due to PL and snap-off displacements. Figure 4a shows the threshold capillary pressure results in a triangular capillary tube with corner half angles of 25.296◦ , 30.688◦ , and 34.016◦ . The capillary pressure ratio (the ratio of oil/water capillary pressure at the end of oil primary drainage to the threshold oil/water capillary pressure needed for the same drainage displacement) is assumed to be ten for this case to enhance the possibility of oil a ≥ π/2 + α requirement layer formations. Three vertical lines in the graph indicate the θow i for the existence of oil layers in each corner. This divides the graphs into four regions. In each region, we present only the threshold capillary pressure of the relevant displacements. For example, any displacements related to configurations with a layer in the wide corner is a < π + α , where α is the corner not relevant in regions of advancing contact angle with θow 3 3 2 half angle of the wide corner. Nine different displacement indices (0–8) have been assigned to different configuration changes. Displacement index (DI) 0 indicates PL displacement with no wetting oil layer in any corner. The next three indices, i.e., 1–3, are used for PL displacements with 1LF in sharp, medium, and wide corners, respectively. Similarly, 2LF could occur in three different pairs of corners, i.e., sharp/medium (DI 4), sharp/wide (DI 5), and medium/wide (DI 6).
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
(a)
(b)
Fig. 4 Two-phase threshold oil/water capillary pressures of the piston-like and snap-off displacements calculated for an irregular triangular capillary tube. a The impact of advancing oil/water contact angle, and b its most favorable displacement(s)
PL displacement with 3LF is represented by DI 7. DI 8 stands for forced snap-off (see Ref. Zolfaghari 2014). As it is observed from Fig. 4a, PL displacement does not lead to formation of oil layers in any corners at lower oil/water advancing contact angle (configuration changes A → C, DI 0). This is attributed to the low advancing oil/water contact angle and relatively wide corner a = 150o , PL displacements angles of the triangle. At higher advancing contact angles, e.g., θow forms an oil layer only in the sharp corner (DI 1). For a given capillary pressure ratio, the minimum advancing contact angle over which DI 1 is favored depends on the sharpness of the corner angle. If the corner angle is small enough, it could form an oil layer for the entire a ≥ π/2 + α , where α is the corner half range of advancing contact angle greater than θow 1 1 angle of the sharp corner. This means that at low capillary pressure ratios, PL displacement does not leave oil layers in the corners unless the corners are extremely sharp. a = 165◦ , DI 4 becomes the most For higher advancing oil/water contact angle, e.g., θow favorable displacement. Displacement 4 corresponds to PL displacement with 2LF in sharp and medium corners. Corner half angles for the sharp and medium corners are relatively similar and hence, as the system becomes more oil-wet, PL displacement with layer formation in two corners becomes favorable over the displacement with layer formation only in the sharp corner. It is also possible not to have displacement 1 through PL displacement at lower advancing contact angles, but instead form 2LF in the sharp and medium corners. This could happen when two corner half angles for the sharp and medium corners are almost the same. The results show that the smaller is the medium corner, the greater is the chance of simultaneous oil layer formation in the sharp and medium corners. This behavior is even more likely under stronger oil-wet conditions. If both corner angles become very small, PL displacement forms an oil layer in both sharp and medium corners for the entire range of
123
A. Zolfaghari, M. Piri
advancing contact angles greater than π/2 + α2 , where α2 is the corner half angle of the medium corner. If two corners become much smaller than the third one, the possibility of oil layer formation in both corners increases particularly at higher advancing contact angles. At intermediate values of advancing contact angle, PL displacement may form an oil layer only in the sharp corner. For the case of two corners having very close half angles, oil layer formation in both corners is the only scenario that takes place at higher advancing contact angles. a = 175◦ , PL displacement leaves an oil layer in For the strongly oil-wet systems, e.g., θow all corners of the capillary tube. The corner half angles of the triangle are not significantly enough different from each other, explaining the reason for configuration change A → B at high advancing contact angles. A comprehensive study has been carried out on the subject of formation and collapse of wetting oil layers in pores with irregular triangular cross sections (Zolfaghari 2014). The corresponding trends have been implemented in our network model (Zolfaghari and Piri 2016). Below, we include a summary of the observed trends:
(1) If a PL displacement forms one oil layer, the layer always forms in the sharp corner, i.e., among displacements 1, 2, and 3, displacement 1 is favored. (2) If a PL displacement forms two oil layers, the layers will always form in the sharp and medium corners, i.e., among displacements 4, 5, and 6, displacement 4 is favored. (3) In triangles with considerably different corner half angles, three-layer formation is only favored if two-layer formation and one-layer formation are favored in lower advancing contact angles (see Fig. 4b). This is also true for PL displacement with two-layer formation, i.e., it takes place only if one-layer formation is favored in lower advancing contact angles. If two corner angles are similar, two-layer formation will be favored without any one-layer formation event taking place at lower advancing contact angles. If the three corner half angles are similar, PL displacement with three-oil-layer formation and PL displacement without oil layer formation are the only possible scenarios. (4) Similar to the displacement in capillary tubes with regular cross sections, as the capillary pressure ratio increases, the possibility of PL displacements with oil layer formation becoming the most favorable displacement, enhances. Higher capillary pressure ratios can make displacement 0 unfavorable for the entire range of advancing contact angle greater than π2 + α1 . (5) The possibility of oil layer formation increases in triangles with low shape factors. This is because the corner half angles of the sharp and medium corners are smaller in low shape factor triangles than those of triangles with higher shape factors. Layer collapse displacements: In this section, we consider the second category, i.e., the collapse of wetting oil layers formed by, for instance, PL displacement. In all examples presented here, the MSP results are included together with the oil LC threshold capillary pressures calculated using the geometric criterion (see Ref. Zolfaghari 2014). We use 10 displacement indices (9–18) in this subsection to cover all possible displacement scenarios. The threshold capillary pressure of displacements 9–15 are obtained using MSP method while those of displacements 16–18 are found from the geometric criterion. Displacement indices of 9, 10, and 11 correspond to single layer collapse events in sharp, medium, and wide corners, respectively. Similarly, simultaneous collapse of two layers could occur in three different pair of corners, i.e., sharp/medium (DI 12), sharp/wide (DI 13), and medium/wide (DI 14). DI 15 stands for simultaneous collapse of three layers in all corners. Displacements 16–18 represent single LC events in different corners of an irregular triangle similar to displacements 9–11.
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
(a)
(b)
Fig. 5 Two-phase threshold oil/water capillary pressures of the layer collapse displacements calculated for an irregular triangular capillary tube. a The impact of advancing oil/water contact angle, and b its most favorable displacement(s)
We analyze LC displacements for the layers formed in Fig. 4. We find the most favorable LC event for each example at different values of advancing contact angle. As we have shown in Fig. 4b, displacements 0, 1, 4, and 7 are favored at different advancing contact angles for the PL invasion of water into oil. The latter three hold oil layers in one, two, and three corners, respectively. Further reduction in capillary pressure will ultimately collapse these layers forming configuration C. Figure 5a presents the threshold capillary pressures associated with those collapse scenarios. Figure 5b shows the events that are the most favored for different ranges of advancing contact angles. For example, PL displacement causes 3LF (through displacement 7) at advancing contact angles greater than 168.82◦ (see Fig. 4b). Decreasing capillary pressure collapses the oil layer in the wide corner first (displacement 11, the most favorable displacement in this range of advancing contact angle). If we continue reducing capillary pressure, the capillary tube will lose the oil layer in the medium corner as the displacement 10 is the next most favorable displacement (see Fig. 5a). One should note that displacement 14, that has threshold capillary pressures greater than those of displacement 10, is irrelevant because it represents simultaneous collapse of oil layers in the medium and wide corners while capillary tube does not have an oil layer in the wide corner, i.e., it was collapsed in the previous step. And finally, the oil layer in the sharp corner collapses (displacement 9), as it is the next most favorable displacement. Similarly, displacements 15, 13, and 12 that have higher threshold capillary pressures than those of displacement 9 are irrelevant because they represent simultaneous collapse of two oil layers while capillary tube has only one layer at this stage. For advancing contact angles between 159.08◦ and 168.82◦ , PL displacement forms two oil layers in sharp and medium corners. The collapse of oil layers in this fluid configuration follows the same sequence discussed above. Piston-like displacement forms only one layer at advancing contact angles ranging from 145.445◦ to 159.08◦ (see Fig. 4b). Further reduction in capillary pressure will lead to collapse of the only layer available through displacement 9 and formation of configuration C (see Fig. 5a).
123
A. Zolfaghari, M. Piri Table 1 Interfacial tensions in mN/m used in Sects. 4.1 and 4.2 (Hui and Blunt 2000; Firincioglu et al. 1999) Fluids
σow
σgo
σgw
Hexane-water-air
48
20
67
It is important to note that the threshold capillary pressures obtained from the geometric criterion for the collapse of the above-mentioned oil layers are significantly lower than those obtained from the MSP method. This means that the layers are stable for shorter span of capillary pressure, which will in turn affect, for instance, the extent to which the residual oil saturation can be reduced during capillary-dominated waterflooding into weakly and strongly oil-wet systems. It may also affect the rate by which the oil relative permeability reduces at low oil saturations in those systems. Our study shows that the threshold capillary pressure associated with oil LC events reduces with the increase in capillary pressure ratio. In other words, the oil layers are stable for a wider span of capillary pressure, which could potentially keep oil, in an oil-wet porous medium, connected at low oil saturations. Furthermore, it is evident that only single oil LC events, as apposed to simultaneous LC events, are favored and take place in corners with a decreasing size of corner half angles, i.e., wide, medium, and sharp corners. This, however, was not the case for the formation of oil layers due to PL displacement, i.e., simultaneous formation of two or three layers were indeed possible and were favored under different scenarios.
4.2 Three-Phase Threshold Capillary Pressures In this section, we study thermodynamically consistent threshold capillary pressures of a three-phase system. Two examples are presented here for gas injection into water and oil in two different capillary tubes with regular and irregular triangular cross sections. The inscribed radius is 13.6 μm for both capillary tubes. We analyze the impact of initial oil/water capillary pressure on the evolution of fluid configurations due to gas injection. The oil drainage capillary pressure ratio and spreading coefficient are 10 and −1 mN/m, respectively. Table 1 lists the three-phase interfacial tensions used in this section. The results are obtained for a weakly water-wet system with receding and advancing oil/water contact angles of 35◦ and 55◦ , respectively. Receding and advancing gas/oil contact angles are 0◦ . Using Bartell and Osterhof (1927) constraint, the receding and advancing gas/water contact angles are calculated to be 27.7◦ and 44.8◦ , respectively. Figures 6 and 7 show the calculated threshold gas/oil capillary pressures for a wide range of initial oil/water capillary pressures for the regular and irregular triangular cross sections, respectively. The x-axes cover the entire range of capillary pressures bounded by oil/water snap-off capillary pressures and the maximum capillary pressures reached during drainages with the capillary pressure ratios of 10. They indicate different extents of waterflooding in the absence of a water MTM in the capillary tubes. This is due to the water invasion process prior to gas injection. Threshold capillary pressures for all relevant displacement scenarios are presented in Fig. 6 for a regular triangular capillary tube. Since corner angles are identical (60◦ ), only displacements with simultaneous three LF, LC, cusp formation, and no LF are included. Other displacements with one or two LF or LC are irreverent in this case. Two geometric constraints proposed by Fenwick and Blunt (1998b) and Hui and Blunt (2000) are also plotted in Fig. 6 to define capillary pressure regions over which spreading oil layers can form. They are called b- and peak constraints in this work. b is the meniscus-apex distance on the pore wall.
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
Fig. 6 Three-phase threshold gas/oil capillary pressures of the layer formation, collapse, and cusp displacements calculated for a regular triangular capillary tube with an oil drainage capillary pressure ratio of 10
b-constraint relates to a situation where both interfaces surrounding an oil layer have similar b values. Peak constraint, however, determines threshold capillary pressure at which oil/water and gas/oil interfaces touch each other in the middle. During gas injection, gas/oil capillary pressure increases (at constant oil/water capillary pressure), and hence the displacement with the lowest threshold gas/oil capillary pressure is the most favorable, if it is relevant. At high oil/water capillary pressures, the most favorable displacement is the layer formation in all corners (configuration changes A → D, see Fig. 1), as it has the lowest threshold capillary pressures. Moving upward at constant oil/water capillary pressure, the next displacement is no LF (configuration changes A → F). This is an irrelevant displacement as configuration A has already changed to D. Further increase in the gas pressure collapses these layers into cusp interfaces (as indicated by the cusp plot in Fig. 6). Cusps shrink as gas/oil capillary pressure increases further. Note that cusps form at much lower threshold capillary pressures than those of the LC displacements (compare cusp and 3LC plots). Without cusp analysis, spreading layers would exist over a wider span of gas pressures. That impacts the predicted oil relative permeabilities and residual saturations during gas injections (Zolfaghari and Piri 2016). Later in this section, we present more results on layer and cusp shrinkage/collapse. At lower oil/water capillary pressures, cusp forms in all three corners without prior layer formation. At very low oil/water capillary pressure, closer to the oil/water snap-off value, three curves of 3LF, 0LF, and 3LC cross at one point beyond which no LF is the only legitimate displacement. That provides a consistent transition of pore fluid occupancies as a function of oil/water capillary pressure. High water pressures in this region do not allow spreading layers to form. This is evident by the small values of gas/oil capillary pressures for the two geometric constraints (see b- and peak constraints in Fig. 6 at very low Pcow values). Note that layers cannot form above geometric criteria marked by b- and peak-constraint curves. Figure 7 illustrates three-phase threshold capillary pressures in an irregular triangle with the corner half angles of 10◦ , 30◦ , and 50◦ . For simplicity, only the most favorable displacements of one LF category (i.e., LF in the sharp corner), and two LF category (i.e., LF in the sharp and medium corners) are plotted. We have included two geometric criteria of LC only for the wide corner. The geometric criterion for each corner is defined as the minimum
123
A. Zolfaghari, M. Piri
Fig. 7 Three-phase threshold gas/oil capillary pressures of the layer formation, collapse and cusp displacements calculated for an irregular triangular capillary tube with an oil drainage capillary pressure ratio of 10
threshold gas/oil capillary pressure between b- and peak constraints. LC geometric criteria for the sharp and medium corners are excluded from Fig. 7. Once again gas pressure increases during gas injection, and as a result, displacements with lower threshold capillary pressures are more favorable. For the results presented in Fig. 7, PL displacement with three LF is the most favorable displacement at high oil/water capillary pressures. As the oil/water capillary pressure reduces, the possibility of cusp formation without a prior layer formation increases. This is consistent with the results of regular capillary tube presented earlier. Even though spreading coefficient is negative, gas injection leaves spreading oil layers or cusps in all three corners of capillary tube. This is consistent with several experimental observations reported in the literature. Layer drainage has been observed experimentally for non-spreading systems during gravity drainage (Sahni et al. 1998; DiCarlo et al. 2000; Hayden and Voice 1993; Zhou and Blunt 1997) and tertiary gas injection (Øren and Pinczewski 1992). Micro-model experiments have also shown the possibility of layer formation in systems with negative spreading coefficients (Dong et al. 1995; Keller et al. 1997). These layers, in essence, are different than those discussed in Sect. 4.1. Formation of oil layers in three-phase systems is highly dependent on the spreading characteristics of the intermediate wetting phase in the presence of the most wetting and non-wetting phases. The above-mentioned parameters also impact the stability of spreading oil layers in threephase systems (i.e., their collapse threshold capillary pressures). Similar to wetting layers, collapse of spreading oil layers could take place through a displacement with an MTM that can be analyzed using MSP method or geometric constraints. For systems with negative spreading coefficient, a PL displacement could collapse a layer completely, or collapse it partially to form a cusp configuration. For systems with positive spreading coefficients, cusp configurations cannot form (see Sect. 3.2), and hence, oil layers could only collapse completely and form gas/water interface. In addition to MSP LC displacements, layers could collapse based on a geometric criterion regardless of the value of spreading coefficient. As
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
mentioned earlier, PL displacement with layer formation in all corners is the most favorable displacement since it has the lowest gas/oil capillary pressure among all possible displacement scenarios shown in Fig. 7. Moving up in this figure at high oil/water capillary pressures, the next set of displacements with lowest gas/oil capillary pressures are displacements with two, one, and no LF. These are all irrelevant displacements since three LF displacement has already taken place at a lower gas/oil capillary pressure. Further increase in the gas/oil capillary pressure causes the layer in the widest corner to collapse into a cusp configuration. The oil layers in the medium and sharp corners follow similar path by collapsing into cusp configurations at higher gas/oil capillary pressures indicated by the cusp (medium) and cusp (sharp) threshold capillary pressure curves in Fig. 7, respectively. One should note that, if cusp formation threshold capillary pressure for a particular corner falls under that of three LF displacement, that corner would have a cusp configuration without prior formation of a spreading oil layer. For example, final configuration of a PL displacement with three LF at oil/water capillary pressure of 14000 Pa in Fig. 7 would be oil layers in the sharp and medium corners and cusp configuration in the wide corner. Cusp configurations eventually shrink to configuration F (see Fig. 1) as gas/oil capillary pressure increases. Our model can handle cusp shrinkages due to increase in gas/oil capillary pressure or decrease in oil/water capillary pressure. All other LC displacements depicted in Fig. 7 are irrelevant, since all oil layers have been collapsed into cusp configurations at lower gas/oil threshold capillary pressures. As it is seen, cusp formations are always favored for the system studied here since their threshold gas pressures are less than those of complete LC events. An additional equation for water area is added to the system of nonlinear equations for each corner with a new gas/water interface for all relevant displacements. These are the displacements in which at least one new gas/water interface is created, i.e., 0LF, 1LF, 2LF, 1LC, 2LC, and 3LC displacements. Contact angle of the newly created gas/water interface is treated as an unknown parameter for each individual corner. As explained in Sect. 3, for each of the displacements with the newly created gas/water interface(s), the threshold values are obtained through an optimization process for the solution of system of nonlinear equations. To illustrate the shrinkage processes, Fig. 8 shows corner fluid configurations with cusps and layers at two different gas/oil capillary pressures. The shrinkage is caused by gas injection into a regular capillary tube of Fig. 6 with and without cusp analysis. Increasing gas/oil capillary pressure from 3124 Pa (Fig. 8a, b) to 4000 Pa (Fig. 8c, d) causes oil areas of the cusp and layer to reduce by 71 and 42%, respectively. This indicates that cusps shrink with a faster rate. Furthermore, cusps hold significantly smaller oil areas in the cross section (compare corner configuration (a) with (b) and (c) with (d) in Fig. 8). To investigate this further, Fig. 9 is constructed to illustrate variation of oil areas for the cusp and layer shown in Fig. 8a, b. Included in this Figure are also MSP and geometric layer collapse threshold capillary pressures consistent with the presented results in Fig. 6 (see 3LC and peak constraint curves). Without cusp analysis, model predicts formation of more conductive oil layers for a given set of capillary pressures (as illustrated in Fig. 8). Cusp configuration significantly reduces oil area in the cross section while maintaining its phase connectivity (see Fig. 9). This behavior impacts conductance of the oil phase (i.e., oil relative permeability) and its residual saturations during gas injection. Cusp configurations eventually shrink to no layer configuration when oil area diminishes to zero. As we have shown in Zolfaghari and Piri (2016), since cusp cross-sectional area is smaller than that of a full layer, incorporation of cusp analysis into the model results in lower oil relative permeability curves while maintaining the oil connectivity throughout the network.
123
A. Zolfaghari, M. Piri
(a)
(b)
(c)
(d)
Fig. 8 Simulated corner configurations in a regular triangular capillary tube used in Fig. 6: a, c with cusp analysis; b, d without cusp analysis. Corner configurations a, b represent oil/water and gas/oil capillary pressures of 4000 and 3124 Pa, respectively. Corner configurations c, d demonstrate cusp and layer shrinkages due to gas injection (Pcgo = 4000 Pa). In this figure blue, red, and green represent calculated oil/water, gas/oil, and gas/water interfaces, respectively
5 Conclusions A state-of-the-art three-phase pore-scale network model was presented. The model incorporates: (1) thermodynamically consistent threshold capillary pressures for different displacements including those that involve oil layer formations and collapses, (2) a new pore fluid configuration with oil cusps, (3) new configuration change routes that allow gas injection into oil producing configurations with oil layers, with oil cusps, or without oil, (4) a new saturation path tracking algorithm, (5) a new set of different layer collapse rules, and (6) a new methodology based on the MSP method to handle newly formed gas/water interfaces in order to eliminate inconsistencies in relation between capillary pressures and pore fluid occupancies. The model’s capabilities were expanded in order to improve its predictive capabilities particularly at low oil saturations where previously developed models and empirical correlations fail to accurately predict flow properties such as relative permeabilities. The model handles the formation and collapse of the wetting and spreading oil layers to capture physics related
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
Fig. 9 Cusp and oil layer shrinkages as a result of a gas injection process
to the experimental observations at low oil saturations during two- and three-phase flow tests. We used MSP method to calculate thermodynamically consistent threshold capillary pressures of piston-like and layer collapse displacements during water or gas flooding. A new fluid configuration involving oil cusps was introduced for three-phase systems with negative spreading coefficients. The formation and collapse of wetting and spreading oil layers were extensively investigated by solving the appropriate thermodynamically consistent threshold capillary pressure equations for capillary tubes with different angular, irregular cross sections. The results of thermodynamically consistent threshold capillary pressures show that spreading oil layers can exist for the non-spreading systems. Cusp analysis significantly reduces oil volume in the pore space while maintaining connectivity of the oil phase. As a result, it impacts oil relative permeabilities and residual saturations particularly at the low oil saturation regions. Acknowledgements We gratefully acknowledge financial support of Total, Saudi Aramco, EnCana, and the School of Energy Resources and the Enhanced Oil Recovery Institute at the University of Wyoming. We extend our gratitude to Hu Dong and Sven Roth of iRock Technologies for generating some of the pore networks for our samples. We thank Pål-Eric Øren of FEI for sharing his Berea network data. We also thank Henry Plancher and Soheil Saraji of Piri Research Group for measuring the interfacial tension and density values and for the details of the IFT experimental procedure used.
Appendix A Cusp Interfacial Contact Lengths and Phase Areas Fluid–Fluid and Fluid–Solid Contact Length Calculations For an interface, the length of an arch can be determined by knowing its two end points and radius. Oil/water and gas/oil interfaces are stretched between P3P contact and bow , and P3P contact and bgo , respectively (see Fig. 2). By using Eq. (3) for P3P’s coordinates, the oil/water and gas/oil contact lengths in a given cusp configuration can be formulated as: ⎫ ⎧ * + ⎪ ⎬ ⎨ + [bcontact cos (α) + H sin (α) − b12 cos (α)]2 ⎪ (40) L 12 = 2 r12 × arcsin (r12 )−1 , ⎪ + [bcontact sin (α) − H cos (α) − b12 sin (α)]2 ⎪ ⎭ ⎩
123
A. Zolfaghari, M. Piri
where subscript 12 represents ow or go interfaces, bcontact is the b-value of the perpendicular line passing through P3P contact to the pore wall, and H is the cusp height. All parameters are clearly shown in Fig. 2. The multiplication by 2 in Eq. (40) is to account for two mirrorimaged interfaces between similar fluids in each corner. The gas/water contact length in a cusp is calculated by: bcontact sin (α) − H cos (α) L gw = 2 r gw × arcsin (41) r gw where r gw is the gas/water curvature. The oil/solid contact length is simply the difference between bgo and bow accounting for cusp configurations on both pore walls. It is expressed as: L os = 2 bgo − bow (42) Since pore’s perimeter is constant, by having oil/solid and water/solid (L ws = 2 bow ) contact lengths, the gas/solid contact length can be calculated by: L gs = Pt − L os − L ws = Pt − 2 bgo
(43)
where Pt is the perimeter of the pore.
Pore Fluid Area Calculations In this section, we present equations for the calculation of pore cross-sectional areas occupied by water, oil, and gas. To calculate different phase areas, all the interfaces and points of contact should be mathematically defined for each corner of an angular pore. Oil area in a cusp configuration is calculated using:
-
max(x P3P , xbgo )
Ao = 2 ×
2 2 dx ycgo + f go r go − x − xcgo
min(x P3P , xbgo )
-
max(x P3P , xbow )
−
2 2 − (x − x ycow + f ow row dx ) cow
(44)
min(x P3P , xbow ) xbgo
−
x dx + eow − ego tan (α)
xbow
One should note that each fluid–fluid interface on the pore cross section is modeled as part of a circle. f go and f ow in Eq. (44) are used to determine the sections to which the gas/oil and oil/water interfaces belong. They could hold −1 and 1 values representing the upper and lower semicircles, respectively. eow and ego are the correcting terms for oil/water and gas/oil interfaces, respectively, when they are part of both lower and upper semicircles, otherwise, they are zero. Their values are calculated geometrically for any interface in a given cusp configuration that requires integration over two semicircles. The first and second terms in Eq. (44) are integrals of the gas/oil and oil/water interfaces, respectively. It is formulated
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
Fig. 10 Cross-sectional areas representing various terms given in Eq. (47) developed for calculation of gas area within the selected region. A g, base is shown by the dotted black area. Areas shaded by the green squares and red lines represent the second term, and the summation of third and fourth terms on the right-hand side of Eq. (47), respectively
based on the premise that gas/oil interface is below oil/water interface between x P3P and xbgo (see Fig. 2). Otherwise, gas/oil interface intersects oil/water interface for x > x P3P , and hence, cusp configuration is not geometrically possible. Equation (44) is used only for geometrically feasible cusp configuration. Another criterion is to have a positive oil/solid contact length (i.e., xbow < xbgo ). This is an important assumption for writing the third term in Eq. (44), which is the integral over the pore line between two intersections of oil/water and gas/oil interfaces (i.e., xbow and xbgo ). For each corner with cusp configuration, there are two cusps that are mirror image of one another with respect to the middle line. And therefore, a multiplication by 2 is used in Eq. (44). Having three fluid phases, i.e., oil, water and gas, in a pore with a constant cross-sectional area requires calculation of two of the fluid phase areas. We select oil and water areas. Although the same concept of integration over two different interfaces could be used for the calculation of water area [as was used for the case of oil area in Eq. (44)], it would be complicated to present a general equation to cover all possible interface curvatures and P3P locations. Instead, gas area is calculated in a predetermined site, called base area, covering the entire cusp configuration in a corner. The water area is then calculated as the total base area minus gas and oil areas within the selected base. Below, we provide more details.
123
A. Zolfaghari, M. Piri
The base area is chosen to encompass the entire cusp configuration. A line with the equation of y = ybase is used to create a triangle by intersecting the pore wall and the middle line (see Fig. 10). The triangle creates the base site for the calculation of water area. The value of ybase is defined by: max ycgw , ybgo , y P3P if r gw > 0 ybase = (45) max ycgw + |r gw |, ybgo otherwise As an example, see Fig. 10, where the base line is drawn for a cusp configuration with positive gas/water capillary pressure, and gas/water center being beneath the P3P contact and gas/oil center. The total area of the base is written as: 1 2 At, base = ybase tan (α) (46) 2 Gas area within the base region (A g,base ), black dotted area in Fig. 10, can be written as the total area of the base minus the areas occupied by water and oil: A g, base = 2 ×
max x P3P , xbgo
1 x 2P3P − At, base − 2 tan (α)
max x P3P , xbgo
-
+ min x P3P , xbgo
-
2 2 − x−x ycgo + c f go r go dx cgo
min x P3P , xbgo
x x P3P dx + bulgego − x P3P y P3P − + A g,int tan (α) tan (α)
(47) where At, base is calculated using Eq. (46). The second term on the right-hand side excludes a triangle at the top of the base area, which is not covered by the third and fourth integrals (see the triangular area depicted by green squares in Fig. 10). The third term is the integral of gas/oil interface between x P3P and xbgo , similar to the first term in Eq. (44). The fourth term is the pore wall’s integral. It covers the same integral bounds as of the previous term. Together, the third and fourth terms deduct the area shaded by the red lines in Fig. 10 from the base area. The sixth term removes an area of a rectangle over the gas/water interface where is mainly covered by the water phase. The last term in Eq. (47) corrects for the gas area that is removed by the rectangle in the previous term. Its value depends on the sign of gas/water curvature and the position of P3P contact with respect to the center of gas/water interface. For a positive curvature, i.e., Pcgw > 0, it is given by:
A g,int
⎧ 2 r x2 x P3P ⎪ ⎪ ⎪ gw × arcsin P3P if y P3P < ycgw − ⎪ ⎪ x P3P 2 r gw ⎪ ⎪ 2 × tan arcsin ⎪ ⎪ ⎪ r gw ⎪ ⎪ 2 ⎪ ⎪ r gw ⎪ ⎪ ×π if y P3P = ycgw ⎪ ⎪ ⎨ 4 = 2 ⎪ r gw y P3P − ycgw π ⎪ ⎪ + × arcsin ⎪ ⎪ ⎪ 2 r gw 2 ⎪ ⎪ ⎪ ⎪ 2 ⎪ if y P3P > ycgw y P3P − ycgw ⎪ ⎪ ⎪ + ⎪ ⎪ y − y ⎪ P3P cgw ⎪ ⎩ 2 × tan arcsin r gw
123
(48)
Pore-Scale Network Modeling of Three-Phase Flow Based on…
By using oil area from Eq. (44), total base area from Eq. (46), and gas area within the base from Eq. (47), water area in the cusp configuration is calculated by: Aw = At, base − A g, base − Ao
(49)
Gas area is then calculated by knowing oil and water areas as well as the total crosssectional area of the pore. At this stage, all the required parameters are determined for invoking the MSP analysis of a cusp formation displacement.
Appendix B Thermodynamically Consistent Threshold Capillary Pressures for Two-phase Flow Displacements in Capillary Tubes With Regular Angular Cross Sections In this section, we present the final equations pertaining to various PL and LC displacements in capillary tubes with regular angular cross sections. They are derived using the MSP method. Since capillary tubes with regular cross sections have the same corner half angels and hinging contact angles in all corners, final threshold capillary pressure equations are relatively simple compared to their counterparts in irregular triangles. One may have three different displacements in regular capillary tubes: – PL displacement with the formation of wetting oil layers (configuration changes A → B in Fig. 1): a π 2 a a cos(θow − α) a 3 row ) + At = 0 (50) + α − θow + cos(θow ) − row L t cos(θow 2 sin(α) – PL displacement without formation of any wetting oil layers (configuration changes A → C in Fig. 1): π a 2 h − α − θow, row L t cos(θow ) − At + 3 row 1 2 h cos(θ ow,1 + α) h a =0 (51) + cos(θow, ) − 2 cos(θ ) ow 1 sin(α) h One should note that subscript 1 in θow, refers to interface one (counted from apex 1 toward center). We do not specify the corner number because all the corners have similar corner half angles. Equations (50) and (51) can be used for the similar displacements in square cross-section pores by replacing the coefficient 3 with 4. Note that α = π/4 for squares. – PL displacement with the collapse of wetting oil layers (configuration change B → C in Fig. 1): a − α) cos(θow sin(α) h cos(θ ow,1 + α) a h =0 ) − cos(θow, ) + 2 cos(θow 1 sin(α)
a h a θow + θow, − π − cos(θow ) 1
(52)
123
A. Zolfaghari, M. Piri
Displacements in Capillary Tubes With Irregular Triangular Cross Sections Here we list the threshold capillary pressure equations for various displacements in capillary tubes that are irregular triangle in cross section. They are divided into two groups: 1) PL displacement of oil by water with and without oil layer formation and 2) oil LC displacements.
5.0.1 PL Displacements With and Without Oil Layer Formation – PL displacement with the formation of one wetting oil layer in corner 1: a 2 h At − row L t cos(θow ) + row θow, 2 ,1 ξ1 + ξ2 + ξ3 h a + θow, − θ + =0 ow 3 ,1 sin(α1 ) sin(α2 ) sin(α3 )
(53)
where a a ) sin(α2 ) sin(α3 ) cos(θow − α1 ) ξ1 = cos(θow h a h ξ2 = sin(α1 ) sin(α3 ) cos(θow,2 ,1 + α2 ) 2 cos(θow ) − cos(θow, ) 2 ,1 h a h ξ3 = sin(α1 ) sin(α2 ) cos(θow, + α3 ) 2 cos(θow ) − cos(θow, ) 3 ,1 3 ,1
(54) (55) (56)
Subscript 1 does not necessarily refer to the sharp corner. It represent any corner in which the oil layer forms. – PL displacement with the formation of two wetting oil layers in corners 1 and 2: a 2 h a At − row L t cos(θow ) + row − 2θow π + θow, 3 ,1 ξ1 + ξ2 + ξ3 =0 (57) + sin(α1 ) sin(α2 ) sin(α3 ) where a a ξ1 = cos(θow ) sin(α2 ) sin(α3 ) cos(θow − α1 )
ξ2 = ξ3 =
a a cos(θow ) sin(α1 ) sin(α3 ) cos(θow − α2 ) h a h sin(α1 ) sin(α2 ) cos(θow,3 ,1 + α3 ) 2 cos(θow ) − cos(θow, ) 3 ,1
(58) (59) (60)
Subscripts 1 and 2 do not necessarily refer to the sharp and medium corners. They represent any pair of corners in which the oil layers form. – PL displacement with the formation of three wetting oil layers (configuration changes A → B in Fig. 1): see Eqs. (21)–(24) in Sect. 4. – PL displacement without formation of any wetting oil layers: see Eqs. (17)–(20) in Sect. 4.
Layer Collapse Events – PL displacement with the simultaneous collapse of two wetting oil layers: a a a − cos(θow ) sin(θow ) −π + θow a h − cos(θow ) sin(θow, )+ 1 ,1
123
1 h h θ + θow, 2 ,1 2 ow,1 ,1
Pore-Scale Network Modeling of Three-Phase Flow Based on…
h h h h + cos(θow, ) sin(θow, ) + cos(θow, ) sin(θow, ) 1 ,1 1 ,1 2 ,1 2 ,1 a h ) sin(θow, )− − cos(θow 2 ,1
−
2 a ) − cos(θ h cos(α1 ) cos(θow ow,1 ,1 )
2 a ) − cos(θ h cos(α2 ) cos(θow ow,2 ,1 ) 2 sin(α2 )
2 sin(α1 ) =0
(61)
Subscripts 1 and 2 do not necessarily refer to the wide and medium corners. They represent any pair of corners in which the oil layers collapse. – PL displacement with the simultaneous collapse of three wetting oil layers (configuration changes B → C in Fig. 1): a a a −π + θow − cos(θow ) sin(θow ) 2 a h h h − cos(θow ) sin(θow, ) + sin(θ ) + sin(θ ) ow,2 ,1 ow,3 ,1 1 ,1 3 1 h h h h h + θow, + θow, + θow, + cos(θow, ) sin(θow, ) 1 ,1 2 ,1 3 ,1 1 ,1 1 ,1 3 h h h h ) sin(θow, ) + cos(θow, ) sin(θow, ) + cos(θow, 2 ,1 2 ,1 3 ,1 3 ,1
− − −
2 a ) − cos(θ h cos(α1 ) cos(θow ow,1 ,1 ) 3 sin(α1 ) 2 a ) − cos(θ h cos(α2 ) cos(θow ow,2 ,1 ) 3 sin(α2 ) 2 a ) − cos(θ h cos(α3 ) cos(θow ow,3 ,1 ) 3 sin(α3 )
=0
(62)
– PL displacement with the single collapse of a wetting oil layer: see Eq. (25) in Sect. 4.
Appendix C Thermodynamically Consistent Threshold Capillary Pressures for Three-phase Flow Threshold capillary pressure equations for various three-phase displacements in irregular triangular cross sections are presented in this section. They are presented in two main groups: 1) PL displacement of oil by gas with and without oil layer formation and 2) oil LC displacements.
PL Displacements With and Without Oil Layer Formation – PL displacement with the formation of one spreading oil layer in corner 1: 3 2 r go C1 + r go C2 + r go C3 + C4 = 0
(63)
where C1 –C4 are given by: C1 =
cos(θ r + α1 ) σow σgo π go r r − cos θgo − α1 − θgo row 2 sin(α1 )
(64)
123
A. Zolfaghari, M. Piri
⎤
Aw,2 + Aw,3 + 2 row θow,2 + θow,3 + α2 + α3 − π ⎥ row σow ⎢ ⎢ ⎥ C2 = ⎢ ⎥ cos θow,3 + α3 cos θow,2 + α2 ⎦ row ⎣ r + σgo cos θgo + L t − 2 row sin (α2 ) sin (α3 ) cos(θ r + α1 ) go 2 π r r − cos θgo + σgo − α1 − θgo 2 sin(α1 ) ⎡ r ⎧ ⎫⎤ r,orig/alt cos(θgw,2 + α2 ) ⎪ ⎪ ⎨ π − θgw,2 + θgw,3 + α2 + α3 − cos θgw,2 ⎬⎥ ⎢ sin(α2 ) ⎢ 2 σgw ⎥ ⎢ r ⎪ ⎪⎥ ⎢ ⎥ ⎩ ⎭ r,orig/alt cos(θgw,3 + α3 ) ⎢ ⎥ − cos θgw,3 ⎢ ⎥ sin(α ) 3 ⎢ ⎥ ⎫ ⎢ ⎥ ⎧ ⎥ cos θow,3 + α3 ⎪ cos θow,2 + α2 + σgw ⎢ ⎪ ⎢ ⎥ ⎨ ⎬ 2 r + ow σow ⎢ ⎥ r,alt sin (α2 ) sin (α3 ) cos θgw ⎢+ ⎥ ⎢ ⎥ ⎪ ⎪ row ⎩ ⎭ ⎢ ⎥ orig orig − L s,2 − L s,3 ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ σow orig r,orig/alt orig r,orig/alt L s,2 cos θgw,2 + L s,3 cos θgw,3 + row ⎡
σow
(65)
⎡ ⎤ 2 Aw,2 + Aw,3 − At ⎥ ⎢ ⎧ cos θow,2 + α2 ⎫⎥ σow σgo ⎢ ⎪ ⎪ ⎢ ⎨ θow,2 + θow,3 + α2 + α3 − π + cos θow,2 ⎬⎥ C3 = ⎢ ⎥ sin (α2 ) 2 ⎥ row ⎢ + row ⎣ ⎦ ⎪ ⎪ cos θow,3 + α3 ⎩ ⎭ + cos θow,3 sin (α3 ) + 2 σow σgo row θow,2 + θow,3 + α2 + α3 − π cos θow,3 + α3 cos θow,2 + α2 2 r + σgo cos θgo + L t − 2 row sin (α2 ) sin (α3 ) orig r,orig/alt orig r,orig/alt + σgw σgo L s,2 cos θgw,2 + L s,3 cos θgw,3 cos θow,3 + α3 cos θow,2 + α2 orig orig r,alt + σgw σgo cos θgw + − L s,2 − L s,3 2 row sin (α2 ) sin (α3 ) ⎧ ⎫⎤ ⎪ ⎪ ⎨ θow,2 + θow,3 + α2 + α3 − π ⎬ ⎥ ⎢ r2 cos θow,2 + α2 cos θow,3 + α3 ⎪⎥ 2 ⎢ ow ⎪ C4 = σgo ⎢ ⎩ cos θow,2 ⎭⎥ + cos θow,3 ⎣ ⎦ sin (α2 ) sin (α3 ) ⎡
(66)
(67)
+ Aw,2 + Aw,3 − At
– PL displacement with the formation of two spreading oil layers in corners 1 and 2: 3 2 C1 + r go C2 + r go C3 + C4 = 0 r go
(68)
where C1 –C4 are calculated from: C1 =
cos(θ r + α1 ) cos(θ r + α2 ) σow σgo go go r r − cos θgo π − α1 − α2 − 2 θgo + row sin(α1 ) sin(α2 ) (69)
123
Pore-Scale Network Modeling of Three-Phase Flow Based on…
⎤
Aw,3 π + 2 row θow,3 + α3 − ⎢ ⎥ row 2 ⎢ ⎥ ⎢ ⎥ cos θ + α ⎥ ow,3 3 σow ⎢ r ⎢ + σgo cos θgo L t − 2 row ⎥ C2 = ⎥ sin (α3 ) row ⎢ ⎢ ⎥ ⎥ ⎢ cos θ + α ⎣ ⎦ ow,3 3 orig r,orig/alt orig r,alt + σgw L s,3 cos θgw,3 + cos θgw 2 row − L s,3 sin (α3 ) r + α ) cos(θ r + α1 ) cos(θgo 2 go 2 r r π − α1 − α2 − 2 θgo + σgo − cos θgo + sin(α1 ) sin(α2 ) r π r,orig/alt cos(θgw,3 + α3 ) 2 + 2 σgw (70) − θgw,3 − α3 − cos θgw,3 2 sin(α3 ) ) ( cos θow,3 + α3 σow σgo π 2 2 Aw,2 − At + row C3 = θow,3 + α3 − + cos θow,3 row 2 sin (α3 ) cos θow,3 + α3 π 2 r + 2 σow σgo row θow,3 + α3 − + σgo cos θgo L t − 2 sin (α3 ) cos θow,3 + α3 orig r,orig/alt orig r,alt + cos θgw 2 row + σgw σgo L s,3 cos θgw,3 − L s,3 sin (α3 ) ⎡
σow
( 2 C4 =σgo
) cos θow,3 + α3 π 2 + Aw,3 − At θow,3 + α3 − + cos θow,3 row 2 sin (α3 )
(71) (72)
– PL displacement with the formation of three spreading oil layers: see Eq. (30) in Sect. 4. – PL displacement without formation of any spreading oil layers: see Eqs. (26)–(29) in Sect. 4.
Three-Phase Layer Collapse Events – PL displacement with the single collapse of a spreading oil layer in corner 3: see Eqs. (31)–(35) in Sect. 4. – PL displacement with the simultaneous collapse of two spreading oil layers in corners 2 and 3: 3 r gw
2 C σow 1 2 2 σow σgw C 2 + r gw + r gw C3 + σgw C4 = 0 2 row row
(73)
where C1 –C4 are expressed by: cos(θgw,2 + α2 ) C1 = θgw,2 + θgw,3 + α2 + α3 − π + cos θgw,2 sin(α2 ) cos(θgw,3 + α3 ) + cos θgw,3 sin(α ) 3 cos θgw,2 + α2 r,orig/alt cos θgw,2 − cos θgw,2 C2 = sin (α2 ) cos θgw,3 + α3 r,orig/alt cos θgw,3 − cos θgw,3 + sin (α3 ) Aw,2 + Aw,3 2 C3 = σow 2 π − θow,2 − θow,3 − α2 − α3 − 2 row
(74)
(75)
123
A. Zolfaghari, M. Piri
cos θ cos θow,3 + α3 ow,2 + α2 + sin (α2 ) sin (α3 ) ⎤ ⎡ orig r,orig/alt orig r,orig/alt L s,2 cos θgw,2 + L s,3 cos θgw,3 ⎥ ⎢ ⎫⎥ ⎢ ⎧ ⎥ σow σgw ⎢ cos θ + α + α cos θ ow,2 2 ow,3 3 ⎪ ⎢ ⎨ 2 row ⎬⎥ − ⎪ + ⎥ r,alt row ⎢ sin (α2 ) sin (α3 ) ⎥ ⎢ + cos θgw ⎦ ⎣ ⎪ ⎪ ⎩ ⎭ orig orig − L s,2 − L s,3 ⎞⎫ ⎧ ⎛ r +α r +α ⎬ ⎨ cos θgo cos θgo 2 3 2 r r ⎠ + σgo 2 θgo + α2 + α3 − π + cos θgo ⎝ + ⎭ ⎩ sin (α2 ) sin (α3 ) ⎧ ⎫ cos θgw,2 + α2 r,orig/alt ⎪ 2 cos θgw,2 − cos θgw,2 ⎪ ⎨ π − θgw,2 − α2 − ⎬ sin (α2 ) 2 + σgw ⎪ ⎪ cos θgw,3 + α3 ⎩ ⎭ r,orig/alt − cos θgw,3 2 cos θgw,3 − θgw,3 − α3 − sin (α3 )
r + 2 σow σgo cos θgo
(76)
Aw,2 + Aw,3 + 2 row θow,2 + θow,3 + α2 + α3 − π row cos θ cos θow,3 + α3 ow,2 + α2 r + − 2 σgo row cos θgo sin (α2 ) sin (α3 ) ⎤ ⎡ orig r,orig/alt orig r,orig/alt L s,2 cos θgw,2 + L s,3 cos θgw,3 ⎥ ⎢ ⎫⎥ ⎢ ⎧ ⎥ ⎢ cos θ + α + α cos θ ow,2 2 ow,3 3 ⎪ ⎨ 2 row ⎬⎥ + σgw ⎢ ⎪ + ⎥ (77) ⎢ r,alt sin (α2 ) sin (α3 ) ⎥ ⎢ + cos θgw ⎣ ⎪ ⎪ ⎩ ⎭⎦ orig orig − L s,2 − L s,3
C4 = σow
– PL displacement with the simultaneous collapse of three spreading oil layers (configuration changes D → F in Fig. 1):
3 r gw
2 C σow 1 2 2 σow σgw C 2 + r gw + r gw C3 + σgw C4 = 0 2 row row
(78)
where C1 –C4 are given by:
cos(θgw,1 + α1 ) C1 = θgw,1 + θgw,2 + θgw,3 − π + cos θgw,1 sin(α1 ) cos(θgw,2 + α2 ) cos(θgw,3 + α3 ) + cos θgw,2 + cos θgw,3 sin(α2 ) sin(α3 )
123
(79)
Pore-Scale Network Modeling of Three-Phase Flow Based on…
cos θgw,1 + α1 r,orig/alt − cos θgw,1 cos θgw,1 C2 = sin (α1 ) cos θgw,2 + α2 r,orig/alt cos θgw,2 − cos θgw,2 + sin (α2 ) cos θgw,3 + α3 r,orig/alt − cos θgw,3 cos θgw,3 + sin (α3 )
(80)
Aw,1 + Aw,2 + Aw,3 2 C3 = σow 2 π − θow,1 − θow,2 − θow,3 − 2 row cos θow,1 + α1 cos θow,2 + α2 cos θow,3 + α3 r + 2 σow σgo cos θgo + + sin (α1 ) sin (α2 ) sin (α3 ) ⎤ ⎡ orig r,orig/alt orig r,orig/alt orig r,orig/alt + L s,2 cos θgw,2 + L s,3 cos θgw,3 cos θgw,1 L ⎢ s,1 ⎥ ⎞ ⎛ ⎢ ⎥ ⎧ ⎢ ⎥ cos θow,2 + α2 ⎫ cos θow,1 + α1 ⎢ ⎥ ⎪ ⎪ + ⎟ ⎜ ⎪ ⎪ ⎢ ⎥ σow σgw ⎢ ⎪ ⎪ sin sin ) ) (α (α ⎟ ⎜ 1 2 ⎪ ⎪ ⎥ ⎨ ⎬ − 2 row ⎜ ⎟ ⎢ ⎥ row ⎢ + cos θ r,alt ⎠ ⎝ ⎥ cos θow,3 + α3 gw ⎢ ⎥ + ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ sin ) (α 3 ⎪ ⎪ ⎣ ⎦ ⎩ ⎭ orig orig orig − L s,1 − L s,2 − L s,3 ⎞ ⎛ r +α r +α r +α cos θgo cos θgo cos θgo 1 2 3 2 r r ⎠ ⎝ 3 θgo − π + cos θgo + + + σgo sin (α1 ) sin (α2 ) sin (α3 ) cos θgw,1 + α1 r,orig/alt ⎧π − θ 2 cos θgw,1 − cos θgw,1 ⎫ gw,1 − ⎪ ⎪ sin (α1 ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎬ cos θgw,2 + α2 2 r,orig/alt + σgw − θgw,2 − (81) 2 cos θgw,2 − cos θgw,2 ⎪ ⎪ sin (α2 ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ cos θgw,3 + α3 r,orig/alt 2 cos θgw,3 − cos θgw,3 − θgw,3 − sin (α3 )
Aw,1 + Aw,2 + Aw,3 C4 = σow + 2 row θow,1 + θow,2 + θow,3 − π row cos θow,1 + α1 cos θow,2 + α2 cos θow,3 + α3 r − 2 σgo row cos θgo + + sin (α1 ) sin (α2 ) sin (α3 ) ⎤ ⎡ orig r,orig/alt orig r,orig/alt orig r,orig/alt + L s,2 cos θgw,2 + L s,3 cos θgw,3 L s,1 cos θgw,1 ⎢ ⎥ ⎞ ⎛ ⎢ ⎥ ⎧ ⎢ ⎥ cos θow,2 + α2 ⎫ cos θow,1 + α1 ⎢ ⎥ ⎪ ⎪ + ⎜ ⎟⎪ ⎪ ⎢ ⎥ ⎪ ⎪ sin sin ) ) (α (α ⎜ ⎟ 1 2 ⎪ ⎪ ⎢ ⎥ (82) + σgw ⎢ ⎨ 2 row ⎜ ⎟⎬ ⎥ r,alt ⎝ ⎠ ⎢ + cos θgw ⎥ cos θow,3 + α3 ⎢ ⎥ + ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ sin ) (α 3 ⎪ ⎪ ⎣ ⎦ ⎩ ⎭ orig orig orig − L s,1 − L s,2 − L s,3
References Al-Dhahli, A., Van Dijke, M.I., Geiger, S.: Accurate modelling of pore-scale films and layers for three-phase flow processes in clastic and carbonate rocks with arbitrary wettability. Transp. Porous Media 98(2), 259–286 (2013)
123
A. Zolfaghari, M. Piri Alizadeh, A.H., Khishvand, M., Ioannidis, M.A., Piri, M.: Multi-scale experimental study of carbonated water injection: an effective process for mobilization and recovery of trapped oil. Fuel 132, 219–235 (2014) Alizadeh, A.H., Piri, M.: Three-phase flow in porous media: a review of experimental studies on relative permeability. Rev. Geophys. 52(3), 468–521 (2014) Alizadeh, A.H., Piri, M.: The effect of saturation history on three-phase relative permeability: an experimental study. Water Resour. Res. 50(2), 1636–1664 (2014) Baker, L.E.: Three-phase relative permeability correlations, Paper SPE 17369. In: SPE Enhanced Oil Recovery Symposium. Society of Petroleum Engineers (1988) Bartell, F.E., Osterhof, H.J.: Determination of the wettability of a solid by a liquid. Ind. Eng. Chem. 19(11), 1277–1280 (1927) Blunt, M.J.: Pore level modeling of the effects of wettability. SPE J. 2(04), 494–510 (1997) Blunt, M.J.: Physically-based network modeling of multiphase flow in intermediate-wet porous media. J. Petrol. Sci. Eng. 20(3), 117–125 (1998) Blunt, M.J.: Flow in porous media pore-network models and multiphase flow. Curr. Opin. Colloid Interface Sci. 6(3), 197–207 (2001) Blunt, M.J., Jackson, M.D., Piri, M., Valvatne, P.H.: Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow. Adv. Water Resour. 25(8), 1069–1089 (2002) Blunt, M.J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P., Paluszny, A., Pentland, C.: Pore-scale imaging and modelling. Adv. Water Resour. 51, 197–216 (2013) Blunt, M.J., Scher, H.: Pore-level modeling of wetting. Phys. Rev. E 52(6), 6387–6403 (1995) Celia, M.A., Reeves, P.C., Ferrand, L.A.: Recent advances in pore scale models for multiphase flow in porous media. Rev. Geophys. 33(S2), 1049–1057 (1995) Chen, S., Doolen, G.D.: Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech. 30(1), 329–364 (1998) Delshad, M., Pope, G.A.: Comparison of the three-phase oil relative permeability models. Transp. Porous Media 4(1), 59–83 (1989) DiCarlo, D.A., Sahni, A., Blunt, M.J.: Three-phase relative permeability of water-wet, oil-wet, and mixed-wet sandpacks. SPE J. 5(01), 82–91 (2000) Dixit, A.B., McDougall, S.R., Sorbie, K.S., Buckley, J.S.: Pore scale modelling of wettability effects and their influence on oil recovery. SPE Reserv. Eval. Eng. 2(1), 25–36 (1999) Dong, M., Dullien, F.A., Chatzis, I.: Imbibition of oil in film form over water present in edges of capillaries with an angular cross section. J. Colloid Interface Sci 172(1), 21–36 (1995) Dria, D.E., Pope, G.A., Sepehrnoori, K.: Three-phase gas/oil/brine relative permeabilities measured under CO2 flooding conditions. SPE Reserv. Eng. 8(02), 143–150 (1993) Ehrlich, R., Davies, D.K.: Image analysis of pore geometry: relationship to reservoir engineering and modeling, Paper SPE 19054. In: SPE Gas Technology Symposium. Society of Petroleum Engineers (1989) Fatt, I.: The network model of porous media I. Capillary pressure characteristics. Trans Soc. Min. Eng. AIME 207, 144–159 (1956) Fatt, I.: The network model of porous media II. Dynamic properties of a single size tube network. Trans. AIME 207, 160–163 (1956) Fatt, I.: The network model of porous media III. Dynamic properties of networks with tube radius distribution. Trans. AIME 207, 164–181 (1956) Fenwick, D.H., Blunt, M.J.: Network modeling of three-phase flow in porous media. SPE J. 3(01), 86–96 (1998) Fenwick, D.H., Blunt, M.J.: Three-dimensional modeling of three phase imbibition and drainage. Adv. Water Resour. 21(2), 121–143 (1998) Ferreol, B., Rothman, D.H.: Lattice-Boltzmann simulations of flow through Fontainebleau sandstone. Transp. Porous Media 20(1), 3–20 (1995) Firincioglu, T., Blunt, M.J., Zhou, D.: Three-phase flow and wetability effects in triangular capillaries. Colloids Surf. A 155(2), 259–276 (1999) Goode, P.A., Ramakrishnan, T.S.: Momentum transfer across fluid–fluid interfaces in porous media: a network model. AIChE J. 39(7), 1124–1134 (1993) Grader, A.S., O’Meara Jr, D.J.: Dynamic displacement measurements of three-phase relative permeabilities using three immiscible liquids, Paper SPE 18293. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers (1988) Grunau, D., Chen, S., Eggert, K.: A lattice Boltzmann model for multiphase fluid flows. Phys. Fluids A 5(10), 2557–2562 (1993) Gunstensen, A.K., Rothman, D.H., Zaleski, S., Zanetti, G.: Lattice Boltzmann model of immiscible fluids. Phys. Rev. A 43(8), 4320–4327 (1991)
123
Pore-Scale Network Modeling of Three-Phase Flow Based on… Hayden, N.J., Voice, T.C.: Microscopic observation of a NAPL in a three-fluid-phase soil system. J. Contam. Hydrol. 12(3), 217–226 (1993) Heiba, A.A., Davis, H.T., Scriven, L.E.: Statistical Network Theory of Three-Phase Relative Permeabilities, Paper SPE 12690. In: SPE Enhanced Oil Recovery Symposium. Society of Petroleum Engineers (1984) Honarpour, M.M., Koederitz, F., Herbert, A.: Relative Permeability of Petroleum Reservoirs. CRC Press, Florida (1986) Hu, X.Y., Adams, N.A.: A multi-phase SPH method for macroscopic and mesoscopic flows. J. Comput. Phys. 213(2), 844–861 (2006) Hui, M.H., Blunt, M.J.: Effects of wettability on three-phase flow in porous media. J. Phys. Chem. B 104(16), 3833–3845 (2000) Keller, A.A., Blunt, M.J., Roberts, A.P.V.: Micromodel observation of the role of oil layers in three-phase flow. Transp. Porous Media 26(3), 277–297 (1997) Khishvand, M., Alizadeh, A.H., Piri, M.: In-situ characterization of wettability and pore-scale displacements during two-and three-phase flow in natural porous media. Advances in Water Resources 97, 279–298 (2016) Kovscek, A.R., Wong, H., Radke, C.J.: A pore-level scenario for the development of mixed wettability in oil reservoirs. AIChE J. 39(6), 1072–1085 (1993) Lago, M., Araujo, M.: Threshold pressure in capillaries with polygonal cross section. J. Colloid Interface Sci. 243(1), 219–226 (2001) Lago, M., Araujo, M.: Threshold capillary pressure in capillaries with curved sides. Phys. A 319, 175–187 (2003) Laroche, C., Vizika, O., Kalaydjian, F.: Network modeling to predict the effect of wettability heterogeneities on multiphase flow, Paper SPE 56674. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers (1999) Laroche, C., Vizika, O., Kalaydjian, F.: Network modeling as a tool to predict three-phase gas injection in heterogeneous wettability porous media. J. Petrol. Sci. Eng. 24(2), 155–168 (1999) Larsen, J.K., Bech, N., Winter, A.: Three-phase immiscible WAG injection: Micromodel experiments and network models, Paper SPE 59324. In: SPE/DOE Improved oil recovery symposium. Society of Petroleum Engineers.(2000) Lenormand, R., Zarcone, C.: Role of roughness and edges during imbibition in square capillaries, Paper SPE 13264. In: the 59th SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers (1984) Lenormand, R., Zarcone, C., Sarr, A.: Mechanisms of the displacement of one fluid by another in a network of capillary ducts. J. Fluid Mech. 135, 337–353 (1983) Lerdahl, T.R., Øren, P.E., Bakke, S.: A predictive network model for three-phase flow in porous media, Paper SPE 59311. In: SPE/DOE Improved Oil Recovery Symposium. Society of Petroleum Engineers (2000) Liu, G.R., Liu, M.B.: Smoothed Particle Hydrodynamics: A Meshfree Particle Method. World Scientific Publishing Co Pte Ltd, Singapore (2003) Lopez, X., Valvatne, P.H., Blunt, M.J.: Predictive network modeling of single-phase non-newtonian flow in porous media. J. Colloid Interface Sci. 264(1), 256–265 (2003) Maloney, D., Brinkmeyer, A.D.: Three-Phase Permeabilities and Other Characteristics of 260-mD Fired Berea. IIT Research Institute, National Institute for Petroleum and Energy Research (1992) Maloney, D.R., Doggett, K., Tomutsa, L.: Three-Phase Relative Permeability Characteristics of a dolomite rock plug. Report NIPER, 329 (1997) Mani, V., Mohanty, K.K.: Effect of the spreading coefficient on three-phase flow in porous media. J. Colloid Interface Sci. 187(1), 45–56 (1997) Mani, V., Mohanty, K.K.: Pore-level network modeling of three-phase capillary pressure and relative permeability curves. SPE J. 3(03), 238–248 (1998) Man, H.N., Jing, X.D.: Network modelling of wettability and pore geometry effects on electrical resistivity and capillary pressure. J. Petrol. Sci. Eng. 24(2), 255–267 (1999) Mason, G., Morrow, N.R.: Capillary behavior of a perfectly wetting liquid in irregular triangular tubes. J. Colloid Interface Sci. 141(1), 262–274 (1991) Mayer, R.P., Stowe, R.A.: Mercury porosimetry breakthrough pressure for penetration between packed spheres. J. Colloid Sci. 20(8), 893–911 (1965) McDougall, S.R., Sorbie, K.S.: The impact of wettability on waterflooding: pore-scale simulation. SPE Reserv. Eng. 10(03), 208–213 (1995) Monaghan, J.J.: An introduction to SPH. Comput. Phys. Commun. 48(1), 89–96 (1988) Morris, J.P., Fox, P.J., Zhu, Y.: Modeling low Reynolds number incompressible flows using SPH. J. Comput. Phys. 136(1), 214–226 (1997)
123
A. Zolfaghari, M. Piri Oak, M.J.: Three-phase relative permeability of water-wet Berea, Paper SPE 20183. In: SPE/DOE Enhanced Oil Recovery Symposium. Society of Petroleum Engineers (1990) Oak, M.J., Baker, L.E., Thomas, D.C.: Three-phase relative permeability of Berea sandstone. J. Petrol. Technol. 42(08), 1–054 (1990) Oliveira, L.I., Demond, A.H.: Estimation of primary drainage three-phase relative permeability for organic liquid transport in the vadose zone. J. Contam. Hydrol. 66(3), 261–285 (2003) Øren, P.E., Billiotte, J., Pinczewski, W.V.: Pore-scale network modelling of waterflood residual oil recovery by immiscible gas flooding, Paper SPE 27814. In: SPE/DOE Improved Oil Recovery Symposium. Society of Petroleum Engineers (1994) Øren, P. E., Pinczewski, W. V.: The effect of wettability and spreading coefficients on the recovery of waterflood residual oil by miscible gasflooding, Paper SPE 24881. In: the 67th SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers (1992) Øren, P.E., Bakke, S., Arntzen, O.J.: Extending predictive capabilities to network models. SPE J. 3(04), 324– 336 (1998) Øren, P.E., Pinczewski, W.V.: Fluid distribution and pore-scale displacement mechanisms in drainage dominated three-phase flow. Transp. Porous Media 20(1–2), 105–133 (1995) Ovaysi, S., Piri, M.: Direct pore-level modeling of incompressible fluid flow in porous media. J. Comput. Phys. 229(19), 7456–7476 (2010) Ovaysi, S., Piri, M.: Pore-scale modeling of dispersion in disordered porous media. J. Contam. Hydrol. 124(1), 68–81 (2011) Paterson, L., Lee, J.Y., Pinczewski, W.V.: Three-phase relative permeability in heterogeneous formations, Paper SPE 38882. In: SPE annual technical conference and exhibition. Society of Petroleum Engineers (1997) Patzek, T.W.: Verification of a complete pore network simulator of drainage and imbibition. SPE J. 6(2), 144–156 (2001) Pereira, G.G., Pinczewski, W.V., Chan, D.Y.C., Paterson, L., Øren, P.E.: Pore-scale network model for drainagedominated three-phase flow in porous media. Transp. Porous media 24(2), 167–201 (1996) Piri, M., Blunt, M.J.: Three-phase threshold capillary pressures in noncircular capillary tubes with different wettabilities including contact angle hysteresis. Phys. Rev. E 70(6), 061603 (2004) Piri, M., Blunt, M.J.: Three-dimensional mixed-wet random pore-scale network modeling of two-and threephase flow in porous media. I. Model description. Phys. Rev. E 71(2), 026301 (2005) Piri, M., Blunt, M.J.: Three-dimensional mixed-wet random pore-scale network modeling of two-and threephase flow in porous media. II. Results. Phys. Rev. E 71(2), 026302 (2005) Princen, H.M.: Capillary phenomena in assemblies of parallel cylinders: I. Capillary rise between two cylinders. J. Colloid Interface Sci. 30(1), 69–75 (1969) Princen, H.M.: Capillary phenomena in assemblies of parallel cylinders: II. Capillary rise in systems with more than two cylinders. J. Colloid Interface Sci. 30(3), 359–371 (1969) Princen, H.M.: Capillary phenomena in assemblies of parallel cylinders: III. Liquid columns between horizontal parallel cylinders. J. Colloid Interface Sci. 34(2), 171–184 (1970) Rowlinson, J.S., Widom, B.: Molecular Theory of Capillarity. Dover Inc., New York (1982) Sahni, A., Burger, J., Blunt, M.J.: Measurement of three phase relative permeability during gravity drainage using CT scanning, Paper SPE 39655. In: SPE/DOE Improved Oil Recovery Symposium. Society of Petroleum Engineers (1998) Saraf, D.N., Batycky, J.P., Jackson, C.H., Fisher, D.B.: An experimental investigation of three-phase flow of water–oil–gas mixtures through water-wet sandstones, Paper SPE 10761. In: SPE California Regional Meeting. Society of Petroleum Engineers (1982) Soll, W.E., Celia, M.A.: A modified percolation approach to simulating three-fluid capillary pressure-saturation relationships. Adv. Water Resour. 16(2), 107–126 (1993) Stone, H.L.: Probability model for estimating three-phase relative permeability. J. Petrol. Technol. 22(02), 214–218 (1970) Stone, H.L.: Estimation of three-phase relative permeability and residual oil data. J. Can. Pet. Technol. 12(4), 53–61 (1973) Suicmez, V.S., Piri, M., Blunt, M.J.: Pore-scale modeling of three-phase WAG injection: prediction of relative permeabilities and trapping for different displacement cycles, Paper SPE 95594. In: SPE/DOE Symposium on Improved Oil Recovery. Society of Petroleum Engineers (2006) Suicmez, V.S., Piri, M., Blunt, M.J.: Pore-scale simulation of water alternate gas injection. Transp. Porous Media 66(3), 259–286 (2007) Suicmez, V.S., Piri, M., Blunt, M.J.: Effects of wettability and pore-level displacement on hydrocarbon trapping. Adv. Water Resour. 31(3), 503–512 (2008)
123
Pore-Scale Network Modeling of Three-Phase Flow Based on… Svirsky, D. S., Van Dijke, M. I., Sorbie, K. S.: Prediction of three-phase relative permeabilities using a porescale network model anchored to two-phase data, Paper SPE 89992. In: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers (2007) Tartakovsky, A.M., Meakin, P.: A smoothed particle hydrodynamics model for miscible flow in threedimensional fractures and the two-dimensional Rayleigh Taylor instability. J. Comput. Phys. 207(2), 610–624 (2005) Tartakovsky, A.M., Meakin, P.: Pore scale modeling of immiscible and miscible fluid flows using smoothed particle hydrodynamics. Adv. Water Resour. 29(10), 1464–1478 (2006) Tsakiroglou, C.D., Payatakes, A.C.: Characterization of the pore structure of reservoir rocks with the aid of serial sectioning analysis, mercury porosimetry and network simulation. Adv. Water Resour. 23(7), 773–789 (2000) Valvatne, P.H., Piri, M., Lopez, X., Blunt, M.J.: Predictive pore-scale modeling of single and multiphase flow. Transp. Porous Media 58, 23–41 (2005) Valvatne, P.H., Blunt, M.J.: Predictive pore-scale modeling of two-phase flow in mixed wet media. Water Resour. Res. 40(7), W07406 (2004) Van Dijke, M. I. J., Piri, M., Helland, J. O., Sorbie, K. S., Blunt, M. J., Skjæveland, S. M.: Criteria for three-fluid configurations including layers in a pore with nonuniform wettability. Water Resour. Res. 43(12) (2007) Van Dijke, M.I.J., Sorbie, K.S., McDougall, S.R.: A process-based approach for three-phase capillary pressure and relative permeability relationships in mixed-wet systems, Paper SPE 59310. In: SPE/DOE Improved Oil Recovery Symposium. Society of Petroleum Engineers. (2000) Van Dijke, M.I.J., Sorbie, K.S., McDougall, S.R.: Saturation-dependencies of three-phase relative permeabilities in mixed-wet and fractionally wet systems. Adv. Water Resour. 24(3), 365–384 (2001) Van Dijke, M.I.J., McDougall, S.R., Sorbie, K.S.: Three-phase capillary pressure and relative permeability relationships in mixed-wet systems. Transp. Porous Media 44(1), 1–32 (2001) Van Dijke, M.I.J., Sorbie, K.S., Sohrabi, M., Tehrani, D., Danesh, A.: Three-phase flow in WAG processes in mixed-wet porous media: pore-scale network simulations and comparison with micromodel experiments. SPE J. 9(2), 57–66 (2004) Van Dijke, M.I.J., Lago, M., Sorbie, K.S., Araujo, M.: Free energy balance for three fluid phases in a capillary of arbitrarily shaped cross-section: capillary entry pressures and layers of the intermediate-wetting phase. J. Colloid Interface Sci. 277(1), 184–201 (2004) Van Dijke, M.I.J., Sorbie, K.S.: Pore-scale network model for three-phase flow in mixed-wet porous media. Phys. Rev. E 66(4), 046302 (2002) Van Dijke, M.I.J., Sorbie, K.S.: Three-phase capillary entry conditions in pores of noncircular cross-section. J. Colloid Interface Sci. 260(2), 385–397 (2003) Van Dijke, M.I.J., Sorbie, K.S.: Existence of fluid layers in the corners of a capillary with non-uniform wettability. J. Colloid Interface Sci. 293(2), 455–463 (2006) Van Dijke, M.I.J., Sorbie, K.S.: Cusp at the three-fluid contact line in a cylindrical pore. J. Colloid Interface Sci. 297(2), 762–771 (2006) Van Kats, F.M., Egberts, P.J.P., Van Kruijsdijk, C.P.J.W.: Three-phase effective contact angle in a model pore. Transp. Porous Media 43(2), 225–238 (2001) Van Kats, F.M., Egberts, P.J.P.: Simulation of three-phase displacement mechanisms using a 2D latticeBoltzmann model. Transp. Porous Media 37(1), 55–68 (1999) Xu, B., Kamath, J., Yortsos, Y.C., Lee, S.H.: Use of pore network models to simulate laboratory corefloods in a heterogeneous carbonate sample. SPE J. 4(03), 179–186 (1999) Zhou, D., Blunt, M.J.: Effect of spreading coefficient on the distribution of light non-aqueous phase liquid in the subsurface. J. Contam. Hydrol. 25(1), 1–19 (1997) Zhu, Y., Fox, P.J., Morris, J.P.: A pore-scale numerical model for flow through porous media. Int. J. Numer. Anal. Methods Geomech. 23(9), 881–904 (1999) Zolfaghari, A., Piri, M., Pore-scale network modeling of three-phase flow based on thermodynamically consistent threshold capillary pressures. II. Results. Transp. Porous Media. doi:10.1007/11242-016-0815-7 (2016) Zolfaghari, A.: Pore-Scale network modeling of two- and three-phase flow based on thermodynamically consistent threshold capillary pressures. Ph.D. dissertation, University of Wyoming (2014)
123
