Generalized network modeling of capillary-dominated two-phase flow

We present a generalized network model for simulating capillary-dominated two-phase flow through porous media at the pore scale. Three-dimensional images of the pore space are discretized using a generalized network -- described in a companion paper (https://doi.org/10.1103/PhysRevE.96.013312) -- that comprises pores that are divided into smaller elements called half-throats and subsequently into corners. Half-throats define the connectivity of the network at the coarsest level, connecting each pore to half-throats of its neighboring pores from their narrower ends, while corners define the connectivity of pore crevices. The corners are discretized at different levels for accurate calculation of entry pressures, fluid volumes and flow conductivities that are obtained using direct simulation of flow on the underlying image. This paper discusses the two-phase flow model that is used to compute the averaged flow properties of the generalized network, including relative permeability and capillary pressure. We validate the model using direct finite-volume two-phase flow simulations on synthetic geometries, and then present a comparison of the model predictions with a conventional pore-network model and experimental measurements of relative permeability in the literature.

using a density functional approach [8,9], have been applied to immiscible flows in porous media [10][11][12][13]. Particle-based methods such as the lattice-Boltzmann method have been used to compute absolute and relative permeabilities [14,15], capillary pressure [16,17], interfacial area [18], and relative permeability for different rock types [19][20][21]. Another approach is to apply mesh-free methods such as smoothed particle hydrodynamics that, in addition to the computation of miscible flow and dissolution [22], have been employed to study twophase flow through porous media [23]. The advantage of direct methods is that the solid and fluid interfacial boundaries can be modeled accurately. Other features include the ability to study the effect of viscous forces, including the impact of flow rate and viscous coupling [21,24,25].
Direct simulations of two-phase flow, however, are computationally expensive. Capturing layer flow through pore-space crevices, for instance, requires a high-resolution mesh. Therefore, layer flow may not be captured accurately when using a coarse mesh that is usually required to make the simulations practical on larger images [13]. Direct simulations are more expensive at lower capillary numbers, where the total duration of the simulations is larger but small time steps should be taken to resolve capillary waves and local instabilities (Haines jumps and snap-off) [1]. For most subsurface processes, flow occurs at very low capillary numbers and the flow domain can have heterogeneity at different scales. Moreover, flow simulations may need to be run many times over a representative elementary volume that is several orders of magnitude larger than the grid resolution to study the effect of input parameters-such as pore structure, contact angle, flow rates, fluid viscosities, and initial conditions-on the macroscopic properties, and to quantify the effect of uncertainties in these parameters. A computationally efficient, and yet accurate, method to perform sensitivity studies to quantify the effect of these parameters is needed.
One solution to this computational challenge is to use a multi-stage upscaling approach. In the first stage, direct highresolution simulations on smaller system sizes can be used to obtain the equations required to describe flow in individual pores and throats. In a subsequent stage, these equations can be used to simulate flow through a coarser-scale network representation of the void space to obtain its macroscopic properties.

Network modeling of two-phase flow
Pioneered by Fatt [26], pore-network models have evolved as an important tool for studying flow through porous media. They have been used to, for instance, investigate the effects of viscous [24,27], wettability, and capillary forces, which control pore-scale configurations of fluids [3,28], discussed in more detail below.

Effect of wettability
Rock adhesion forces, which give rise to a contact angle between a fluid-fluid interface and a solid wall, have a major impact on pore-scale displacement mechanisms. Flow of the wetting phase through crevices of the void space (wetting layers), although slow, can lead to snap-off and hence trapping of nonwetting phases residing in the centers of the void space. This is in contrast to nanometer-thick wetting films that are stabilized by molecular forces in strongly water-wet media, which can have a significant impact on apparent contact angle and capillary pressure [29,30] but their contribution to fluid conductivity is negligible [31]. Wetting layers, on the other hand, provide a connected conduit for the wetting-phase flow down to low saturations [31,32].
Wetting layers can be modeled when using network elements with angular cross sections, including fractal roughness models [33], grain boundary pore shapes [34,35], squares [36,37], and triangles [31,38]. In most current network models, a shape factor, G, defined as the ratio of the crosssectional area to the perimeter length squared, is used to assign the shape of pores and throats. The shape of the cross section-a circle, square, or scalene triangle-is chosen such that it has the same shape factor as the 3D image of the porous medium [38,39].
The wettability of rock surfaces can change due to the adhesion of surface-active components of the oil to the solid surface [40]. The degree of wettability alteration depends on the composition of the oil and water, the mineralogy of the solid surface, and the capillary pressure imposed during primary drainage [40][41][42]. The average wettability of a fluid rock system can be measured using core-flood experiments [43]. More recently, contact angles have been measured directly in situ using micro-CT imaging [44][45][46][47]. Pore-network modeling has been used to link the pore-scale description of wettability to bulk measurements as well as to study the trend in recovery with wettability [48,49]. Network models allow the incorporation of these in situ measurements of contact angle and its history dependence on a pore-by-pore basis.

Effect of viscous forces
Most quasi-static network models impose a single capillary pressure over the entire network. This is used to define the fluid configuration in each element and the corresponding volumes of each phase. Fluid interfaces, residing between pores and throats or in their crevices, move and change their configuration as the capillary pressure changes. Filling of individual elements is assumed to take place over a much shorter time than the duration of the displacement process: this occurs in the form of Haines jumps, piston-like filling, layer collapse, or snap-off events. This is a valid assumption in most cases since typical macroscopic capillary numbers (C A = μU D σ , where μ is the fluid viscosity, U D is the Darcy velocity, and σ is the interfacial tension) in petroleum reservoirs are very low, in the range of 10 −6 to 10 −10 [2]. This represents a significantly large ratio of capillary pressure to viscous pressure drop across a single pore, usually more than two orders of magnitude. A single pore-filling event normally occurs in fractions of a second, as observed using fast x-ray imaging or acoustic measurements [50][51][52]. In contrast, it may take several days to years for a displacement process to be completed at a given location in a natural setting.
The viscous pressure drop ( Φ) can play a significant role when the flow rate is high (for example, in near-well-bore flow in hydrocarbon reservoirs), in near-miscible displacements with a low interfacial tension, or when the length ( x) of the representative elementary volume needed to obtain averaged properties becomes large. The viscous pressure drop scales as Φ = U D μ x/K, where K is the effective rock permeability, while the capillary pressure, P c , is not related to flow rate, U D , or system size, x. Capillary pressure scales as P c = 2σ/r c ∼ σ √ φ/K, where σ is the interfacial tension, r c is the mean radius of interfacial curvature, and φ is the porosity of the rock [1]. When the ratio of viscous pressure drop to capillary pressure is high, macroscopic flow properties, such as relative permeability, can be functions of flow rate, leading to a Darcy law where flow rate nonlinearly depends on the pressure gradient [31].
To model two-phase flow when the viscous pressure cannot be ignored, the quasi-static assumption can be relaxed to consider the viscous pressure drop as a perturbation to the local capillary pressure along the length of the system. A higher level of sophistication can be achieved by using a dynamic network model that also incorporates the effect of changes in fluid volumes on the viscous pressure. In these dynamic models, the volume of each phase in each element is updated using the flow rates from the computed pressure field. Usually a very small time step is required so that the configuration of fluid interfaces does not change significantly [53][54][55][56][57]. On the other hand, considering the viscous forces as a perturbation to the local capillary pressure will not cause a significant compromise in the computational efficiency. In this approach the pressure fluctuations associated with rapid changes in fluid configurations during filling events are ignored, but the viscous pressure drop is used as a perturbation to assign a nonuniform local capillary pressure across the system. We adopt this perturbative approach in our model. constructing the element connectivity and shapes from the analysis of micro-CT images [38,62,63]. However, Bondino et al. [64] showed that the algorithm used to generate the network and assign its parameters using current network extraction algorithms can have a significant impact on the predicted macroscopic properties. In fact, improving the predictability of pore-network extraction and flow modeling remains an active research topic [65][66][67][68]. The problem is that the algorithm used to assign the conventional network model parameters-specifically pore and throat lengths, shape factors, and volumes-cannot be independently verified and hence the model primarily relies on calibration to predict experimental measurements. Therefore, there is little confidence that a network calibrated using a small set of benchmark experiments can predict the properties of different rock types or different displacement scenarios reliably. To avoid this problem, the geometry of individual pores in the network model should be as close as possible to the real system.
We have presented a generalized network representation of the void space that eliminates the need for intermediate parameters, specifically shape factors and pore and throat lengths, used to describe the network elements [69]. Instead, the parameters required for network modeling are extracted using a medial-axis analysis of the void space and direct singlephase flow simulation. Each pore is subdivided into half-throats and further into corners. The parameters approximating each corner-corner angle, volume, and conductivity-are directly extracted from the underlying micro-CT image at different discretization levels and exported as tabulated data to the network flow simulator.
The emphasis of this paper is on the formulations used to describe capillary-dominated two-phase flow through this generalized network representation of the void space, and on the prediction of multiphase flow properties-capillary pressure and relative permeability. We present the algorithms used to track the fluid-fluid interfaces and incorporate the effect of contact angle that describes the wettability of the fluid-rock system. Gravity and viscous forces are treated as a perturbation to compute the local capillary pressures throughout the network.

II. GENERALIZED NETWORK FLOW MODELING
The network model can be viewed as an upscaled alternative to direct two-phase flow models. It is comprised of three main components: (a) parametrization and tracking of fluid interfaces in each element as the local capillary pressure (P c ) changes during a displacement cycle, (b) tracking fluid phase distribution and connectivity, and (c) computation of fluidphase saturations and conductivities for individual elements and for the whole system.
The flow simulations are designed to represent typical coreflood experiments. Initially, the whole void space is assumed to be filled with water, and a second fluid, oil, is injected from one side of the network to obtain the fluid configurations at the end of this primary oil-invasion. The simulations are continued for a second cycle, water flooding, in which the water pressure is increased to replace the oil. The flow simulations can be continued for a third cycle by injecting oil to obtain, among other properties, wettability indices. 1. An illustration of two fluids (highlighted in red and light blue) occupying a void space that is discretized into (a) pores and half-throats and further into (b) corners. The thick solid black lines show the boundaries between the pores (the throat surfaces). The solid blue lines show the partitioning of the pore space into half-throats. The dashed blue lines show the boundary between the corners of each half-throat.
The following sections present the details of our network flow model, which involves tracking the fluid configurations within the void space (see Fig. 1) and upscaling their flow properties. An overview of how the void space is represented using half-throat corners is discussed in Sec. II A. Fluid configurations in a corner are described using their interfaces with other fluids in the corner and their connectivity to fluids in the surrounding pores and throats; this is explained in Sec. II B. Section II C presents the details of how fluid interfaces are represented and tracked. In Sec. II F, we describe the computation of entry pressures required for fluid interfaces to change configuration. These include filling of centers of pores and throats, and growth and collapse of water and oil layers in edges of void space corners. In Sec. II G, we discuss how displacement events (changes in fluid configurations) affect fluid phase connectivity: forming disconnected phases and newly connected phases (coalescence), and exposing new invasion paths for subsequent displacement events. Finally, in Sec. II H, we discuss the assignment of fluid volumes and conductivities that are then used in the computation of the fluid saturations and relative permeabilities of the network. Figure 1 shows how the void space is divided into pores, half-throats, and their corners. We use a watershed segmentation of the distance map (the distance of any point in the void space to the nearest solid) of the underlying image to divide the void into pore regions bounded by throat surfaces [69]. In this algorithm, the voxels located at local maxima of the distance map are assigned as the pore centers and the remaining voxels are successively assigned to the pore to which one of their adjacent voxels (the voxel with the largest distance map) is assigned. We further divide these pore regions into half-throats and their associated corners. Every point in the void space is uniquely assigned to a pore, a throat, and a half-throat corner. We refer to the half-throat corners simply as corners in what follows.

A. Description of the network elements
The corners are parametrized at different discretization levels, i = 1-3, which are obtained during network extraction. Each discretization level contains the corner void space outside maximal spheres (spheres inside the void space that are tangent to the solid boundary at two or more points) with radius larger than R i (see Fig. 2). In this paper, we choose R 1 = R p , R 2 = R t , and R 3 = 0.7R t .
We define a local coordinate, (x,y), for each corner: x represents the distance from the throat surface, measured along the throat line, and y represents the distance from the throat line measured along the corner's medial axis (Fig. 2). The term sagittal plane is used to represent planes parallel to the corner medial axis. The coronal plane is used to represent the planes perpendicular to a corner side (boundary between the corner and the solid walls; see Fig. 3) that passes through the pore or throat center. The corner edge is defined as the line where the two sides of the corner meet, and the vector along the corner edge that connects the throat surface to the corner's axial cross section at the pore center is called the edge vector, e.
The maximal-sphere radii and cross-sectional areas are used to compute discretization level depths (H i ) measured along the corner sides, and half angles (γ i ), as illustrated in Fig. 3 and discussed in Appendix A.
The extracted corner parameters for each discretization level consist of the maximal-sphere radius (R i ), cross-sectional area in the corner's axial plane (A i ), volume (V i ), and flow (g q ) and electrical (g e ) conductivities. These properties come from the analysis of the underlying pore-space images [69].
R i , H i , and γ i define the geometry of the corners and are used to semianalytically track the location of fluid interfaces through the pore space as the local capillary pressure changes. V i and g i define the volumes and conductivities of the corners, which are used to compute fluid saturations and relative permeabilities of the network.

B. Connectivity and fluid configurations
How we define fluid configurations and assign their connectivity is shown schematically in Fig. 4. The fluid configurations during a flow simulation are modeled by considering four flow paths for each corner: (1) in the throat center, (2) in the pore center, (3) in the corner edge, and (4) sandwiched between the corner edge and corner center. Different fluid configurations are then constructed by marking each flow path as filled by oil (α = o) or water (α = w), subject to the following rules: (i) Corners share the same throat center, and the fluid occupying a throat center occupies all its corner centers. Similarly, (ii) the fluid in the pore center is shared between all of the pore corners. (iii) Only oil layers are considered to be sandwiched between water in the edge and in the center of a corner. If two adjacent flow paths are occupied by different fluids, a fluid-fluid interface is allocated between them. A flow path can grow or shrink in size if the interface separating it from other fluids in adjacent paths moves between different corner discretization levels, or along the same discretization level.
Corner edges and corner centers are considered adjacent flow paths. However, if there is an oil layer sandwiched between a water layer and water in the center, the two water phases are considered disconnected from each other. The interface between a pore center and a throat center is called a piston-like interface. An interface separating two fluid layers, or a fluid layer and a fluid in the center, is called a layer interface or simply a layer. Sections II C and II E discuss how these two interface configurations are represented.
Another factor controlling a fluid configuration is its connectivity with fluids in the surrounding pores. At the throat surface, the corners on either side, belonging to neighboring pores, are connected by definition. At a pore center, each corner does not necessarily connect to every other corner associated with the pore's half-throats, discussed next. This is where our approach differs from conventional network models (e.g., [39,49,70]), where it is assumed that at a pore all the corners are connected to each other.
At a pore center, we only allow a corner to be connected to one or two other corners belonging to different throats, based on their proximity. We first find the adjacent throats, up to two throats, whose throat lines have the smallest angle with the corner's y axis (see Fig. 2). Then, among corners of each of these throats, we find the adjacent corner that has the smallest angle between its y axis and the corner's y axis. If the angle is less than 60 • , the two corners are assumed to be connected to each other and are called adjacent corners. We only allow connections to corners in different throats: fluids in the corner edges of a throat are not directly connected to each other.
The water layers residing in each corner are considered to be connected to the water inside their adjacent corners. Similarly, oil layers are considered to be connected to oil phases in their adjacent corners. These adjacent fluids can themselves be in layer, piston-like, or single-phase configurations.
During network extraction, throats that are connected to the left and right sides of the image boundary are considered boundary throats. In this paper, for the sake of simplicity, we call the left side boundary the inlet and the right side the outlet. The connectivity of each fluid to outlet and inlet throats is obtained using a graph search through the adjacent flow paths containing the fluid. Fluids that are not connected to the outlet throats are considered trapped and do not change configuration. Similarly, the invading fluid connectivity, which is injected from the inlet throats, affects the order that the new interface configurations are formed, as discussed in Sec. II G. In addition, fluid connectivities affect the computation of their conductivity (discussed in Sec. II H) and computation of interface curvatures and entry pressures for different displacement events, which are discussed next.

C. Layer configurations
The location of a layer interface (h l )-that can be oil or water or a sandwiched oil layer-in a corner is defined as the distance of its interface contact line from the center of the pore or throat, which is computed along the side (solid wall) of the corner as illustrated in Fig. 5. A sandwiched oil layer, which has two interfaces (see Fig. 4), is described using its interface with the water in the center. The location of the other interface, between the oil layer and the water in the edge, is referred to as the water layer interface. A layer, if present, is tracked using its interface location, which at the throat surface 5. An illustration of the parameters used to describe the location, h l [Eqs. (1)-(3)], of a layer (light blue) located at level 2 of a corner at the throat surface, and its relationship with contact angle (θ l ), radius of interface curvature in the throat's axial plane, r l [Eq. (5)], and the interface distance, y l [Eq. (7)], from the corner center along its medial axis.
can reside at discretization levels 2 and 3. Near the pore center, however, layer interfaces can reside additionally at level 1; this can happen in a piston-like configuration as discussed in Sec. II E or due to the variations in a layer interface curvature in its sagittal plane (Sec. II D).
To uniquely describe a layer configuration, in addition to h l , we need to know the contact angle, θ l , between its interface and the solid wall. For simplicity, we define the layer contact angle as the angle between the layer interface and solid wall that is measured through the layer (l) that is closer to the corner edge. This definition is not the same as the conventional contact angle (θ w ) that is always measured through the denser (water) phase (θ w = θ l for water layer interfaces and θ w = π − θ l for oil layers). The contact angle is history dependent and can vary between a receding and an advancing contact angle (θ r and θ a , respectively) that are input into the network flow model.
The location, h l , of a layer interface (Fig. 5), residing in a corner discretization level i (Fig. 3) and for a given contact angle, θ l , is where H i [Eq. (A2)] is the corner level depth and γ i [Eq. (A1)] is its half angle. r l is the interface radius of curvature in the corner's axial plane. r l can be obtained from the local capillary pressure P c : Here κ l is the interface total curvature and r s is its radius of curvature in the sagittal plane of the corner, discussed in Sec. II D in detail.
To allow a unique assignment of the contact angle and interface location, we need to record the initial interface location as the layer configuration forms and track it as the simulation progresses. We initially set h l = H 4 − (near the corner edge; see Fig. 3) if the interface forms due to layer growth, and h l = H 1 + if the interface is left behind during a piston-like invasion of the corner center (see Fig. 4). is a small number, set to 10 −9 m in this paper.
We then use a multistage computation to find a unique solution for θ l and h l from Eq. (1) for a given capillary pressure or interface curvature. First we find the discretization level, i, in which the interface has been previously residing (see Fig. 6): At the throat surface, however, i can be 2 or 3 only. If Eq. (1), with this i and h l fixed to its previous location, gives a θ l in the range [θ r , θ a ], it is assumed that the location of the interface remains pinned at its original position and the computed θ l is accepted as the so-called hinging contact angle. Otherwise, the contact angle is fixed to θ r when the layer pressure decreases, or to θ a when the layer pressure increases, and Eq. (1) is used to compute the new interface location.
If the computed interface location using Eq. (1) does not correspond to its previous level based on Eq. (4), the interface will move to a new discretization level that satisfies Eq. (4), or will be pinned at their boundary. If (a) γ 2 < γ 3 , the interface is assumed to jump to the new corner level and Eq. (1) is invoked with the new discretization level corner angle to recompute the new interface location. However, if (b) γ 2 > γ 3 , the interface gets pinned at h l = H 3 before moving to the next discretization level. In this case, we first check if the interface has a stable position at h l = H 3 assuming γ = γ 2 with a contact angle between [θ r − γ 2 + γ 3 , θ a ]. If Eq. (4) has a solution in this range, the interface will be pinned at h l = H 3 ; otherwise it will move to the new corner level and Eq. (1) is invoked with the new level corner angle to recompute h l .
The same algorithm is used for a piston-like interface to find the location of its tailing layers that reside at the pore center but not at the throat surface (see Fig. 8), when they move between different discretization levels i = 1-3.
h l and the contact angle, θ l , can be used to obtain other interface parameters (see Fig. 5), including the width of the layer interface, W l , and the distance of the interface center from the throat center, y l : To obtain the layer cross-sectional area, A l , we first compute the cross-sectional area for a hypothetical layer located at h l but with a contact angle of zero, A h , and then add the effect of contact angle, A θ : where θ * = θ + γ i , r * l = r l cos θ * , and R h is the radius of the maximal sphere tangent to the solid wall at the interface location h l (see Fig. 5), Note that for the case of oil layers sandwiched between water in the edge and in the center, A l , obtained using Eq. (8), includes the cross-sectional area of the oil as well as the water in the corner.

D. Layer description in a corner's sagittal plane
The layer interface radius of curvature in a corner's sagittal plane (r s ) is needed in Eq. (3) to relate the radius of curvature in the corner's axial plane (r l ) to the local capillary pressure. We estimate it from the tangent vectors, s 1 [Eq. (B1)], to the interface in the corner's mid-sagittal plane at a distance x = x 1 = x p /2 from the throat surface (see Fig. 7).
At the throat center, r s is estimated from Similarly, r s at the pore center is estimated from  12) and (13)] in the sagittal plane of a corner (highlighted using horizontal yellow stripes). The sign of r t s is negative while r p s is positive in this case. j is the adjacent corner in the same pore and p 2 is the adjacent pore. e is the tangent vector to the corner edge and s is the unit vector tangent to the layer interface in the corner's sagittal plane.x andŷ are the unit vectors along the corner local coordinates, x and y.

E. Piston-like configuration
A piston-like configuration, Fig. 8, refers to a fluid-fluid interface separating a pore center from a throat center; these interfaces are often called terminal menisci, since they block the center of the pore space [1].
The curvature of a piston-like interface, κ pl , is controlled by the interface contact angle with the solid walls and the configuration of its tailing layers residing near the corner edges. It is obtained by writing a force balance equation on the interface: Corner's Coronal plane where A t x is the throat total cross-sectional area at a distance x from the throat surface, Eq. (A7). The summations ( ) are performed over all the throat corners, c = 1-n c , where n c is the total number of corners in the throat. l stands for the layer in the corner, c, that is adjacent to the fluid in the center. If the layer does not exist in the corner, its area (A l ) and arclength W l are set to zero and h l = H 4 . Finally, β is the angle between the corner side plane and the throat line-the line connecting the throat center to the pore center (see Fig. 8), computed using Eq. (A6).
We compute the piston-like curvature at three points along the throat line, x, at the throat center (x 0 = 0), at the pore center (x 2 = x p ), and in between at x 1 = x p /2. The interface curvature is assumed to vary linearly from x = x t to x 1 , and from x = x 1 to x p : The sign of κ pl is considered positive when the interface is curved toward the pore center, hence A location is assigned to each piston-like interface, defined as the distance of the interface contact line from the throat surface, x pl . Similar to the layer interface location h l , x pl is assigned as soon as a pore or a throat center is filled by a fluid, forming the piston-like configuration, and tracked as the simulation progresses. The interface is assumed to remain pinned (fixed x pl ) as long as the invading phase pressure is below P e (x pl ) [Eq. (16)]. Once the invading phase pressure surpasses P e (x pl ), the interface location is updated by solving Eqs. (15) and (16) for x pl (replacing x with x pl ).
The computed interface curvature and location are used to assign the pore and throat entry pressures and to compute the fluid volumes and conductivities, which are described in the following sections in more detail.

F. Computation of entry pressures
As discussed in Sec. II B, we consider four flow paths for each corner: in the throat center, in the pore center, in the corner edge, and sandwiched between the corner edge and its center. A displacement is a change in the occupancy of a flow path and involves a reconfiguration of fluid interfaces. The entry or threshold pressure, P e , is the pressure required to overcome the interfacial force as the interface passes through the flow path during this reconfiguration. P e is defined as the maximum relative pressure, the invading phase (subscript inv) pressure minus the receding phase pressure, encountered during each displacement:

Throat center entry pressure
For throats there are (a) one entry path from each of the two neighboring pores (piston-like displacement) and (b) one from each corner (snap-off displacement), if they contain the invading phase and are connected to the inlet.
a. Piston-like displacement. The throat entry pressure by piston-like invasion is computed by solving Eqs. (14) and (16) at the throat center, by iteratively changing the interface curvature, κ pl , using the Newton-Raphson method.
b. Snap-off entry pressure. The throat snap-off pressure is approximated from the lowest interface curvature (measured toward the center or receding phase) that the interface experiences as it moves from its initial location toward the throat center meeting the interface of other throat corners. It is approximated as the minimum of the k l , Eq. (3), obtained from (a) solutions of Eq. (5) with h l varying from its initial location to h l = 0 and (b) the solution of Eq. (7) for y l = 0.
Note that both h l and y l are history dependent due to contact line pinning (see Fig. 6). The effect of contact line pinning is considered in the calculation of h l and y l , in all the equations that require the computation of interface location using Eqs. (1)-(4).

Pore-filling pressure
A pore center can be invaded from any of its throats that contain the invading phase. The pore center entry pressure, P p e , is the smallest of all the threshold pressures (σ κ t max ) that are obtained for the throats, t, from which the invading phase can fill the pore. κ t max is the largest curvature (positive toward the throat center) encountered by the interface in throat t, as it moves toward the pore center.
To obtain κ t max , we compute the interface curvature using Eq. (14) at two locations and choose the maximum: (a) when the interface is midway between the throat and the pore and (b) when it reaches the pore center. Note that, except in an unstable configuration where the phase in the throat is the nonwetting phase (see Fig. 4), the maximum curvature is expected to occur at the pore center. However, the interface can also be pinned between the throat and the pore center due to the expansion of the throat, especially for contact angles close to 90 • where cos(θ + β) can have its minimum value between the pore and throat centers [see Eq. (14) and Fig. 8].
The pore-filling pressures take into account the configuration of wetting fluid in adjacent throats and hence provide an accurate representation of pore filling for complex geometries. The effect of the fluids in the adjacent throats in Eq. (14) is accounted for through the incorporation of layer interface tangent vectors, s l , which depend both on the corner connectivity to corners of its adjacent throats, and on their fluid occupancy (see Appendix B). This contrasts with current network models that use empirical formulas to compute pore-filling pressures [49,70], which may be inaccurate and lead to, for instance, poor predictions of the amount of trapping [71], although recently more accurate models have been developed [68].

Oil layer collapse and growth pressures
When water invades a throat center, an oil layer can be left behind in the corner if it has a stable configuration-its collapse pressure is lower than the local capillary pressure. The oil layer will collapse once the local capillary pressure falls below its threshold collapse pressure. The threshold oil layer collapse pressure is obtained using a geometric criterion when it is continuous, connected to the oil phase from all its 9. An illustration of the parameters used to estimate oil layer collapse and growth pressures. W is the interface length in the corner's axial plane. A o and A w are the oil and water layer areas. adjacent corners. It is obtained by iteratively increasing the invading phase (water) pressure until the interfaces on either side of the oil layer join [36], either from the center, or from the sides, where h o l is the location of the interface between the oil layer and the water in the center of the throat and h w l is the location of the water-layer interface, both measured along the sides of the corner; y o l and y w l are the locations of the oil and water layer interfaces measured along the center line of the corner, as shown in Fig. 9.
However, if an oil layer is not continuous from one side (i.e., it is adjacent to a corner that does not contain oil, either in its center or in its crevice) it is expected to collapse more easily. The entry pressure, P ol e = σ max κ ol , for such scenarios is estimated based on a thermodynamic criterion [72], by writing a force balance on the interface in the normal direction to the corner's axial plane: This equation is solved using the Newton-Raphson method to obtain κ ol . Oil layers can grow in oil-wet corners if the entry pressure required for their growth is lower than the entry pressure for the invasion of the throat center by oil. Equation (20) is used to obtain the oil layer growth pressure, but with θ r o replaced by θ a o . This thermodynamic criterion is expected to make the layer growth more difficult compared to the geometric criteria [Eqs. (18) and (19)].

Water layer collapse pressure
Water layers form if at the maximum depth of the corner (h l, max ) the condition for their formation is satisfied (γ 3 + θ < π/2) and the center of the throats is filled by the oil phase, but they are assumed not to collapse.

G. Displacement sequence
During each flooding cycle, the invading phase pressure at the inlet, P inlet a , is increased relative to the receding phase pressure, P inlet r . The local invading phase pressure relative to the receding phase pressures, P a-r , throughout the flow domain is assumed to change proportional to the difference in advancing and receding phase viscous pressures, Φ a − Φ r : where Φ is the viscous pressure obtained by solving the mass conservation equation for each phase, discussed in Sec. II H. A displacement happens when P a-r surpasses the entry pressure, P e , for a receding fluid that is adjacent to an invading fluid connected to the inlet. Once a fluid is displaced, all the adjacent flow paths that contain the receding phase and are not trapped (i.e., are part of a cluster that is connected to the outlet) are considered for subsequent displacement events and their threshold entry pressures are (re)computed.
Any receding fluid that is part of a cluster not connected to the outlet is considered trapped. The curvature of trapped fluid interfaces are assumed to be preserved. This implies that the capillary pressures of trapped ganglia remain independent of the imposed capillary pressure at the inlet boundary.
Furthermore, if an adjacent flow path contains the invading phase and has been previously marked as part of a trapped ganglion, the ganglion is removed from the trapping list and is brought into capillary equilibrium with the invaded flow path (coalescence event). Coalescence events may require running a mini-imbibition cycle if the ganglion's local capillary pressure is higher than the system capillary pressure, or a mini-drainage cycle if the ganglion has a lower capillary pressure than the capillary pressure of the invaded flow path. First, the entry pressures for all the formerly trapped fluids, for filling by the mini-cycle's invading phase, are computed. In a mini-imbibition cycle, if the entry pressure is higher than the system capillary pressure, the flow paths are filled with the mini-cycle's invading phase. Once a flow path is invaded, the connectivity of the adjacent receding fluids in the mini-cycle is checked; if they are not connected to the inlet (invasion front), they are marked as trapped, forming smaller trapped ganglia.
The main displacement cycle is continued by successively increasing P a-r and displacing the receding fluids that are not trapped and are adjacent to the invading fluids, if their entry pressure is less than P a-r . This process is continued until the network saturation or capillary pressure reaches a desired userdefined limit.

H. Computation of relative permeability
To compute fluid saturations, permeabilities, and electrical resistivity of the network, we first need to calculate the volume and conductivities of fluids residing in the corners constituting the network.
The layer volumes and conductivities are computed by interpolating the tabulated corner volume and conductivities that are obtained during network extraction. The interpolations are done both at the pore and at the throat centers and then their average (harmonic for conductivities and arithmetic for volume) is taken.
The layer volume (V l ) and electrical conductivity (g e l ) in each corner are assumed to scale linearly with the corner crosssection area: where A l is given by Eqs. (8), (9), and (10).
To compute the layer flow conductivity in each corner, g q i , we assume that it scales with cross-sectional area squared: where . When there is a piston-like interface separating the fluids in the pore and throat centers, the volume and conductivity of the fluid in the throat are obtained by linear interpolation between the fluid volumes and conductivities of the levels 1 and 2, using the distance of the interface from the throat center, x pl , as the interpolation parameter: The equations above [Eqs. (22)-(24)] lead to cumulative fluid volume and conductivities, the volume and conductivity of all the fluids below the interface-toward the throat surface for piston-like configurations and toward the corner edge in layer configurations. To compute the area, volume, and conductivity of the fluids above the interface-toward the pore center for piston-like configurations and toward the throat center line in layer configurations-we simply subtract these volumes and conductivities from the single-phase corner volumes and conductivities. The flow conductivities are further scaled by the new fluid area relative to the area before this subtraction: where g q SP and A SP are the level 1 (single-phase) flow conductivity and cross-sectional areas. The exponent C A is set to 0.4, which is chosen such that the center-phase conductivity closely matches the direct simulation results presented in Sec. III A.
We compute a single flow rate for each throat and a single pressure for each pore. The conductivities of the fluids in the corners of each pair of half-throats are averaged to assign a conductivity to each throat. The averaging is performed by grouping the fluids in the corners into two categories, (a) those that are connected together through the fluid occupying the throat center, and (b) those that are not. For group (a), we first use an arithmetic sum of the conductivities of fluids in each half-throat, and then take the harmonic average of the two half-throat conductivities. For group (b), we first compute the harmonic average conductivity for each corner pair on opposite sides of the throat surface, to obtain the full-corner conductivity, and then add them together. Finally the two group conductivities are added together to compute the throat conductivities, g α t , for each phase (α = o,w) and for electrical current, e.
In addition to the rules above, for adding the corner conductivities to obtain a single value for the throat conductivity, we assume that layers that are not continuous from both sides do not contribute to the conductivity of the throat. In other words, if a layer is not continuous, including the tailing layers of piston-like configurations, from one side or from either side, its conductivity is not added to the throat conductivity. Our results, presented in Sec. III B, show that this exclusion of discontinuous layer conductivities is essential to predict the correct behavior of the relative permeability.
Once the individual throat conductivities are computed, the relative permeabilities of each phase, α, are obtained by solving for mass conservation in each pore, p, in the network: where t counts over all the throats connected to the pore p and g α t is the conductivity of the throat connecting the pore to its neighboring pore (subscript nei).
These equations are solved assuming a dimensionless pressure of 1 at the inlet and zero at the outlet nodes. The flow rate is then obtained by summing the flow rates entering the flow domain, which is the same as the sum of flow rates in the throats adjacent to the outlet. The relative permeability of each phase is obtained by dividing its flow rate by the single-phase flow rate, which is computed before starting the first cycle using the same boundary conditions for pressure. Finally, the computed pressures are scaled to correspond to the flow rate or capillary number assigned for each phase. In the results presented in this paper, for the sake of simplicity, we assume a low capillary number such that the effect of flow rate on the displacement sequence can be ignored.
The electrical resistivity of the network can be obtained using Eq. (26), but with g α replaced by g e and Φ α replaced by the electrical potential, Φ e .

III. VERIFICATION AND VALIDATION
In the following, we first use analytical approximations and direct simulation of two-phase flow through simple pore geometries to validate the network model on a pore-by-pore basis. Then, we compare the generalized network model predictions with a conventional network and core-flood experiments on sandstones. The aim of this comparison is to demonstrate the improvements achieved by the generalized network model in predicting relative permeabilities from micro-CT images of porous rocks.

A. Pore-by-pore validation
This section evaluates the accuracy of the network model in calculating the capillary entry pressures, fluid volumes, and conductivities on a pore-by-pore basis. Figure 10 shows the synthetic geometries used for this purpose, which include starand triangular-shaped geometries with different corner angles (described by Raeini et al. [69]), and different pore-throat contraction (R p /R t ) and aspect (L t /R t ) ratios, where L t is the distance between the two pore centers. The images are converted to 3D images similar to micro-CT scans with a voxel size of 1.6 μm, corresponding to a resolution of R p /δx = 18.75, where δx is the voxel size, which is also equal to the average grid-block size used in the direct two-phase flow simulations, described below. for all geometries except (e) and (f), for which L t /R t = 12.5 and 6.25, respectively. The pore-to-throat contraction ratio (R p /R t ) is 2.5 for all geometries except for (g) and (h), for which R p /R t = 1.5 and 3.5, respectively. Water is shown in light blue and red represents the oil phase. The interface is in a piston-like configuration in (a) and (b) and in a layer configuration in the rest.
The threshold capillary entry pressures for filling the middle throat by piston-like and snap-off events, and for the pore in the right side of the middle throat, are shown in Fig. 11. The generalized network model (GNM) results are compared to a conventional network model [62,70] and analytical approximations.
The analytical approximations (Anl) are obtained using the same equations presented in this paper but using the corner angles of the original geometry. However, the effects of interface curvatures in the corner sagittal planes are ignored in their calculation. The conventional network model (CNM) uses, in essence, the same equations as the GNM and analytical approximations. However, the number of corners in the CNM and the corner angles is not the same as in the original geometry [69]. The interface curvatures in the corner sagittal planes are also ignored in the CNM.
The GNM matches the analytical approximations for the entry pressures closely. The CNM underpredicts the pore entry pressures and overpredicts the snap-off entry pressures. This is because corner angles in the CNM, estimated using shape factors, differ from the original geometry [69]. Consequently, in the CNM, when the contact angle is 30 • the interface in the middle throat at the second cycle snaps off before forming a piston-like configuration and the throat piston-like entry pressures are not computed in the second cycle; the throat piston-like entry pressures for the CNM shown in Fig. 11 are taken from the first cycle entry pressures, with the receding contact angles set equal to the advancing contact angles. The GNM predicts a lower snap-off pressure compared to the analytical approximation, which is partly because the sagittal curvature is included in the GNM but not in the analytical approximations. To predict relative permeabilities accurately, in addition to entry pressures, the computations of fluid volumes and conductivities should be accurate too. We validate our computation of fluid volumes and conductivities by comparing the network model predictions with direct twophase flow simulations on these synthetic geometries, using the same contact angles as the network model.
The direct numerical simulation (DNS) is a volume-offluid-based finite-volume method. It uses the improved surface tension algorithm described in [73], and the pressure-velocity coupling and filtering algorithms described in [74]. An unstructured mesh, with grid blocks roughly the same size as the voxels used in the network extractions (1.6 um), is used to discretize the flow domain. The grid blocks away from the solid walls are cubic. However, near the solid walls they are deformed to align with the solid boundary, using the snappyHexMesh meshing tool from OPENFOAM [75]. An additional cell layer is added adjacent to the solid walls, inside the flow domain, so that wetting layers can be captured more accurately. The fluid densities and viscosities chosen for both fluids are the same, 1000 kg/m 3 and 0.001 Pa s, respectively. The interfacial tension is 0.03 N/m. The simulations are performed by first initializing the water layers from the images of corner discretization levels obtained during network extraction. In addition to the three discretization levels discussed in this paper, we run the simulations for a level halfway between discretization levels 1 and 2, called level 1.5, resembling a piston-like configuration (see Fig. 10).
A no-slip boundary condition is used on the solid walls of the direct simulations. The outlet boundary condition is zero gradient for velocity and indicator function, and a constant value for the dynamic pressure [74]. At the inlet, we have used a zero-gradient boundary condition for all the variables except velocity, for which a constant flow rate for each phase is assigned. The inlet boundary velocities are initially set equal to their adjacent internal cell velocities and then corrected so that the flow rate of each phase converges toward a desired value [13]. The chosen apparent velocities (flow rate divided by image cross-sectional area), for the simulations initialized with images of corner discretization levels 1, 1.5, 2, and 3, are q w = 0.6, 0.06, 0.06, and 0.03 mm/s for the water phase, and q o = 0, 0, 0.6, and 0.6 mm/s for the oil phase, respectively. Note that oil is not present in the system at level 0, and at level 1.5 it only occupies the pore centers. For the system sizes studied here, these velocities are sufficiently small that the flow can be considered capillary dominated. The capillary numbers based on the average pore-scale velocity (μu α /σ , where u α is the average velocity of each phase α = o,w), obtained from the simulation results, are all less than 9 × 10 −5 and 5 × 10 −5 for the oil and water phases, respectively (see Appendix D).
In summary, four steady-state direct two-phase flow simulations are run for each case, corresponding to four fluid saturations. Visualizations of the fluid configurations for a set of these simulations, when the flow is steady state, are shown in Fig. 10. In the following, we upscale the direct simulation results to compute the saturation and conductivities, for the voxels comprising the middle throat [13], and compare them with the GNM and the CNM saturation and conductivities.
A sensitivity study on the effects of mesh resolution, viscosity ratio, and flow rates-compared to the base case parameters discussed here-is presented in Appendix E, for the star-shaped geometry with the corner angle of 45 • and a contact angle of 30 • . Furthermore, Appendix F presents a visualization of the predicted interface location in DNS and the GNM for this geometry. The interface location is an intermediate parameter that relates the assigned capillary pressure at each stage of simulations to the computed saturation and conductivities that are the subject of the validation study that follows. Figure 12 shows a comparison between the GNM, CNM, and DNS results for the volume fraction of water in the middle throat as a function of curvature radius (r c ), as the system capillary pressure decreases during the second (water flooding) cycle.
The GNM results show a good match with the direct simulations for the computed water saturations. In the CNM simulations, with contact angles of 0 • and 30 • , all the throats are filled by snap-off, trapping the oil and leading to a fixed saturation as the imposed capillary pressure varies. This saturation is lower than the DNS results when the water phase is initialized from the level 2 corner images (the middle point in the DNS results in Fig. 12), which predict a stable layer configuration for these contact angles. This implies that the CNM underpredicts the water saturation of the wetting layers at a given capillary pressure or curvature radius. Note that the DNS results are only presented for contact angles less than 60 • . To keep the water layer stable for higher contact angles in the DNS, contact angle hysteresis similar to network models, should be implemented; this is considered a subject for future work.
Comparisons between the GNM, CNM, and DNS conductivities for oil and water in the middle throat are given in Figs. 13 and 14, respectively. The results, presented for the second (water flooding) cycle, are normalized by singlephase flow conductivities and plotted as a function of water saturation.
We have also compared our results with a generalized network model but with conductivities obtained using correlations presented in Appendix C; this method is called GNMCrl.
The GNM results originally underpredicted the water phase conductivities by a factor of 2. In all the results presented in this paper, we have multiplied the GNM conductivities at levels 2 and 3 by a factor of 2 to resolve this difference. Although the GNM conductivities are computed using direct single-phase flow simulations on corner images, the corner images are obtained ignoring the effect of sagittal interface curvature of the wetting layers. In addition, this difference can be explained, partly, by the differences in the boundary conditions used in DNS and the single-phase flow simulations used to compute the GNM conductivities. In the single-phase simulations, we have used a slip boundary condition-no mass and no momentum transfer between the corner voxels and the voxels in the center of throats. The DNS, however, incorporates a continuous velocity across the interface of the two fluids, which implies that there is a drag force between the two fluids that can increase the apparent conductivity of the water layer [25]. Further work is needed to incorporate the effect of viscous coupling in the network model [76,77]. The CNM water conductivities follow the trends of the DNS results, but they underpredict the water saturation at their end points, as discussed above.

B. Micro-CT image of porous rocks
In this section we present a set of flow simulations on micro-CT images of porous rocks to demonstrate the improvements achieved by the incorporation of the generalized network model parameters and to show that we can have a reasonable estimation of macroscopic properties of relatively simple rocks with straightforward choices of input parameters.
We simulate capillary-dominated two-phase flow through networks of a Berea and a Bentheimer sandstone, based on 1000 3  throats. Other network properties are given in [69]. During network extraction, the corner images obtained for discretization level 3 were not connected from the inlet boundary to the outlet, so their conductivity could not be obtained using direct simulations. Higher resolution images are needed to obtain the level 3 conductivities using direct simulation. To overcome this problem, we have extrapolated the discretization level 2 conductivities to obtain the conductivities at level 3: g q 3 = (R 3 /R 2 ) 4 g q 2 , and ϕ 3 = (R 3 /R 2 ) 2 ϕ 2 , where R 3/2 = R 3 /R 2 and ϕ = A,V ,g e . This implies that we effectively use two levels to discretize the corners of the micro-CT images.
A uniform intrinsic contact angle of 45 • is used in all simulations and Morrow's hysteresis model III [78] is used to compute the receding (oil-injection) and advancing (waterinjection) contact angles: 3 • and 46 • , respectively.
Berea sandstone is known to contain 7-8% clay minerals [79][80][81]. Bentheimer sandstone is reported to have 1-3% clay minerals by weight [82][83][84]. In this paper we add clay porosities of 6.4% and 1.5% of the total image volume to each network, respectively. Note that there is uncertainty involved with these values, but these can be considered reasonable choices based on the available data from the literature.
The Berea simulation results are presented in Fig. 15 and are compared to the CNM [62,70] predictions and to experimental measurements of relative permeability for oil-water flow by Fulcher et al. [85] and Oak and Baker [86] and for CO 2 -water flow by Akbarabadi and Piri [87]. FIG. 15. Drainage (top) and imbibition (bottom) relative permeabilities for Berea sandstone, obtained using different network modeling approaches and compared with the experimental measurements indicated [85][86][87]. A uniform intrinsic contact angle of 45 • is used in the simulations.
The GNM relative permeabilities show a reasonable match with the experimental data for both rocks. The oil relative permeability, however, is slightly overpredicted by the GNM. Although the quality of the match for the oil relative permeabilities can be improved, for instance, by adjusting the area exponent used to compute oil conductivities in Eq. (25), a more rigorous investigation of the effect of various parameterssuch as viscous coupling, assignment of contact angles, and uncertainties in experimental measurements-is needed to fully resolve this anomaly. The CNM clearly underpredicts the water saturation, which is also the case for the synthetic geometries presented in the previous section. Further adjustments to contact angles are needed to obtain a better match for the water relative permeability in the imbibition cycle of the CNM. Moreover, these results show that the GNMCrl, which uses conductivities based on correlations obtained from direct two-phase flow simulations, produces a similar behavior as the GNM formulation. Further work is needed for a rigorous validation of the generalized network model, with contact angles measured from the micro-CT imaging of two-phase flow [45,53] and an FIG. 16. Drainage (top) and imbibition (bottom) relative permeabilities for Bentheimer sandstone computed using GNM, GNMCrl, and CNM, and compared to experimental measurements [38,88,90]. assessment of the effect of uncertainties, for instance, in image segmentation and unresolved and clay porosity.

IV. CONCLUSIONS AND FUTURE WORK
We have presented a generalized network model for simulating two-phase flow through micro-CT images of porous media. The network represents a coarse discretization of the pore space with properties obtained from upscaling of direct simulation of single-phase flow through the corners of the underlying image and correlations based on direct two-phase flow simulations.
This workflow allowed us to model pore-scale events considering the complexity of the pore geometries encountered in natural porous media and to validate our computations using direct simulation of two-phase flow and experimental measurements of relative permeability. The results show that accurate computations of corner volumes, conductivities, and assignment of corner connectivity is critical in the prediction of relative permeabilities from micro-CT images of porous media.
However, there are other sources of uncertainty in the predictions of pore-scale models, related to image resolution, clay volume, and corner connectivity, for instance. To fully resolve these sources of uncertainty, the network model can be extended by incorporating viscous coupling and dynamic effects. Its validation can be extended by considering more complex geometries and different wettability distributions. The validation can include, in addition to fluid saturation and flow conductivity, interface location, wetting layer area, and electrical conductivity, for instance. Overall, the results presented in this paper show that accurate network modeling, informed by direct two-phase flow simulations and potentially by multiphase pore-scale imaging [50,51,92] and core-scale measurements of flow and electrical properties, offer a predictive framework for linking the pore-scale processes and fluid and rock properties to their macroscopic counterparts, helping to answer open questions that cannot be addressed by these methods individually.

ACKNOWLEDGMENTS
The authors are grateful to TOTAL for the financial support, fruitful exchanges, and permission to publish this work.

APPENDIX A: CORNER SHAPES AND LOCAL COORDINATES
The void space in the generalized network model is reconstructed from the parameters extracted during network extraction for each discretization level, i = 1-3, as shown in Figs. 3 and 17.
First, the maximal-sphere radius and cross-sectional areas are used to compute discretization level depths (H i ), measured along corner sides and half angles (γ i ) [69]: where δ is the difference operator: δA = A i − A i+1 and δR = R i − R i+1 . A 1 and A 2 are defined at the throat surface and are assumed not to change between pore and throat centers. A 1 is defined at the pore center.
FIG. 17. An illustration of the parameters used to reconstruct the shape of half-throat corners. Note that the cross-sectional areas at levels 2 and 3 are assumed constant along the lengths of the corner.
The discretization level lengths, L i (see Fig. 17), are estimated as follows: The line connecting the throat center to the pore center is called the throat axis and is considered the x axis of the corner. Therefore, the local x coordinate of the pore center is We also compute the vector along the edge of the corner, e i (see Fig. 7): The pore distance, x p , and pore and throat radii, R i , are used to compute the angle, β, between the throat axis and the corner sides (Fig. 8). β is computed at the throat center (x t = 0), at the pore center (x p ), and at distance x = x 1 = x p /2 from the throat surface: The total throat cross-sectional area at any point along the throat line is computed as where the summation is performed over all the throat corners and A c x is the corner cross-sectional area at a distance x from the throat surface, where A h is obtained using Eqs. (9) at h = x p −x x p H 1 + x x p H 2 , located in the discretization level i = 1.

APPENDIX B: INTERFACE TANGENT VECTORS
The interface tangent vector, s 1 , at x = x 1 , between the pore and throat centers, is assumed to be parallel to the corner edge vector: The interface tangent vectors at the pore and throat centers are also needed to compute the interface curvature in the pistonlike configuration. The interface tangent vector at the throat center, s t , is obtained by averaging the unit edge vectors, e [Eq. (A5)], at either side of the throat surface: For a case that the adjacent corner interface is in a layer configuration, the interface tangent vector at the pore center, s p , is obtained by averaging the unit edge vector with the adjacent corner's unit edge vector: whereẑ is the unit vector normal to the throat center line (x) andŷ. However, if the adjacent corner's interface is in a pistonlike configuration, the interface is expected to bend toward the throat center line. When the piston-like interface is at the pore center, we assume s p is parallel to the throat center line (x): Equations (B3) and (B4) account for cooperative pore body filling [32,68], which implies that the curvature of the terminal menisci, positive toward the throat center, decreases as more of its surrounding throats are filled with the invading phase, making filling more favorable.

APPENDIX C: CORRELATIONS FOR COMPUTING CORNER CONDUCTIVITIES
In this section we present a set of correlations that we have used to estimate the corner conductivities in the GNMCrl formulation. These correlations are considered a faster but approximate alternative to direct single-phase flow simulations. The effect of viscous coupling is ignored in these equations.
The corner electrical and flow conductivities (g e i and g q i , respectively), at discretization levels i = 2 and 3, are obtained using the following equations, which are approximations to the correlations used by Valvatne and Blunt [70] for corners of throats with triangular cross-sections: To estimate the conductivity of corner centers, we first estimate them using the corner parameters at the throat surface: Then, these are corrected for the effect of the expansion of the half-throat cross-sectional area between the throat surface and the adjacent pore centers, assuming that the maximal-sphere radius changes linearly: where R p/t = R p R t and δR p/t = R p −R t R t are the expansion ratio and the relative expansion of the maximal-sphere radius from the throat center to the pore center. Finally, the corner conductivities at level i = 2 are added to these conductivities to obtain the level 1 (single-phase flow) conductivities: g e 1 = g * * e 1 + g e 2 , (C7) (C8)  Figure 18 presents the capillary numbers (N c ) obtained from direct simulations for the flow rates discussed in Sec. III A. The definition of capillary number based on the Darcy velocity, U D = Q/A, is not applicable here since the cross-sectional area, A, is not constant along the system. Therefore, the capillary numbers are presented in terms of the average pore velocity, u α , of each phase, α = o,w, in the middle throat: N c = μu α /σ . These capillary numbers are sufficiently low that viscous effects do not affect the shape of the interface (see Appendix F).

APPENDIX E: SENSITIVITY OF DIRECT SIMULATIONS AND THE GENERALIZED NETWORK MODEL
In this section, we evaluate the effect of fluid properties, flow rates, and the convergence of DNS and the GNM with mesh resolution. The simulations are performed on the star-shaped geometry with a corner angle of 45 • and for a contact angle of 30 • . The results are presented for a base case with the same resolution (R t /δx = 18.75) and same fluid and rock properties as discussed in Sec. III A. Two other mesh resolutions R t /δx of 12.5 and 25 are considered, using both the GNM and DNS (Fig. 19). To obtain an estimation of the dependence of the flow conductivities on the flow rate, the DNS results are additionally  19. Convergence of the middle throat radius of curvature, r c = σ/P c , and water and oil conductivities normalized by singlephase flow conductivity, g rw and g ro , respectively. The suffixes 12, 18, and 25 represent the simulations with mesh resolutions of R t /δx = 12.5, 18.75, and 25. The suffix 2q w shows the simulation whose water flow rate is twice the base case value. The suffixes μ w /2 and 2μ show the simulations with water viscosities twice and half the base case, respectively.
presented for a case with a water flow rate twice the base case. Finally we investigate the effect of viscosity ratio, by presenting one simulation in which the water phase viscosity is twice the base case, and one with a water viscosity that is half the base case value.
The difference between the results for the three mesh resolutions is less than 5%. These results together with the validations of the direct simulation presented by Shams et al. [73] and Raeini et al. [74] show that DNS captures the flow conductivities and the interfacial curvature accurately. Although the GNM results provide an acceptable match with the DNS, they show slightly larger variations with image resolution.
Increasing the water flow rate by a factor of 2, in DNS, increases the oil conductivity by less than 5%. The water-oil viscosity ratio has a larger effect on the oil conductivity: halving the water phase viscosity increases the oil conductivity by roughly 5% for this case. The effect of the flow rate and viscosity is less significant for the water phase saturation and conductivity. Larger corner angles or higher relative flow rates and viscosity ratios are expected to change the fluid conductivities more significantly [76,77]. A more rigorous quantification of the effect of these parameters is a subject for future work. Nevertheless, these results show that when the fluid viscosities and the driving force are in the same range for both fluids, the assumption of a single fluid conductivity for each phase (hence ignoring the viscous coupling effect) does not introduce a significant error in the calculations. Figure 20 presents a visualization of the interface location obtained from the GNM and compares it with the fluid phase locations obtained using DNS for the star-shaped geometry with a corner angle of 45 • , contact angle of 30 • , and a mesh resolution of R t /δx = 18.75. The direct simulation is initialized with the level 3 corner image that corresponds to the results with S w = 0.11 in Fig. 19. The GNM results are obtained by assigning a final capillary pressure of P c = 2450 Pa for the second cycle, chosen to be equal to the DNS capillary pressure. The predicted water saturation in the GNM is 0.097, a discrepancy of 12%.

APPENDIX F: VISUALIZATION OF INTERFACE LOCATION
The interface location [Eq. (1)] in the GNM, for a given contact angle, is computed analytically from the corner angle and axial interface curvature [Eq. (3)]. This implies that if FIG. 20. A visualization of the interface location, the transition between water (blue) and oil (red), in DNS (top) and axial slices through the pore throat centers in DNS (middle) and GNM (bottom). the axial curvature and corner angle are accurate the interface location will be accurate too. The accuracy of assigned corner angles is studied in Raeini et al. [69]. The accuracy of the interface curvatures are presented in this paper in the plots of radius of curvature as a function of wetting phase saturation (Figs. 12 and 19). Figure 20 shows that the water layer area is smaller at the throat surface compared to the pore centers, in both DNS and the GNM, which is due to the variations in the sagittal interface curvature. This quality check is not possible for the CNM since there is no one-to-one correspondence between corners of the CNM and the underlying geometry.