A Discrete-Domain Approach to Three-Phase Hysteresis in Porous Media

An inherent challenge with multiphase flow processes in porous media, such as water-alternate-gas invasions for oil recovery or CO2 injection for geologic storage, is that fluid displacements are irreversible, that is, they are history-dependent and exhibit hysteresis (Spiteri & Juanes, 2006). Hysteresis can be rate-dependent or rate-independent. Here, we examine rate-independent hysteresis in three-phase systems based on an energy landscape exhibiting metastability and barriers. Thus, we approach fluid displacements quasi-statically as a series of small changes of pressure or saturation. This is a reasonable approximation for slow displacement at pore scale where capillary forces typically dominate over viscous forces and gravity (Hilfer & Øren, 1996). At darcy (or, core) scale, hysteresis in two-phase systems emerges as the difference between drainage and imbibition capillary pressure curves (Pαβ(Sβ)-curves), where capillary pressure is the difference in phase pressures p between nonwetting (α) and wetting (β) phases (i.e., Pαβ = pα − pβ), and Sβ is wetting-phase saturation. Collectively, the shape of Pαβ(Sβ)-curves depends on porous structure, pore-scale fluid displacements, saturation history, and wetting state. Morrow (1970) interpreted Pαβ(Sβ)-curves thermodynamically as a sequence of reversible and irreversible pore-scale fluid displacements, termed “isons” and “rheons,” respectively, that describe transitions between capillary equilibrium states (i.e., local energy Abstract We present a discrete-domain approach to three-phase displacements and hysteresis in porous media. In this method, constrained energy minimization leads to evolution equations for local saturations that describe a wide range of three-phase displacements, including pressureand saturationcontrolled displacement with or without preservation of one of the defending phases. Under action of global saturation constraints, irreversible displacements lead to significant fluid redistribution, as well as abrupt fluctuations of both the three-phase saturation paths and the corresponding capillary pressures. These features are a consequence of Haines jumps with cooperative behavior that occur at pore scale in three-phase systems. The method is a fast and convenient way to investigate hysteresis behavior of three-phase displacement in porous media. As free energy is an extensive property, the framework links pore and core scales and provide a means to achieve upscaled three-phase displacements for higher-order hysteresis loops, which rarely is obtained in time-consuming three-phase measurements or pore-scale simulations.

For three-phase flow, hysteresis and irreversible displacements are more complex and much less explored. This is due to several reasons. Displacement behavior vary significantly with immiscible and near-miscible fluid systems, the spreading behavior of oils on gas/water interfaces, and different wetting orders of the three phases (Alhosani et al., 2019;Hui & Blunt, 2000;Keller et al., 1997;Khishvand et al., 2016;Scanziani et al., 2020; van Dijke & Sorbie, 2002). At the pore scale, double displacements (Helland & Jettestuen, 2016;Helland et al., 2019;Keller et al., 1997;Khishvand et al., 2016;Øren & Pinczewski, 1995), where a continuous phase displaces a second isolated phase that displaces a third continuous phase, and multiple displacement chains (Jettestuen et al., 2021;Scanziani et al., 2020; van Dijke & Sorbie, 2003), where the continuous phases displace a chain of several isolated fluid ganglia, make three-phase flow different to two-phase flow. Further, three-phase experiments are time-consuming and challenging, while three-phase pore-scale simulations are computationally demanding (Helland et al., 2019;Jettestuen et al., 2021). They both rely on having established the relevant two-phase saturation history prior to the investigation. Finally, a two-phase displacement occurs with saturation change in one of two directions, while a three-phase saturation path can take infinitely many directions. Thus, measuring a sufficient amount of P αβ (S β )-curves for main and nested hysteresis loops in three-phase systems at relevant conditions is hardly achievable.
The standard approach to deal with hysteresis of P αβ (S β )-curves is to use correlations, equipped with a hysteresis loop logic model and a trapping model relating endpoint saturations (e.g., Helland & Skjaeveland, 2004;Land, 1968;Lenhard, 1992;Lomeland & Ebeltoft, 2013;Skjaeveland et al., 2000;Spiteri et al., 2008, and others). While this is practical, it requires several fitting parameters that lacks a firm physical basis. Other approaches seek to eliminate hysteresis by adding new state variables based on Minkowski functionals (like interfacial area) in the incomplete P αβ (S β )-function (Hassanizadeh & Gray, 1993;McClure et al., 2018), whereas emerging theories circumvent P αβ (S β )-curves and hysteresis completely (Hansen et al., 2018). Cueto-Felgueroso and Juanes (2016) developed a two-phase discrete-domain model for hysteresis that describes Haines jumps as irreversible transitions among metastable states in a rugged energy landscape. The method envisions the porous medium as a discrete set of compartments, each with their own given energy function and local saturation. This facilitates simulations of pressure-controlled displacement, which calculates saturations for a prescribed series of capillary pressures, and saturation-controlled displacement, which mimics displacement at infinitesimally low rate and calculates capillary pressures for a prescribed series of global saturations. The latter displacement mode enforces a global saturation constraint that leads to fluid redistribution among compartments (cooperative behavior) and pressure fluctuations, consistent with saturation-controlled pore-scale simulations (Helland et al., 2017). Discrete-domain methods can be viewed as upscaling procedures from pore to core scales, as compartment energies are extensive properties (Cueto-Felgueroso & Juanes, 2016). Similar methods have also been used in other scientific fields (Dreyer et al., 2010(Dreyer et al., , 2011Fraternali et al., 2011;Puglisi & Truskinovsky, 2005).
In this work, we present a discrete-domain approach to characterize three-phase displacements with hysteresis in porous media. We accomplish this by formulating the interfacial free energy as a sum of energy contributions from the respective phases and derive evolution equations for each phase that relates each individual energy contribution to its phase saturation and phase pressure. By imposing constraints that couple the equations together and ensure the sum of all saturations is one, we obtain the desired conditions for equilibrium that relate the capillary pressures to the appropriate change in interfacial free energy. Hence, the developed approach is easily adaptable to an arbitrary number of fluid phases. As the two-phase method of Cueto-Felgueroso and Juanes (2016) employs an evolution equation for the wetting-phase saturation that directly relates capillary pressure and interfacial free energy, our approach is a required generalization of their method to investigate hysteresis in porous media containing three (or more) fluid phases. With the developed method, we show that metastability and the inclusion of various saturation constraints gives rise to irreversible three-phase displacements with cooperative behavior and fluid redistribution, accompanied by fluctuations of capillary pressures and saturation trajectories.

Method
We consider a porous medium saturated with three continuous fluids, gas (g), oil (o), and water (w), that are connected to their respective fluid reservoirs. Assuming the multiphase system is isothermal with incompressible fluids of constant composition, the thermodynamic potential for Gibbs free energy is where p α is pressure and V α is volume of phase α. Here, we introduce the phase-specific Helmholtz free energy, F α , which models the energy contribution from phase α to the fluid/fluid and fluid/solid interfacial free energies. This splitting strategy borrows ideas from pore-scale modeling of three-phase displacement where interfacial properties, like interfacial tension and contact angle, are formulated as the sum of individual phase contributions (Helland & Jettestuen, 2016;Helland et al., 2019;Jettestuen et al., 2021). We also assume this Helmholtz energy is a function of its own phase volume only, that is, Following Cueto-Felgueroso and Juanes (2016), let us now divide the porous medium into N compartments with equal pore volumes V p . The global (or, average) phase saturation is  where s α,i is local saturation of phase α in compartment i. We introduce phase-specific compartment energy densities f α,i as functions of their own saturations, Minimization of G with respect to s α,i yields an equation per saturation per compartment. The resulting set of 3N equations is uncoupled and neglects that the sum of phase saturations must equal one in all compartments. To handle this problem, still by treating the compartment saturations s α,i as independent variables, we introduce N constraints that couple the three saturations together in each compartment: In the numerical model, we set the left-hand side of Equation 3 equal to a small positive number ϵ. The method of Lagrange multipliers is a suitable approach to solve such constrained minimization problems. Thus, we introduce the function where α = g, o, w and i = 1, …, N. We solve Equation 5 for given pressures p α and energies f α,i to find equilibrium states of the compartment saturations, using the explicit Euler method for the iteration-time steps. We assume Equation 5 has been made dimensionless by scaling with a characteristic pressure. As we calculate equilibrium states, this scaling will only impact the rate of convergence while the solutions remain the same.
The compartment saturations reach equilibrium states when . Equation 5 implies that this occurs when the capillary pressures for all fluid pairs αβ = go, ow, gw satisfy When the third phase is absent or immobile, ds α,i = −ds β,i , and we can write Equation 6 as This is the equilibrium states for a two-phase system, consistent with the single saturation approach of Cueto-Felgueroso and Juanes (2016).
The energy functions f α,i (s α,i ) contain concave and convex segments (peaks and valleys) to describe metastability in a rugged energy landscape. Stable and metastable states satisfy Equation 6. In addition, the interfacial free energy must belong to convex segments with local minima of Gibbs free energy. In the two-phase case these valleys satisfy A peak in f i represents an energy barrier that must be overcome in a transition from the current to the next, adjacent, metastable state. In the three-phase system, while the second derivatives of the phase contributions f α,i to the interfacial free energy must satisfy the following conditions (see Appendix A): ,  2  2  2  2  2  2  ,  ,  ,  ,  ,  ,  2  2  2  2  2  2  , , , Equation 7 is the equilibrium locus, that is, the states that correspond to metastable compartment saturations for any phase pressure combination.
The interfacial free energies f i describe stability of multiphase displacements in compartments representing irregular pore structures. At the pore level, interfaces in two-phase systems stabilize at the entrance of pore throats during drainage and in wider pore openings during imbibition according to Laplace's equation, P αβ = 2σ αβ C αβ , where σ is interfacial tension and C is mean interface curvature. These metastable states belong to valleys in the interfacial energy f i that describe reversible displacements. Fast and irreversible interface displacements (Haines jumps) occur through wide pore openings during drainage and through narrow pore throats during imbibition, representing transitions across barriers separating adjacent energy valleys in f i . Consequently, the local saturation intervals for both reversible and irreversible displacements will in general differ for drainage and imbibition, leading to hysteresis. Thus, we can interpret the local saturation interval of an energy barrier in f i as the saturation interval shared by irreversible jumps in opposite directions (drainage and imbibition).
Three-phase systems achieve metastability when Equations 6 and 7 are satisfied. The behavior of gas invasion through pore spaces containing oil in pore centers and water in pore corners is very similar to the two-phase scenario described above. However, for double displacements, where gas displaces an oil slug that in turn displaces water, metastable states can also correspond to configurations where either the trailing gas/oil interface or the leading oil/water interface is located in a wider pore opening while the other interface is located at the entrance of a pore throat, depending on the size of the oil slug (or, local oil saturation) (Helland et al., 2019). Therefore, the local saturation intervals of reversible displacements (energy valleys) and interfacial jumps (energy barriers) described by f i is in general different for two-phase systems and three-phase systems. We also note that for two-phase systems, the way we split f i into individual phase contributions f α,i and f β,i is somewhat arbitrary as it is sufficient to use evolution equations for one of the phases in this case (Cueto-Felgueroso & Juanes, 2016). However, in three-phase systems, the splitting of f i into phase contributions f g,i , f o,i and f w,i must be done with care to ensure the model describes the desired case. In future work, we will construct the energy contributions f α,i for various three-phase systems and fit the model to pore-scale simulations and pore-scale experiments.

Three-Phase Displacement Protocols
Equation 5 offers the opportunity to investigate various three-phase displacement modes, like invasion under preservation of a constant global saturation, or invasion controlled by prescribed global saturations. This is an important capability because many experimental setups for core-scale three-phase flow facilitate injection of one phase and production of a second phase, while the saturation of the third phase is constant (Bradford & Leij, 1995Lenhard & Parker, 1988). Such experiments can also examine how the presence of a disconnected phase (e.g., fluid ganglia) impacts storage or recovery of other fluids in porous rock.
Other applications, like three-phase saturation-path tracking, require control of two of the saturations. This feature is particularly useful for predicting three-phase capillary pressure curves along a saturation trajectory obtained from core-scale measurements of three-phase relative permeability curves (Alizadeh & Piri, 2014;Zolfaghari & Piri, 2017).
Introducing a global saturation constraint for one or more phases α, will couple the saturation Equation 5 for the different compartments together. Despite global saturations are fixed under such constraints, significant redistribution of local saturations among compartments can occur, possibly accompanied by three-phase capillary pressure fluctuations. Here we apply Equation 5 to a range of different three-phase displacement protocols.

Gas-Pressure Controlled Displacement (PG)
This case prescribes all phase pressures p g , p o , and p w . Gas invasion (retraction) occurs by increasing (decreasing) p g in small specified steps Δp after each converged state of Equation 5, consistent with standard methods for measuring static capillary pressure curves in rock samples. The N compartment constraints (3) must hold for all iteration times. Thus, Combining Equations 5 and 9 yields expressions for the Lagrange multipliers that we update in every iteration step:

Gas-Pressure Controlled Displacement With Constant Global Oil Saturation (PGSO)
This case prescribes phase pressures p g and p w . As for the previous case, gas invasion (retraction) occurs by increasing (decreasing) p g stepwise after each equilibrium state. We preserve global oil saturation S o by enforcing constraint (8) on the oil phase, while treating the oil pressure p o as the associated Lagrange multiplier. Together with Equation 3 we now have N + 1 constraints that must hold for all iteration times.
By inserting Equation 5 into Equations 11 and 9, we obtain and respectively. Equations 12 and 13 constitute N + 1 equations that we solve for the Lagrange multipliers p o and λ 1 , …, λ N in every iteration step. Displacement protocol PGSO conforms to standard methods for measuring static three-phase capillary pressure curves under circumstances where one of the defending phases (oil) cannot exit the rock sample.

Gas-Saturation Controlled Displacement (SG)
This case prescribes phase pressures p o and p w . Gas invasion (retraction) occurs by increasing (decreasing) global gas saturation S g in small specified steps Δs after each converged state of Equation 5, while we calculate p g . For each target gas saturation, we initialize the compartment gas saturations as s g,i ±Δs and adjust the compartment oil and water saturations correspondingly, s o,i ∓Δs/2 and s w,i ∓Δs/2. Obviously, the compartment saturations at equilibrium will differ from this initial configuration. Equation 8 enforces the target global gas saturation, so that p g becomes the associated Lagrange multiplier. Thus, .
Equations 13 and 14 constitute N + 1 equations that provide solutions for p g and λ 1 , …, λ N in every iteration step. The saturation-controlled displacement mode SG mimics three-phase flow where gas invasion or withdrawal occurs at infinitesimally low rate.

Gas-Saturation Controlled Displacement With Constant Global Oil Saturation (SGSO)
This case prescribes water pressure p w and calculates p g and p o . Gas invasion (retraction) occurs by increasing (decreasing) the global gas saturation S g in small specified steps Δs after each equilibrium state, while global oil saturation S o is fixed. For each target gas saturation, we initialize compartment saturations as s g,i ±Δs and s w,i ∓Δs. We enforce Equation 8 on both the gas and oil phases, so that both p g and p o are Lagrange multipliers. Together with Equation 3 we now have N + 2 constraints. Equations 12-14 constitute the N + 2 equations from which we calculate the Lagrange multipliers p g , p o , λ 1 , …, λ N in every iteration step. Displacement protocol SGSO is an approximation to three-phase flow where gas invasion or withdrawal occurs at infinitesimally low rate while the oil saturation is preserved. We note that this strategy can also be used to predict three-phase capillary pressure curves along any saturation path given a set of target oil and gas saturations.

Results
We apply the method to elucidate three-phase displacements and hysteresis behavior for the previously described displacement protocols, using synthetic and oscillatory energy density functions to capture energy barriers and metastability. To this end, we make advantage of the energy functions employed in the single-saturation, two-phase approach (Cueto-Felgueroso & Juanes, 2016): where s i is the wetting-phase compartment saturation, ω represents the strength of the wetting state, K reflects the number of minima in each compartment, and c i , i = 1, …, N, describe the barriers between energy minima that are responsible for the hysteresis. If ω varies among compartments, the system has a fractional wetting state, whereas a constant and uniform ω (as we use in this work) describes a uniform wetting state. Cueto-Felgueroso and Juanes (2016) showed that refining Equation 15 into a higher number of domains N, which increases the number and density of metastable states, will dampen pressure fluctuations in saturation-controlled displacement and result in smoother capillary pressure curves while the level of capillary pressure and hysteresis persists. Our numerical examples demonstrate the same behavior in three-phase systems.
In the multiple saturation approach, we construct phase-specific energy contributions f α,i (s α,i ) that ensure the equilibrium solutions (6) for fluid pair αβ are consistent with the solutions from the single-saturation approach for two-phase systems using Equation 15. We accomplish this with where the parameters for each fluid pair αβ satisfy K = K α = K β , ω = ω β − ω α , and c i = c α,i + c β,i (for even integers K). As we explore displacements constrained by stepwise saturation changes Δs, we specify a saturation interval [s min , s max ] that we require all compartment saturations to be located within. If a compartment saturation reaches one of these limits, that compartment is static and excluded from the further saturation evolution based on Equation 5 and Lagrange multipliers. Thus, the number of possible metastable states decreases toward the end of the displacement process, possibly accompanied by more severe pressure fluctuations. Here, we use s min = 0.001 and s max = 0.999. This approach also offers a means to include residual saturations that may vary with compartment. The simulations presented here achieve equilibrium states when all compartment saturations change by less than 10 −8 during an iteration step.
First, we demonstrate the discrete-domain method on a small number of compartments for a water-wet state where oil is intermediate-wet phase and gas is nonwetting phase (Hui & Blunt, 2000; van Dijke & Sorbie, 2002). This implies ω w > ω o > ω g and c w,i > c g,i > c o,i > 0, i = 1, …, N. We set N = 3, K = 4, ω w = 11, ω o = 6, ω g = 1, . We simulate a two-phase oil/water hysteresis loop consisting of primary drainage and imbibition by applying Equation 5 on the two phases, using Δp = 0.1 for pressure-controlled displacement and Δs = 0.005 for saturation-controlled displacement. Figure 1 shows that pressure-controlled Haines jumps create irreversible saturation jumps at constant pressure. Such an event is a jump in one of the compartment saturations across an energy barrier separating two metastable states of Gibbs free energy, while the other compartment saturations stay in their energy valleys. Overall, this generates a stepwise staircase shape on the P ow (S w )-curves, where steep and smooth curve segments represent reversible displacements within the same energy valleys. Haines jumps in saturation-controlled mode occur as abrupt, irreversible pressure drops (drainage) or increases (imbibition) at constant global saturation. While these events also create a jump across an energy barrier in one compartment, the other compartment saturations move in the opposite direction to maintain the global target saturation. Hence, fluid redistribution among compartments occurs. Reversible displacements follow the same branches of the P ow (S w )-curves as the pressure-controlled displacements. Overall, the P ow (S w )-curves obtain nonmonotonic sawtooth shapes which we also observe in the evolution of compartment saturations. The peaks (during drainage) and valleys (during imbibition) of the capillary pressure fluctuations in saturation-controlled displacement coincide with the capillary pressure levels in pressure-controlled displacement, which indicates high precision in the simulations of both displacement modes.
We proceed with corresponding three-phase displacements by introducing gas at S w = 0.4 after the pressure-controlled imbibition. Figure 2 presents results for a gas invasion-retraction hysteresis loop using both PG and SG displacement protocols. The P go (1 − S g )-and P gw (1 − S g )-curves have identical shapes and differ only by the constant P ow . They display similar behavior as the two-phase case. Energy barriers connecting the equilibrium locus of Gibbs free energy represent instabilities. In three-phase systems we cannot construct these barriers because compartment-saturation paths are not known a priori. This differs from two-phase systems where saturation changes in one of two directions. A Haines jump in SG mode occurs when the saturation state in one compartment jumps across a barrier to another state in the neighbor energy valley, while compensating saturation changes occur in the other compartments to maintain the global, target gas saturation. For the transition of states during gas invasion marked by magenta and cyan circles in Figure 2, s g,1 increases at the expense of lower s o,1 and s w,1 in compartment 1, while in the other compartments s g,2 and s g,3 decrease at the expense of higher s o,2 , s w,2 , s o,3 , and s w,3 . The global saturation paths show that the overall result for this change of states is slightly higher S g and S w and lower S o . This emerges as a characteristic feature for Haines jumps in SG mode: Gas invasion (retraction) occurs with accompanying water invasions displacing oil, or oil invasions displacing water, leading to abrupt directional changes of the saturation path (that is, oil and water saturation jumps). Reversible displacements within the same energy valleys aim at driving the saturation path back on track so that equilibrium states in PG and SG modes coincide between Haines jumps. The main observation is that saturation paths in SG mode display a characteristic stepwise staircase shape, while in PG mode they consist of gas-saturation jumps intertwined with HELLAND ET AL. 10.1029/2021WR029560 8 of 16 Figure 1. Demonstration of the discrete-domain approach for a two-phase oil/water system using N = 3 and K = 4. Top, left: Capillary pressure curves for primary drainage and imbibition. Top, right: Evolution of local saturations for compartment 1 (blue), 2 (red), and 3 (green). The large circles show the compartment saturations for the corresponding three capillary pressures highlighted on the P ow (S w )-curves. Bottom, left: Gibbs free energy for compartment 1 (blue), 2 (red), and 3 (green), as functions of compartment water saturations for the three P ow shown by circles on the P ow (S w )-curves. The convex segments are the equilibrium locus (bold curves), the concave segments are energy barriers (dash-dotted curves), and the circles show the three metastable compartment states for the three P ow shown on the P ow (S w )-curves. Bottom, right: Phase-specific energies f o,i and f w,i used in the simulations for the three compartments, and total compartment energies (f o,i + f w,i , i = 1, 2, 3), as functions of compartment water saturations. The convex segments are the equilibrium locus (bold curves). more steady saturation changes. Differences in compartment-saturation evolution confirm the behavior: SG mode shows significant fluid redistribution among compartments and saturation fluctuations, while PG mode shows stepwise, monotonically increasing or decreasing compartment saturations.
HELLAND ET AL.
10.1029/2021WR029560 9 of 16 Figure 2. Demonstration of the discrete-domain approach with N = 3 and K = 4 for a three-phase system. Top row, left: Three-phase capillary pressure curves for a gas invasion-retraction hysteresis loop, using PG and SG displacement modes. Gas invasion occurs at S w = 0.4 after two-phase oil/water imbibition (Figure 1). Top row, right: Three-phase saturation paths, where the large circles show the saturations for the corresponding three states on the P αβ (1 − S g )curves. Middle row, left: Equilibrium locus of Gibbs free energy, obtained from Equation 7, for compartment 1 (blue), 2 (red), and 3 (green), as functions of compartment liquid saturations, for the three capillary pressures shown by large circles on the P αβ (1 − S g )-curves. Middle row, right: Phase-specific energy functions f α,i (s α,i ), α = g, o, w, i = 1, 2, 3, used in the simulations. Bottom row: Evolution of local saturations for compartment 1 (blue), 2 (red), and 3 (green) for the gas invasions with displacement modes SG (solid curves) and PG (dashed curves). The large circles show compartment saturations for the states highlighted on the P αβ (1 − S g )-curves.
The parameters ω α and c α,i in Equation 16 make it possible to explore hysteretic three-phase displacements at other realistic wetting states. A weakly oil-wet system (Hui & Blunt, 2000;van Dijke & Sorbie, 2002), where water is intermediate-wetting phase and gas nonwetting phase, requires ω o > ω w > ω g , c g,i > 0 > c o,i > c w,i , and c g,i > |c w,i |, i, …, N. In a strongly oil-wet system (Hui & Blunt, 2000;van Dijke & Sorbie, 2002), where gas is intermediate-wetting phase and water nonwetting phase, we must have ω o > ω g > ω w , c g,i > 0 > c o,i > c w,i , and |c o,i | < c g,i < |c w,i |, i, …, N. To explore three-phase hysteresis for these oil-wet states, we use a slightly larger system with N = 6 and K = 8. Further, ω o = 11, c w,i is evenly spaced between c w,1 = −0.15 and c w,6 = −0.01, and c o,i = −0.001, i = 1, …, 6. Then, for the weakly oil-wet state, ω w = 6, ω g = 1, and c g,i = −1.2c w,i , and for the strongly oil-wet state, ω w = 1, ω g = 6, and c g,i = −0.5c w,i , i = 1, …, 6. We add a water-wet case for comparison, with parameters ω w = 11, ω o = 6, ω g = 1, c w,i evenly spaced between c w,1 = 0.01 and c w,6 = 0.15, c o,i = 0.001, and Figure 3 presents two-and three-phase results from simulations of these wetting states, using Δp = 0.05 in PG mode and Δs = 0.0025 in SG mode. These steps are sufficiently small to capture all metastable states. The three-phase hysteresis loops begin at S w = 0.5 on the pressure-controlled two-phase imbibition curve. The capillary pressure curves capture the expected differences in pressure levels for the different wetting states, and they all show the same characteristic Haines jump behavior. In SG mode, the three-phase capillary pressure curves and saturation paths fluctuate due to Haines jumps. For intervals with reversible displacements both the capillary pressure curves and saturation paths coincide in PG and SG mode. The compartment-saturation profiles show fluctuations in SG mode, particularly for the water-wet state, which indicates fluid redistribution among compartments, while for PG mode, these curves are stepwise and monotonic curves. Further, the compartment-saturation profiles exhibit hysteresis and show that fluid distributions vary with wetting state.
Along with the hysteresis loops, we also simulate gas invasion in PG mode from S w = 0.3 and S w = 0.7 to investigate saturation-dependencies of the three-phase capillary pressures in the presence of hysteresis. Simple capillary-tube bundle models, where hysteresis is absent, predict that two of the capillary pressures are functions of one saturation, as follows (Hui & Blunt, 2000;van Dijke & Sorbie, 2002): P go (S g ) and P ow (S w ) (water-wet state), P ow (S o ) and P gw (S g ) (weakly oil-wet state), and P go (S o ) and P gw (S w ) (strongly oil-wet state). In each case, the remaining capillary pressure is a function of two saturations. However, our simulated results deviate from this idealized description, which we attribute to the inherent hysteresis of the method and chosen energy functions. Instead we find that all three capillary pressures depend on two saturations. For example, the simulated saturation paths at constant P ow in Figure 3 shows that P ow depends on two saturations although significant oil displacement (especially in the water-wet system) indicates strong S w -dependency. Further, the P go (1 − S g )-curves (water-wet state) and P gw (S w )-curves (strongly oil-wet state) for the three PG-based gas invasions do not coincide, implying that also these capillary pressures are functions of two saturations.
Finally, we explore three-phase PGSO and SGSO displacements in a large water-wet system with constant global oil saturation S o = 0.3, established after a two-phase pressure-controlled imbibition. We use N = 24, , i = 1, …, 24. To capture all metastable states of this large system, we use Δs = 0.001 (SGSO) and Δp = 0.05 (PGSO). Figure 4a shows P ow (S w )-curves for two-phase hysteresis loops, as well as results from the subsequent PGSO-and SGSO-simulations of main and nested three-phase hysteresis loops for gas invasions and gas retractions. For these displacement protocols all three capillary pressures vary during the displacement. Even in PGSO mode, P go and P ow display nonmonotonic variations, implying significant redistribution of the globally preserved oil phase. A sharp increase (or, drop) of P go leads to a corresponding drop (or, increase) of P ow for prescribed P gw . For gas invasion in PGSO mode P ow increases almost monotonically, while P go vary nonmonotonically around the same pressure level. Three-phase pore-scale simulations of the same process explain this observation (Helland et al., 2019): A double displacement in which gas displaces a disconnected oil phase through wide and narrow pore channels, that in turn pushes oil/water interfaces further into narrower pore spaces, leads to nonmonotonic P go and increasing P ow . The three-phase capillary pressure curves from SGSO simulations display a more irregular sawtooth structure than the two-phase curves due to significant oil redistribution among compartments. The s o,i (1 − S g )-curves show that this redistribution generally is different in PGSO and SGSO simulations so that the corresponding three-phase capillary pressure curves also deviate from each other.
For this system we also simulate main and nested hysteresis loops with PG and SG mode, starting at S w = 0.3 from the two-phase pressure-controlled imbibition. In this case, the capillary pressure curves for PG and SG modes are more aligned with each other and SG mode also displays a more regular sawtooth structure than SGSO mode, see Figure 4b. Because the saturation paths for gas invasion and gas retraction vary, nested hysteresis loops are not fully contained within the main hysteresis loop. This is in contrast to the PGSO and SGSO cases where the saturation path is enforced. The SG and SGSO simulations in Figure 4 display severe pressure fluctuations near the end of the displacement processes. This is an effect of excluding compartments where saturations have reached their limits s min or s max , which significantly reduces the number of possible metastable states in subsequent calculations. Still, SG simulations of this large system exhibit less significant capillary pressure fluctuations in most of the saturation interval, while the stepwise, staircase features of the saturation paths in SG mode are less prominent, compared with smaller systems (Figures 1-3). Thus, increasing the system size (i.e., N and K) yields smaller fluctuations because the density and number of metastable states increases, while hysteresis persists. This shows the potential of the discrete-domain approach as a method for upscaling three-phase capillary pressure curves with hysteresis from pore scale to larger scales where such fluctuations vanish.
Another feature which makes the discrete-domain approach attractive for upscaling is that the compartment interfacial free energies are extensive properties. Therefore, energy functions defined on a highly resolved compartmentalization of a porous medium (e.g., where each compartment describes only a few pores) can easily be constructed for coarser compartment resolutions. This allows simulations with the discrete-domain method on different scales. Coarsening of the energy functions could to some extent eliminate or moderate the energy valleys and barriers, resulting in smoother capillary pressure curves on coarse scales. Pore-scale simulations on segmented rock images (e.g., Helland et al., 2019Helland et al., , 2017Zacharoudiou & Boek, 2016) or flow experiments with microtomographic imaging of pore-scale fluid distributions on rock samples (e.g., Alhosani et al., 2019;Berg et al., 2013;Khishvand et al., 2016;Paustian et al., 2021;Scanziani et al., 2020) provide valuable data to validate and test such workflows for upscaling. The coarsening may result in a single compartment for the entire porous sample. This allows for utilizing the discrete-domain approach on larger cores by stitching together several such porous samples, where each sample corresponds to a compartment with a statistically representative energy function.

Conclusions
This work presents a discrete-domain approach to three-phase displacements with rate-independent hysteresis in porous media. Constrained energy minimization leads to evolution equations for compartment saturations that collectively describe a wide range of three-phase displacements, including reversible and irreversible pressure-and saturation-controlled displacements with or without preservation of one of the defending phases. Saturation-controlled invasion, which mimics flow controlled by low rate, shows that three-phase irreversible displacements occur with significant fluid redistribution among compartments. This leads to abrupt capillary pressure jumps as previously seen in two-phase systems, and in addition, abrupt saturation changes and fluctuating saturation paths, which we find is a new characteristic feature for three-phase displacements. Both pressure and saturation jumps decrease with increasing system size. Three-phase displacements with a preserved saturation exhibit the most substantial fluid redistributions, and in these cases both pressure-and saturation-controlled displacement generate nonmonotonic capillary pressure curves, consistent with quasi-static pore-scale simulations (Helland et al., 2019;Jettestuen et al., 2021). The presented discrete-domain approach is a fast and tractable method to investigate threephase hysteresis in porous media. As the method uses saturation evolution equations phase by phase, it can easily be adapted to an arbitrary number of phases. A current drawback is that the compartments are not HELLAND ET AL.
10.1029/2021WR029560 12 of 16 Figure 3. Application of the discrete-domain approach to water-wet (left column), weakly oil-wet (middle column), and strongly oil-wet (right column) states, using N = 6 and K = 8. First row: Two-phase oil/water capillary pressure curves for primary drainage and imbibition from pressure-and saturation-controlled simulations, representing the displacement history before gas invasion. Second row: Saturation paths for hysteresis loops of gas invasion (from S w = 0.5) and gas retraction from PG-and SG-simulations, as well as saturation paths from PG-simulations of gas invasions at S w = 0.3 (dashed curves) and S w = 0.7 (dash-dotted curves). Third row: Capillary pressure curves for the hysteresis loops and the gas invasions from S w = 0.3 (dashed curves, P go in green and P gw in black) and S w = 0.7 (dash-dotted curves, P go in green and P gw in black). Fourth, fifth, and sixth row: Evolution of local saturations in a compartment subset {1, 2, 3} during the hysteresis loops with SG (solid curves) and PG (dashed curves) displacement protocols. spatially correlated and hence they are all accessible for invasion. In future work we will use spatially dependent energy functions from pore-scale experiments and pore-scale simulations to validate the approach and to evaluate effects of wetting state on three-phase irreversible displacement. As compartment energies are extensive quantities, establishing a direct link between pore scale and the compartmental description will provide a reliable upscaling method. Finally, compartment energies tailored to match pore-scale data will reveal whether the energy functions themselves exhibit hysteresis and if they require other state variables (like interfacial area) to describe uniquely drainage and imbibition. Figure 4. Three-phase hysteresis behavior during alternate gas invasions and gas withdrawals for the case with N = 24 and K = 14. (a) Top, left: Two-phase oil/water capillary pressure curves for primary drainage and imbibition representing the displacement history before gas invasion. Top, right: Three-phase saturation paths from simulations with constant global oil saturation (PGSO and SGSO). Second row: Three-phase capillary pressure curves for the main hysteresis loop and a nested hysteresis loop from PGSO and SGSO simulations. Third row: Evolution of local saturations in a subset of compartments {1, 4, 10, 14} during the first gas invasion for displacement protocols SGSO (solid curves) and PGSO (dashed curves). (b) Saturation paths (left) and capillary pressure curves (right) from PG-and SG-simulations of two gas invasion/withdrawal cycles, demonstrating nested hysteresis-loop behavior for three-phase displacements. Initial global water saturations for the first gas invasions are S w = 0.7 (PGSO and SGSO) and S w = 0.3 (PG and SG).