Pore-Scale Modelling of Three-Phase Capillary Pressure Curves Directly in Uniformly Wet Rock Images

State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan, China School of Engineering, University of Aberdeen, King’s College, Aberdeen, UK AB24 3UE Petroleum Engineering Department, University of Houston, 5000 Gulf Freeway, Technology Bridge, Building 9B, Room 158, Houston, Texas 77204-0945, USA International Research Institute of Stavanger, P.O. Box 8046, 4068 Stavanger, Norway State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing, China


Introduction
Pore-scale modelling has been used extensively to quantify properties related to two-phase flow in porous media, namely, capillary pressure curves and relative permeability curves. Pore-scale modelling includes pore network models [1] and direct simulation techniques proposed to simulate multiphase in segmented 3D rock images, such as the Lattice Boltzmann methods [2,3], morphology-based methods [2,4], and level set methods [5,6]. Capillary bundle models, despite their simplistic nature, do account for some of the important capillary phenomena that occur in porous media; for example, described mixed wettability uses an angular tube cross-section [7], formulates physical-based capillary pressure correlations [8], and provides analysis of spontaneous imbibition [9].
Three-phase flow in porous media is relevant to describe oil displacement processes, such as Water-Alternate-Gas (WAG) injection and CO 2 injection to improve oil recovery from mature oil reservoirs. Compared with the two-phase flow, the three-phase flow in porous media is much more challenging and complicated as it contains three fluid-fluid interfaces, gas-oil, oil-water, and gas-water, and thus results in more displacement scenarios [10][11][12][13]. The fluid invasion, fluid configurations, and flow functions strongly depend on the interfical properties of the fluids and porous fluids and solid surfaces of the porous medium [14][15][16][17]. Quasistatic displacement is a reasonable approximation of the threephase flow at the pore scale, and pore-network models can be used to simulate the capillary-dominated threephase flow in porous media [10,12,15,[17][18][19]. The pore network modeling requires first to extract the pore network from 3D images of porous media with the pore space being represented by an interconnected system of pores and throats having an idealized geometry and constant cross-section [20]. Pore network modeling of the three-phase flow also adopts an analytical entry capillary pressure together with invasion percolation algorithms to describe displacement in the predetermined pore network [12,21]. The advantages of pore network models for the three-phase flow compared with the most recently developed direct methods, such as the level set method [11], pore-morphology method [22], and Lattice Boltzmann method [23,24], are their efficiency and lack of dependence on grid resolution effects used to represent thin oil layers in the idealized geometry.
Due to the simplified pore space representation using interconnected tubes with analytical geometries, the pore network model fails to capture the complexity of the pore space geometry. Recently, a generalized pore network model was proposed aimed at providing more realistic pore geometries [25,26]. In this model, the pore throat space is approximated with semianalytical polygons; thus, the local roughness of the pore throat and its impact on entry capillary pressure and fluid configurations are not well accounted for. In addition to the pore network and other direct pore-scale models, capillary tube bundle models have also proven effective in analyzing systematically the effects of wettability and interfacial tension on three-phase relative permeability and capillary pressure-saturation relationships [8,13,15]. Capillary tubes with analytical pore geometry, such as circular or triangular cross-sections, have been adopted for most of the capillary tube models. In our previous publications, we have developed a semianalytical model to overcome the limitations of a simplified capillary tube model that was constructed with analytical cross-sections; we have applied this model directly on 2D rock images to successfully calculate two-phase capillary pressure curves and fluid configurations in porous media [16,27] and three-phase entry pressure and fluid configurations [13].
In this work, we extend the semianalytical model to simulate three-phase capillary pressure curves and the corresponding three-phase fluid configurations during gas invasion. The model requires a 2D segmented rock image as input and treats the identified pore spaces as crosssections of straight tubes. The capillary entry pressures and corresponding cross-sectional fluid configurations in these pore spaces are determined at any capillary pressure and wetting condition by combining an arc menisci-determining procedure with free energy minimization by means of generalizing the so-called Mayer and Stowe-Princen (MS-P) method [28][29][30]. Using this model, we simulate three-phase capillary pressure curves and fluid configurations in SEM images of Bentheim sandstone for gas invasion at uniformly water-wet conditions.

Semianalytical Model for Computing Fluid Configurations and Capillary Pressure Curves
Pore geometries are considered as straight capillary tubes with constant cross-sections extracted directly from the 2D segmented rock images. The smooth pore/solid boundaries Γ are calculated by a Euclidean path method based on the discrete boundaries that are identified in the binary image [27]. Orientation angles α are then computed for all boundary points b ∈ Γ. The angle α is defined between the boundary tangent at b and a line parallel with the abscissa. In such media, the primary drainage and initial oil-water configurations prior to gas invasion are determined by the semianalytical method [16,27]. The semianalytical model is then implemented to compute three-phase capillary pressure curves and associated fluid configurations in these extracted realistic pore spaces. For each prescribed set of three-phase capillary pressures values, (a) consider each of the three fluid pairs separately and move a circle in opposite directions along the pore boundary to determine the three sets of geometrically permissible interfaces, (b) combine these interfaces systematically and group them into boundary elements of different region types, (c) combine these regions in all geometrically permissible ways to identify the possible candidates for gas and oil configurations, (d) calculate the capillary entry pressure curvature for all these configuration candidates based on the three-phase MS-P equation and compare their entry values to determine the valid displacement, and finally, (e) apply an iterative procedure for the most favorable region combinations invaded by gas in consecutive iteration steps until no more configuration candidates satisfy the criterion (proposed capillary pressure is smaller than the entry pressure of configuration candidate); at this point, the fluid configuration is determined. Steps (a) to (e) are repeated for each prescribed gas-oil capillary pressure. This procedure assumes that for each gas-oil capillary pressure, the oil-water configuration determined at the end of the primary drainage exists initially in the pore space and that no gas is present. However, for gas-oil capillary pressures higher than the capillary entry pressure for gas invasion, the correct three-phase configurations are obtained by applying the iteration procedure in step (e).

Primary Oil Invasion.
Primary drainage at uniform water-wet conditions is simulated first with a contact angle θ ow using the algorithm developed by Frette and Helland [27] and Zhou et al. [16]. In this process, a sequence of oilwater capillary pressures is given, and for each capillary pressure, the pore space is assumed to be occupied initially by water. The method computes geometrically permissible arc menisci (AMs) for each capillary pressure by moving a circle with radius r ow = σ ow /P cow in opposite directions along the pore boundary such that the contact angle θ ow is defined at the front arcs. Circle intersections determine the geometrically allowed AMs. These interfaces are grouped into boundary elements of regions in the pore space that are classified as either BR or LR regions (see the appendix for these definitions). All geometrically permissible region combinations 2 Geofluids are generated, and the capillary entry pressure curvatures for each of these configuration candidates are computed based one the two-phase MS-P equation [16,31]. The valid oilwater configuration at the prescribed capillary pressure is determined iteratively by executing in each step the displacement among all the allowed region combinations which is associated with the most favorable entry pressure being smaller than the considered capillary pressure. An example of geometrically allowed oil-water AMs and determined bulk regions (BRs) and layer regions (LRs) in an individual pore space for r ow = 8 pixel lengths and θ ow = 0°is shown in Figure 1. The associated oil-water configuration before gas invasion for this data is presented in Figure 2.

Gas Invasion
Modeling. The modeling of gas invasion in porous rocks requires a prescribed sequence of gas-oil capillary pressures P cgo . The oil-water capillary pressure P cow is determined from the end of the drainage process, and the gaswater capillary pressure P cgw is obtained from P cgw = P cgo + P cow . The interface radii r ij for each of the three fluid pairs can then be calculated from r ij = σ ij /P cij . We consider uniform water-wet conditions, which imply that gas is the nonwetting phase, oil is the intermediate-wetting phase, and water is the wetting phase [15]. The oil-water contact angle θ ow and the set of interfacial tensions are specified, and the gas-oil and gas-water contact angles, θ go and θ gw , are calculated from the Bartell-Osterhof equation [15]. The next subsections describe the various procedure steps required to determine the threephase fluid configurations in an individual Bentheim sandstone pore space at a given gas-oil capillary pressure.

Possible AMs for Gas Invasion.
For each of the three fluid pairs ij = go, ow, gw, a circle with radius r ij = σ ij /P cij is moved along the pore boundary in opposite directions such that the front arcs form the given contact angle θ ij in the boundary points. The loci of the circle centers constitute two closed curves which are given by where ðx, yÞ are the coordinates of the pore boundary point b.
In Equation (1), v cc ij and v c ij represent the circle centers obtained based on the circle movement along the pore boundary in counterclockwise and clockwise directions, respectively. Geometrically allowed arc menisci for each of the three fluid pairs occur at locations where the circles moving in opposite directions intersect, i.e., v cc ij = v c ij , and the coinciding arcs point toward pore-space constrictions or corners, as shown in Figure 1 for oil-water interfaces in a pore geometry extracted from a Bentheim sandstone. Figure 3 shows the corresponding gas-oil and gas-water interfaces obtained for r ow = 8, r go = 12, and r gw = 9 pixel lengths using θ ow = θ go = θ gw = 0°and interfacial tensions representing a spreading system (C S,o = 0 N/m) that are given by σ go = 0:01 N/m, σ ow = 0:02 N/m, and σ gw = 0:03 N/m.    Figure 3: Example illustrating the gas-oil and gas-water AMs (black curve: pore boundary; red curve: gas-oil; green curve: gas-water) and the extracted regions in a pore geometry extracted from Bentheim sandstone (r go = 12 and r gw = 9 pixel lengths and θ go = θ gw = 0°).

Possible Fluid Configurations for Gas Invasion.
In irregular pore geometries, which might include several wide pore areas connected by narrow constrictions, several fluid configurations are geometrically possible for a given capillary pressure. We determine the possible configuration candidates by combining the interfaces systematically and grouping them into boundary elements of different region types. Individual regions or combinations of regions constitute the configuration candidates which are invaded by gas sequentially in consecutive iteration steps based on which a candidate configuration is identified in each iteration that has the lowest entry curvature in each iteration. In the first iteration, the initial oil-water configuration exists initially in the pore geometry, while in the second iteration, the most favorable gas configuration candidate from the first iteration has invaded the pore, which possibly gives rise to a three-phase configuration with both oil, water, and gas present in the same pore. Therefore, displacements of oil layers surrounded by bulk gas initially and corner water are new displacements that are only considered in iteration steps after the first one. At the end of the second iteration step, gas might have invaded two wider bulk areas that are separated by a layer region (LR) occupied by oil or water in a narrow constriction, and hence, gas invasion into the layer region occurs for the first time in the third or subsequent iteration steps.
The following procedure is common for all iterations: First, the possible gas configuration candidates are determined based on the gas-oil and gas-water systems only, without honoring either the oil-water interfaces or the existing oil-water configuration. Second, the oil-water interfaces are added to the system to determine the possible oil configuration candidates to each of the obtained gas configuration candidate. This step is required for layer displacements (in the second or subsequent iterations), but also to account for situations in which gas invasion occurs with simultaneous invasion of oil that separates the invading gas and the displaced water. Such scenarios are relevant for spreading conditions [30]. The final step is to add oil layers at all gas-oil interfaces to each of the determined gas-and oil-configuration candidates. We make use of the region types defined in the appendix.
(1) Possible Gas Configurations. Possible gas configurations are determined based on the gas-oil and gas-water AMs without honoring the oil-water interfaces. The pore-boundary tracking procedure presented by Frette and Helland [27] is utilized to identify the different regions with boundaries composed of alternate AMs and pore-boundary segments, which subsequently are classified as bulk regions (BRs), layer regions (LRs), and mixed regions (MRs) based on the definitions provided in the appendix. Figure 3 highlights the identified regions in the sandstone pore geometry used as an example.
All allowed combinations of these regions that have not yet become invaded by gas are generated as part of the procedure to determine the valid three-phase fluid configuration. As shown in Figure 3, five BR regions, four LR regions, and 12 MR regions are extracted based on the computed geometrically allowed gas-oil (red) and gas-water (green) interfaces.
Initially (first iteration), none of these regions are occupied by gas. The valid gas configuration is determined among all the geometrically possible combinations of regions that are determined as follows: (i) Individual BR configurations. The number of this type of gas configurations equals the number of bulk gas regions, N BR . In the numerical example shown in Figure 3, N BR = 5, and the five possible gas configurations are given by fBRð1Þg, fBRð2Þg, fBRð3Þg, fBRð4Þg, and fBRð5Þg.
If the number of LR regions N LR ≥ 1, all the geometrically allowed combinations of LR regions are generated and merged with their neighboring BR regions, by also including the cases where the layer and bulk regions are separated by MR regions. There are 2 N LR − 1 different combinations. However, this number may be reduced as we only include the possibilities in which the LR combinations together with their neighboring MRs and BRs form a distinct, extended gas region. As shown in the example of Figure 3, there are ten different connected BR-LR-MR gas regions. These are given as follows: (1) fBRð1Þ, LRð1Þ, BRð2Þg,   In total for all the ten connected BR -LR-MR candidates, there exists 2 0 + 2 7 + 2 9 + 2 8 + 2 7 + 2 9 + 2 8 + 2 10 + 2 9 + 2 3 − 10 = 3327 possible merged BR-MR gas configurations if none of these regions are invaded by gas.
The above procedure determines all possible gas configuration candidates for pore spaces that are not yet invaded by gas. Thus, in the numerical example presented in Figure 3, there are five bulk gas regions and ten connected BR-LR -MR regions, which, together with the corresponding 270 and 3327 merged MR region combinations added, yields a total of 3612 different possible gas entry configurations.
At the present iteration step, if some of the extracted regions are invaded by gas, other displacement scenarios can occur, which must also be examined. These include the displacement of oil layers located between bulk gas and corner water (in the second or subsequent iteration steps) as well as displacement of water or oil layers that are described by LR regions and located in constrictions (in the third or subsequent iteration steps). Such displacement possibilities are evaluated by adding the determined oil-water interfaces, including the initial oil-water configuration, to the pore space.
(i) Displacement of oil layers surrounded by bulk gas and corner water. Gas-oil AMs separate oil layers neighboring the bulk gas configuration. Displacement of such oil layers are considered only if the required oil-water and gas-water AMs can be determined. This is done by first examining whether the gas-oil AM belongs to a mixed region MR based on Definition 6 given in the appendix. If this is the case, then there exists a supporting gas-water AM which also is part of the mixed region boundary. Subsequently, a pore-boundary tracking routine [27] is used on all the gas-oil AMs that belong to both the boundaries of the gas region and a mixed region. This procedure is used to identify supporting oil-water AMs, and it tracks the pore-solid boundary starting with each of the gas-oil AMs and continues along alternate poresolid segments and oil-water AMs until the initial gas-oil AM is reached and closed region boundaries form. Thus, gas invasion into such oil layers is made possible by the existence of the three types of supporting interfaces. We account for situations in which the amount of water in the corners increases when gas-water AMs form. For example, as shown in Figure 4, BRð3Þ constitutes the gas configuration and eight different neighboring oil layers surrounded by gas and corner/neck water. The eight associated MRs, shown in Figure 3, are also identified. In this case, the eight oil layers could be displaced by the invading gas in separate displacements.
(ii) Displacement of layers in constrictions. If a layer region LR is neighboring the invaded gas configuration directly or via a mixed region MR (oil layer), it is considered (together with the mixed region) as a possible candidate for gas invasion if it is included in a combination with a bulk region that has yet to be invaded by gas. An example of such a gas invasion candidate is the region combination fMRð11Þ, BR ð4Þ, MRð9Þ, LRð4Þ, MRð6Þg shown in Figure 3 if region BRð3Þ is already invaded by gas as indicated in Figure 4. However, in subsequent iterations, the layer region may be completely surrounded by invaded bulk gas regions, either directly or via one or more mixed regions (MRs), and in this case, the layer region (together with the mixed regions) are considered to be a candidate for gas invasion alone. For example, if gas has invaded region combination fMRð11Þ, BRð4Þg in the second iteration step as indicated in Figure 5, it is possible that gas will invade the constriction, i.e., combination fMRð9Þ, LRð4Þ, MRð6Þg, in the third iteration step, resulting in an extended merged gas region. Such layer displacements seem most relevant if the LR region is occupied by oil.
(2) Possible Oil Configurations. The next step is to add the geometrically allowed oil-water AMs to the pore space, including the AMs belonging to the existing oil-water configuration before gas invasion, in order to determine possible oil configurations for each of the obtained gas configuration candidates. This step is needed to account for cases where BR (1) LR (1) LR(2) BR (2) BR (4) LR(3) Figure 4: The favorable configuration of gas (yellow area) which invaded the pore space in the first iteration step is shown along with the geometrically allowed oil-water AMs (blue) in a pore space extracted from Bentheim sandstone (r ow = 8, r go = 12, and r gw = 9 pixel lengths and θ ow = θ go = θ gw = 0°).

Geofluids
oil invades the pore together with gas, which is assumed to occur either in the form of spreading oil layers, through oil invasion into combinations of bulk regions (BRs) and layer regions (LRs) that are in contact with the invading gas, or via oil invasion into bulk region combinations that are separated from the invading gas by water located in pore constrictions. Zhou et al. [13] described these features in detail.
Finally, oil layers are added to the bulk gas and oil configuration candidates. The possible oil layers are identified by applying the pore-boundary routine [27] to all gas-oil AMs belonging to both the boundary of mixed regions and the boundary of the considered gas configuration candidate. This procedure tracks the pore-solid boundary, starting with each gas-oil AM in the gas configuration, and continues along alternate pore-solid segments and oil-water AMs until the initial gas-oil AM is reached and closed region boundaries form. Two generic oil-layer configurations are considered in this process. The first type represents oil layers forming between bulk gas and water in the pore corners. The second type represents oil layers forming between bulk gas and water layer regions located in pore constrictions which in turn connect bulk oil regions. In the latter case, while we still account for the oil layer located between bulk gas and the water layer, we also construct an oil layer which occupies both layer regions and merges with the bulk oil configuration, leading to an expanded oil layer adjoining the bulk gas. All oil layers generated this way are added to the bulk oil and gas configurations. Obviously, gas-water AMs that are part of the gas configuration boundary do not support oil-layer formation.
Summarizing the procedure, the possible three-phase configurations are generated by (i) determining the different gas configuration candidates (including layer displacements in the second and subsequent iteration steps), (ii) determining the possible bulk oil configurations for each gas configuration, and (iii) determining the possible oil layer configurations for each bulk oil-and gas-configuration combination. The gas-oil capillary entry pressure curvatures for all determined fluid configuration candidates are calculated by the three-phase MS-P equation given by Equation (2).

Three-Phase Entry
Pressure. Capillary entry pressures for three-phase displacements, including oil-layer displacements, in tubes with idealized and angular cross-sections have been calculated using the MS-P method, The MS-P method was initially proposed to calculate the two-phase entry pressures in straight tubes with constant cross-sections of idealized shapes [1,28,29] and, later, in tubes with arbitrary, yet relatively convex, polygonal crosssections by making use of the relation between the entry fluid configuration and the medial axis of the pore space [32]. Previously, we used to calculate three-phase capillary entry pressures and associated fluid configurations for gas invasion under uniform water-wet conditions [13]. In this work, we further developed the semianalytical method to simulate complete three-phase gas invasion processes.

Computational Procedure.
Fluid configuration at a given gas-oil interface radius r go is determined by the aforementioned iterative procedure. In each iteration step, the most favorable entry pressure curvature F * among all avail-able displacement scenarios N c is determined by and the corresponding region combination is invaded by gas. All allowed combinations of regions that are not yet occupied by the invading gas are generated, and the associated entry pressure curvatures are calculated for each region combination i by the three-phase MS-P equation given by Equation (2). The iterations are terminated when F * ðr go Þ cannot be determined, and the algorithm proceeds with the next gasoil capillary pressure. An individual pore geometry extracted from Bentheim sandstone is adopted to illustrate the iteration procedure we Figure 5: Numerical example illustrating the three-phase configuration of gas (yellow region), oil (red region), and water (blue region) after the second iteration step in the extracted pore space from Bentheim sandstone.
have developed to compute the three-phase fluid configurations for a given set of capillary pressure radii. The computed geometrically allowed oil-water interfaces are shown in Figure 1, and the valid initial oil-water configuration before gas invasion is shown in Figure 2. For the gas-oil and gaswater systems, all the geometrically allowed interfaces and extracted BR, LR, and MR regions that are used to determine the possible gas configuration candidates are shown in Figure 3. As described previously, the number of possible gas invasion candidates is 3612 for the first iteration. However, the corresponding bulk-oil and oil-layer configurations for each of these gas invasion candidates are also determined, which lead to an increased number of possibilities. The most favorable fluid configuration change as determined by Equation (3) in the first iteration is simply gas invasion into BRð3Þ, which occurs at the gas-oil entry curvature F * ðr go Þ = 0:0624 in pixel units. The resulting three-phase fluid configuration in the pore space after the first iteration is shown in Figure 6.
In the second iteration, the initial fluid configuration is updated with the gas occupying BRð3Þ. Thus, there are four individual bulk regions, BRð1Þ, BRð2Þ, BRð4Þ, and BRð5Þ, considered for gas invasion. The possible connected BR-LR -MR configuration candidates are given by (1) fBRð1Þ, LRð1Þ, BRð2Þg (2) fBRð1Þ, LRð1Þ, BRð2Þ, LRð2Þ, MRð1Þg   Candidates 1 and 6 were also considered in the first iteration step, whereas candidates 2-5, which all involve layer regions that adjoin the existing gas configuration as well as the other bulk regions, are new candidates in the second iteration step because of the updated gas configuration. Mixed regions can only be added to one of the individual bulk regions, BRð4Þ, which leads to 2 4 − 1 = 15 new possibilities. Adding the possible mixed region combinations to the above six connected BR-LR-MR configuration candidates lead to the following new configuration candidates in each case: ð1Þ 0, ð2Þ0, ð3Þ0, ð4Þ2 3 − 1 = 7, ð5Þ2 2 − 1 = 3, and ð6Þ2 3 − 1 = 7. Further, eight oil layers that are surrounded by the existing gas and corner water are associated with eight mixed regions, which result in seven new possible displacements (because the layer displacements involving MRð1Þ are already included in the BR-LR-MR combination candidates 2 and 3). Hence, a total number of 50 gas configuration candidates are identified in the second iteration, which represents a significant reduction in the number of possibilities as compared to the first iteration, even though new displacements have become possible due to the present gas region in the pore space. It turns out that the favorable displacement in the second iteration is obtained by merging a mixed region with a bulk region, fBRð4Þ, MRð11Þg, as shown in Figure 5. The corresponding gas-oil entry curvature is F * ðr go Þ = 0:0628 given in pixel units. In the third iteration step, none of the identified gas invasion candidates satisified the entry condition as stated in Equation (3). At this stage, the procedure is terminated, and the algorithm proceeds with the next r go .

Comparison of Simulated and Analytical Solution Results in Regular Star-Shaped Pores
We validate the method developed to compute three-phase capillary pressure curves and associated fluid configurations for gas invasion by comparing simulated results against analytically determined gas-oil capillary pressure curves in a straight tube with a regular star-shaped cross-section that contains three corners. The inscribed radius of the star shape is 40 pixel lengths, and the chosen half-angle of all corners is π/30 radians. We consider both spreading and nonspreading fluid systems. The interfacial tensions in the spreading system are given as σ go = 0:01 N/m, σ ow = 0:02 N/m, and σ gw = 0:03 N/m, whereas in the nonspreading system (C S,o = −0:005 N/m), the corresponding values are σ go = 0:01 N/m, σ ow = 0:02 N/m, and σ gw = 0:025 N/m. The oilwater contact angles considered are θ ow = 0°(strongly water-wet system) and θ ow = 40°(weakly water-wet system), and the corresponding gas-oil and gas-water contact angles are computed from the Bartell-Osterhof equation. The analytical solution for the three-phase capillary pressure curves in this type of pore space was derived by van Dijke et al. [21] and Helland et al. [8], and it is adopted for the comparison against our developed combinatorial, semianalytic method. As illustrated in Figure 7, the simulated gasoil capillary pressure results (data points) for gas invasion processes originating from different constant oil-water interface radii r ow exhibit an excellent agreement with the analytical solution (lines). The simulated results also demonstrate 7 Geofluids some basic displacement scenarios. In the spreading system, as shown in Figures 7(a) and 7(b), oil layers form if the pore space is filled with oil and water initially. Oil layers may also form in the nonspreading system if the oil-water interface radius is small enough, such as the cases of r ow = 2 and 8 pixel lengths for the strongly water-wet condition and for r ow = 2, 8, and 16 pixel lengths for the weakly water-wet condition. However, three-phase displacements in which gas displaces both oil and water simultaneously are most common within the nonspreading system, such as for r ow = 16 pixel lengths for the strongly water-wet condition and for r ow = 28 pixel lengths at the weakly water-wet condition.

Three-Phase Capillary Pressure Curves in 2D
Rock Images Simulations of three-phase capillary pressure curves and associated three-phase configurations of oil, water, and gas are performed for gas invasion processes originating from different initial oil and water saturations by applying the developed combinatorial, semianalytic method directly in a segmented 2D SEM image of Bentheim sandstone. The resolution of the image is 0.204 μm and it contains 1134 × 761 pixels and 293 separate pore geometries. The computed porosity and permeability values are 0.17 and 0:76 × 10 −12 m 2 , respectively. The image permeability is computed by solving Poisson's equation directly in the identified pore spaces using the MATLAB Partial Differential Equation Toolbox. In the simulations, we employ the same set of interfacial tensions and contact angles as given in the model validation part of this work.

4.1.
Three-Phase Distribution of Gas, Oil, and Water. Figure 8 shows three-phase fluid configurations at selected gas-oil capillary pressures during gas invasion at a constant oilwater capillary pressure for the strongly water-wet spreading  : Gas-oil capillary pressure as a function of total liquid saturation in a tube with star-shaped cross-section for gas invasion originating from different constant oil-water interface radii (given in pixel lengths) r ow = 40 (yellow), 28 (black), 16 (green), 8 (blue), and 2 (red). Solid curve: analytical solution; open circle: semianalytical solution. 8 Geofluids system. Oil layers between bulk gas and corner water form at small gas-oil capillary pressures because gas displaces only oil; these oil layers, despite becoming thin, continue to exist at low oil saturations. However, in the majority of the pore geometries occupied by both oil and water, gas invasion results in the formation of oil layers in some of the pore corners while gas-water interfaces form in other pore corners. Figure 9 shows corresponding results for the strongly water-wet nonspreading system, in which oil layer formation in most cases is prohibited due to the larger gas-oil contact angle. Therefore, gas invasion occurs mainly by three-phase displacements in which gas displaces both oil and water in the individual pore geometries, thus leading to the formation of gas-water interfaces. However, while a few thick oil layers are apparent in Figures 9(c) and 9(d), for gas invasion processes at larger oil-water capillary pressures in the nonspreading systems, oil layers form more frequently and may continue to exist in a larger range of gas-oil capillary pres-sures. The simulated fluid configurations appear to agree qualitatively with three-phase distributions of gas, oil, and water recently obtained experimentally in water-wet Bentheim sandstone for spreading and nonspreading systems using microtomographic imaging techniques [14]. Selected three-phase fluid configurations from simulations carried out for the weakly water-wet (θ ow = 40°) spreading and nonspreading systems are presented in Figures 10  and 11, respectively. The results show that the pore surface wettability affects the three-phase configurations significantly. In the weakly water-wet systems, many more oil layers are observed compared with the ones obtained in the strongly water-wet systems. Figure 10(a) indicates that gas invasion into water-filled pore geometries occurs with simultaneous invasion of oil layers separating the invading gas and the displaced water. Such displacements occur more frequently for spreading oil under weakly water-wet conditions in pore geometries for which the oil-water capillary entry (e) r go = 11 (f) r go = 2:5 Figure 9: Three-phase fluid configurations at selected gas-oil capillary pressure radii (given in pixel lengths) during gas invasion at a constant oil-water capillary pressure radius, r ow = 5 pixel lengths (P cow = 1:9608 × 10 4 Pa), in a strongly water-wet nonspreading system. Yellow region: gas; red region: oil; blue region: water.
(a) r ow = 14, r go = 9:5 (b) r ow = 5, r go = 4:5 Figure 10: Three-phase fluid configurations for two gas invasion processes originating from different initial oil and water saturations in a weakly water-wet spreading system. The gas-oil and oil-water capillary pressure radii are given in pixel lengths. Yellow region: gas; red region: oil; blue region: water.
10 Geofluids pressure is slightly higher than the oil-water capillary pressure. As shown in Figure 10(b), more thin oil layers continue to exist at low oil saturations than in the strongly water-wet spreading system. Many thin oil layers also continue to exist after gas has started to invade water filled pore geometries. In the results for the weakly water-wet nonspreading system shown in Figure 11, more oil layers are also observed than in the strongly water-wet nonspreading system.

4.2.
Three-Phase Capillary Pressure Curves. The saturation paths for the simulated gas invasion processes in the seg-mented rock image are presented in Figure 12. These results show that gas displaces all the oil phase first at the early stages of gas invasion and that the presence of gas-water interfaces do not contribute to a significant displacement of water, except in the regions of low oil saturation. Hence, in all the cases, the oil-water capillary pressure can be described as a function of only the water saturation, unless the oil saturation is low. This threshold oil saturation value depends on the rock wettability (strong/weak water-wet) and spreading/nonspreading system. Figure 12 suggests that for the strongly water-wet formations, there is practically no water (a) r ow = 4, r go = 14 (b) r ow = 4, r go = 11:5 Figure 11: Three-phase fluid configurations for two gas invasion processes originating from different initial oil and water saturations in a weakly water-wet nonspreading system. The gas-oil and oil-water capillary pressure radii are given in pixel lengths. Yellow region: gas; red region: oil; blue region: water. 11 Geofluids displacement by the pore space-invaded gas regardless the spreading characteristics of the system. However, in the nonspreading systems, shown in Figures 12(c) and 12(d), the water saturation decreases slightly during the gas invasion process due to the formation of gas-water interfaces, whereas in the spreading systems, shown in Figures 12(a) and 12(b), the saturation paths are almost constant isowater-saturation lines due to the oil layer formation. In the weakly water-wet systems shown in Figures 12(b) and 12(d), an increase of oil saturation is observed during the initial stages of the gas invasion process originating at relatively high water saturations. In these cases, the favorable displacement is determined as simultaneous invasion of bulk gas surrounded by accompanying oil layers in the pore spaces that are initially completely water-filled. However, these oil layers cease to exist at relatively low gas-oil capillary pressures. The wettability also affects the saturation paths at low oil saturations. In the weakly water-wet systems (Figures 12(b) and 12(d)), the saturation paths are affected by gas invasion into water-filled pore geometries while thin oil layers still exist in the spreading system and few small pores are still occupied by oil in the nonspreading system.
The two-phase oil-water drainage capillary pressure curves for the same sandstone Bentheim image are shown in Figure 13, and the three-phase gas-oil capillary pressure curves for the gas invasion processes are shown in Figure 14. For the three-phase case, the gas-oil capillary pressure curves collapse practically to a unique curve as a function of the total liquid saturation in the major part of the three-phase region. Therefore, the gas-oil capillary pressure curves can be described with high accuracy as a function of only the gas saturation. However, for low to intermediate gas saturations in the nonspreading systems, Figures 14(c) and 14(d), the simulated gas-oil capillary pressure curves are located within a narrow band. This behavior is attributed to the formation of gas-water interfaces when gas invades pore shapes which are saturated with both oil and water. Hence, the presence of water in the pore corners, along with the initial oil-water capillary pressure, affects the gas-oil capillary entry pressure. Such three-phase fluid displacement is associated with a different capillary entry pressure from the oil layer formation process in which gas displaces only oil at a two-phase gas-oil capillary entry pressure. The presence of corner water may also change the invasion order of the pore space compared to that of the two-phase gas-oil displacements.
At low oil saturations, the gas-oil capillary pressures seem to depend strongly on two fluid saturations in all cases, except for the strongly water-wet spreading system. This behavior occurs because gas has started to invade water-filled pores while oil still exists in the porous medium. For the nonspreading systems, the oil is located in few small pore spaces at low oil saturations, and gas invasion results in the formation of gas-water interfaces as described above (see Figures 14(c) and 14(d)). However, for the weakly water-wet spreading system, this changed saturation dependency of the gas-oil capillary pressure is most significant, as shown in Figure 14(b). In this case, oil layer formation occurs most frequently, and thin oil layers are observed to remain present at very low oil saturations even after the injected gas has started to invade pores occupied by water, thus leading to the nonunique gas-oil capillary pressure curves for low oil saturations presented in Figure 14(b).

Summary and Conclusions
A semianalytical pore-scale model is developed to simulate three-phase capillary pressure curves and fluid configurations for gas invasion processes in 2D segmented rock, which are initially fully saturated with oil and water. The pore spaces identified in the rock images are represented as cross-sections of straight capillary tubes. The three-phase fluid configurations occurring in the highly irregular pore spaces during gas invasion are modeled under uniformly water-wet systems. The model computes gas-oil, oil-water, and gas-water interfaces for each set of the three capillary pressures by moving a circle for each fluid pair in the two opposite directions along the pore/solid boundary, such that the contact angle is defined at the front circular arcs. Geometrically allowed interfaces are determined for each fluid pair based on circle intersections in the three separate cases. The resulting three sets of arc menisci are then grouped into boundary elements of regions, which are classified and combined systematically to allow for all possible gas invasion scenarios, both at the gas-oil capillary entry pressure and higher pressures. In all cases, the three-phase MS-P equation is used to identify the favorable displacements among all possible scenarios.
The semianalytical method is validated against analytically determined three-phase capillary pressure curves in an individual star-shaped pore space, and an SEM image of Bentheim sandstone is taken as input to the developed model for simulating the three-phase fluid configurations and capillary pressure curves under uniformly water-wet conditions and three-phase fluid systems with nonspreading and spreading oils. Overall, the simulated fluid configurations and capillary pressure curves depend strongly on the wetting condition and the spreading coefficient. The specific conclusions obtained in this study are summarized as follows: (1) Spreading systems clearly display oil-layer behavior, especially in a high oil-water capillary pressure state. Nonspreading systems display less oil layers, although some thick oil layers exist when the gas-oil capillary pressure is relatively low (2) Pore surface wettability affects the fluid configurations significantly. Numerous oil layers are observed in the weakly water-wet spreading system, even at low oil saturations, and more oil layers are also identified in the nonspreading system compared with the strongly water-wet systems (3) For both spreading and nonspreading systems, the simulated three-phase fluid configurations exhibit the same qualitative behavior with recently published three-phase fluid distributions obtained experimentally by microtopographic imaging techniques in Bentheim sandstone [14] (4) The simulated three-phase gas-oil capillary pressure depends strongly on two saturations at low oil saturations. This is because gas displaces oil and water simultaneously in this region.
(5) The gas-oil capillary pressure is described well as a function of only the gas saturation in the major part of the three-phase region, despite the formation of gas-water interfaces at low gas saturations in the nonspreading systems.
The proposed semianalytical model allows us to look at three-phase capillary entry pressure and its associated 13 Geofluids three-phase fluid configuration in realistic pore geometries under arbitrary wetting conditions, and thus, one possible future extension of this work is to combine it with the most recent generalized pore network modeling [25] to overcome the simplification assumed in its invasion percolation.