Three‐Phase Fractional‐Flow Theory of Foam‐Oil Displacement in Porous Media With Multiple Steady States

Understanding the interplay of foam and nonaqueous phases in porous media is key to improving the design of foam for enhanced oil recovery and remediation of aquifers and soils. A widely used implicit‐texture foam model predicts phenomena analogous to cusp catastrophe theory: The surface describing foam apparent viscosity as a function of fractional flows folds backwards on itself. Thus, there are multiple steady states fitting the same injection condition J defined by the injected fractional flows. Numerical simulations suggest the stable injection state among multiple possible states but do not explain the reason. We address the issue of multiple steady states from the perspective of wave propagation, using three‐phase fractional‐flow theory. The wave‐curve method is applied to solve the two conservation equations for composition paths and wave speeds in 1‐D foam‐oil flow. There is a composition path from each possible injection state J to the initial state I satisfying the conservation equations. The stable displacement is the one with wave speeds (characteristic velocities) all positive along the path from J to I. In all cases presented, two of the paths feature negative wave velocity at J; such a solution does not correspond to the physical injection conditions. A stable displacement is achieved by either the upper, strong‐foam state, or lower, collapsed‐foam state but never the intermediate, unstable state. Which state makes the displacement depends on the initial state of a reservoir. The dependence of the choice of the displacing state on initial state is captured by a boundary curve.


Introduction
Catastrophe theory, initially founded by Thom in 1960s and further developed by Zeeman in 1970s, is a branch of bifurcation theory for dynamical systems and of singularity theory in geometry (Arnold et al., 1999;Thom & Zeeman, 1974;Wiggins, 2013;Zeeman, 1977). The cusp catastrophe, among the seven elementary catastrophes described by Zeeman (1977), is schematically illustrated in Figures 1 and 2. The theory of cusp catastrophe states that equilibrium behavior controlled by two independent quantities forms a smooth surface folding on itself. Theoretically, all states between the two edges of the fold in Figure 1 are unstable and cannot be observed in nature. At each edge of the fold, the system makes a sudden and dramatic jump between the two states, illustrated in Figure 2, upon a small change of the controls.
In principle, any event meeting the conditions of the theory would exhibit behavior similar to that in Figures 1 and 2. This is found in many models in physics and engineering, for example, in wave propagation, minimum surface area, nonlinear oscillations, and elasticity (Golubitsky & Keyfitz, 1980;Guckenheimer, 1986;Holmes & Rand, 1976;Kravtsov & Orlov, 1990, Purohit & Bhattacharya, 2003. Some theoretical predictions have been verified experimentally or in practice, with others not yet observed. Tang, Ansari, and Rossen (2019) find similar phenomena for foam flow with oil in porous media. Specifically, the widely used foam simulation model in STARS (Computer Modeling Group (Calgary, Alberta, Canada), 2015) predicts, for some combinations of oil, water, and gas fractional flow, three different sets of saturations. Tang, Vincent-Bonnieu, and Rossen (2019) show that this model gives realistic representation of steady state foam behavior with oil. Either this simulator or a foam model similar to that in the simulator has been used to represent a variety of coreflood studies and field applications of foam Chalbaud et al., 2002;Ma et al., 2013;Norris et al., 2014;Rognmo et al., 2018;Sharma et al., 2017;Spirov et al., 2012). Thus, it is essential to understand the model's behavior, including how it represents displacements with multiple possible displacing states. The study here presents physical insights into these phenomena and their implications for the dynamics of oil displacement by foam.
Numerous processes involve foam flow through porous media, for example, oil displacement in reservoirs (Rossen, 1996;Kovscek et al., 1997;Li et al., 2010), removal of nonaqueous phase liquid contaminants in aquifers and soils (Geistlinger et al., 2009;Jeong & Corapcioglu, 2005;Johnson et al., 2001;Kao et al., 2008) and CO 2 storage (Adebayo, 2018;Adebayo, 2019;Iglauer, Paluszny, et al., 2011;Iglauer, Wülling, et al., 2011;Juanes et al., 2006, Vitoonkijvanich et al., 2015. Foam is not a separate phase but a dispersion of gas in liquid such that gas bubbles are separated by interconnected liquid films, called lamellae. The applications above mainly rely on the fact that foam can reduce gas mobility considerably by trapping gas bubbles via those lamellae (Rossen, 1996). Most oils are detrimental to foam, and this affects significantly the effectiveness of foam for gas mobility control. The interaction between foam and oil is complex and is not yet fully understood (Farajzadeh et al., 2012).
Currently, there exist two groups of models describing foam dynamics in porous media: population-balance models (Kam & Rossen, 2003;Kovscek et al., 1995) and implicit-texture (IT) models (Cheng et al., 2000;Computer Modeling Group (Calgary, Alberta, Canada), 2015). The two groups of models capture different physics involving foam generation and destruction as well as foam behavior at steady state (Lotfollahi et al., 2016).
Population-balance foam models capture the dynamics of foam generation and destruction through a dynamic calculation of bubble density (number of bubbles per unit volume). In the absence of oil, the model of Kam and Rossen (2003) predicts behavior analogous to the cusp catastrophe of Figure 2, where pressure gradient ∇p is plotted as a function of gas (u g ) and water (u w ) superficial velocities (Afsharpoor et al., 2010;Kam & Rossen, 2003). This behavior has been experimentally confirmed in foam corefloods in two ways: first, by showing a sudden jump in ∇p upon a small increase in superficial velocity and, thereafter, a hysteresis with the velocity decreasing; second, in experiments with ∇p fixed, by revealing the entire S-shaped curve illustrated in Figure 2 (Gauglitz et al., 2002). Nevertheless, few current population-balance models represent the effect of oil on foam. The model of Myers and Radke (2000) accounts for the effect of oil by reducing bubble-generation rate in that oil occupies part of pore space and reduces the number of sites where lamellae could be created. This model does not capture the impact of oil on foam stability (e.g., oil condition for foam collapse). Ma et al. (2018) in a recent study  (left) Illustration of the cusp catastrophe: a smooth surface describes behavior at equilibrium as a function of two independent controls; adapted from Zeeman (1977).
attempt to represent the effect of oil by increasing bubble destruction rate in a population-balance model. In separate ways, that model would be expected to alter the low-and high-quality foam regimes, respectively, in a way similar to that of the model described below. The details of behavior of their model and the validity of the model to represent foam flow with oil is not yet experimentally justified.
IT models assume local equilibrium, meaning that foam everywhere immediately reaches a state where bubble generation rate matches destruction rate. Foam texture, that is, bubble size, is not represented explicitly in IT models but is reflected implicitly through a gas mobility-reduction factor. The IT foam model in the STARS simulator defines gas mobility in foam as a function of water (S w ) and oil (S o ) saturations (i.e., fraction of pore volume occupied by water or oil). The approximations of local equilibrium and implicit texture simplify the analysis of foam flow dynamics, in particular in a complex interaction with oil. The model and parameters used below are based on a fit of model parameters to foam behavior without oil (Cheng et al., 2000) and is consistent with foam behavior with oil reported by Tang, Vincent-Bonnieu, and Rossen (2019).
A modeling study by Tang, Vincent-Bonnieu, and Rossen (2019) using the STARS foam model suggests the behavior illustrated in Figures 1 and 2, for steady state foam flow with oil. The behavior was evaluated through foam apparent viscosity μ app , that is, the inverse of total relative mobility (see Foam Model section below for its specific definition). With oil, μ app predicted by the model, when plotted on a ternary diagram of oil, water, and gas fractional flows, appears as a surface folding on itself. Figure 3 illustrates this behavior in terms of ∇p as a function of fractional flows, with more details given in the Foam Model section in terms of μ app . The folded region means that, as in Figure 3, there are three possible steady states (associated with different saturations), which fit the same injection condition J defined by a given set of fractional flows. The middle state is intrinsically unstable and therefore not seen in nature. The existence of multiple steady states has not been directly confirmed in the laboratory, but it is physically plausible and consistent with observations of foam behavior (Tang et al., 2018). At the same injection rates, a strong foam, with large μ app and ∇p, might displace oil, resulting in a low oil saturation and thus stable foam; nevertheless, a collapsed foam state with low μ app and ∇p might leave oil saturation high and foam unstable. This raises an essential question concerning the effectiveness of foam displacements: Which of the multiple steady states, with different μ app and ∇p, all fitting to a same set of injected fractional flows, actually occurs in a given displacement?
Tang , Ansari, and Rossen (2019) performed simple 1-D numerical simulations with fixed injection rates corresponding to cases with multiple possible states. Their simulation results suggest a displacement by either the upper state, with large μ app , or the lower state, with small μ app , but never the (unstable) one in the . Pressure gradient ∇p as a function of water (f w ) and gas (f g ) fractional flows at total superficial velocity u t = 3.53 × 10 −5 m/s, predicted by the wet-foam representation in Appendix A with a fixed ratio of (f o /f w ) = 0.25; adapted from Tang, Ansari, and Rossen (2019). The trace follows a sequence: a-e. The dashed line indicates a case of multiple steady states fitting same set of fractional flows, caused by the portion of the curve between c and d folding toward the lower right corner. Model parameters used for the illustration are referred to Table A1 in Appendix A. middle. These observations are consistent with catastrophe theory. However, a numerical simulation does not explain why a particular (strong or weak) foam is seen in the displacement.
We here present an analysis of a displacement with multiple possible injection states fitting the same injected fractional flows, using an analytical approach, the wave curve method (WCM) for three-phase flow. The WCM is also referred to, in petroleum engineering, as fractional-flow theory or the method of characteristics. Fractional-flow theory excludes the numerical artifacts that can afflict foam simulation (Rossen, 2013).
A computer-assisted design package, RPn (n-dimensional Riemann Problem) that applies the WCM has been developed in the group of fluid dynamics at IMPA (Instituto Nacional de Matemática Pura e Aplicada) in Brazil (Azevedo et al., 2010;. This program efficiently determines Riemann solutions for 1-D three-phase flow, including the complications of foam. The Riemann solutions provide a mathematical criterion that distinguishes states observed or not observed in the displacement. Using this tool, we explore the dependence of the displacement on initial condition.
The significance of our findings for field applications is discussed below. Some suggestions are provided for the direct experimental verification of multiple steady states through foam corefloods in laboratory.

Three-Phase Fractional-Flow Theory
The core of fractional-flow theory concerns modeling transport of fluids in porous media, specifically, by solving for phase saturations as a function of position x and time t (Lake et al., 2014). The solutions are usually illustrated as a saturation profile as a function of position at a given time or saturation history at a given position.
In a system of foam flow with oil, the presence of three phases means there are two independent variables to be determined, water (S w ) and oil (S o ) saturations as functions of (x, t). Gas saturation is For the purpose of our study, we have made the following assumptions: 1. 1-D flow in a homogeneous porous medium. The model does not represent fingering that would occur in 2-D or 3-D, but it does identify regions or fronts with viscous instability where fingering would be likely. 2. incompressible fluids and rock, 3. immediate attainment of local-equilibrium behavior, 4. no dispersive processes (diffusion, dispersion, or flow driven by capillary-pressure gradients), 5. gravity can be ignored, 6. no mass exchange between phases (immiscible phases), 7. isothermal process, 8. Newtonian rheology for all phases, and 9. uniform concentration of surfactant in the aqueous phase throughout the medium, which implies that adsorption has already been satisfied.
The local-equilibrium assumption means that foam generation and destruction rates come to equilibrium instantaneously. In local-equilibrium modeling of a displacement, this assumption means foam at any location is modeled as a function of local phase saturations. This approximation is justified by foam-displacement experiments without oil (e.g., Kovscek et al., 2010) and by simulations in which such effects are explicitly introduced (Kam et al., 2007).
Capillary diffusion (flow driven by local gradients of capillary pressure) is modeled as a function of saturations, which affects phase distributions in a displacement. Multiple steady states as represented by the model in Appendix A are a result of combined effects of S w and S o . With or without localized capillary diffusion, in a sufficiently long displacement, there are some combinations of S w and S o in ternary saturation space that correspond to different foam viscosities but the same fractional flows. The existence of the multiple steady states is, therefore, independent of capillary diffusion.
With foam, gas mobility can be reduced by a factor, for example, of thousands (Tang, Vincent-Bonnieu, & Rossen, 2019). Compared to viscous forces in foam injection, gravitational forces are negligibly small. That is supported by CT foam corefloods where saturations are measured (e.g., Tang, 2019).
The modeling of foam-oil flow described below assumes a homogeneous porous medium. Some other physical factors not considered here include formation heterogeneity, surfactant adsorption, gradients of capillary pressure, and gravity separation. These would greatly complicate conditions in geological formations for foam generation and/or propagation. For instance, Shah et al. (2018) experimentally show that low-to-high permeability boundaries assist in foam generation through capillary snap-off. Al Ayesh et al. (2017) in a modeling analysis show that permeability affects foam stability and strength, thus resulting in diversion of flow from high to low permeability layers. Our goal is first to understand behavior in a simplified system before extending the analysis to these complications.
The system is governed by two mass-conservation equations for x > 0 and t > 0: where φ is porosity, a rock property (a volume fraction of rock that is pore space) and u is the total superficial velocity of the three phases, that is, u ≡ u w +u o +u g . f w and f o in equations (1) and (2) are the fractional flows of water and oil, defined as the fraction of phase volumetric flux to the total volumetric flux: where subscript j = w, o, or g represents water, oil, or gas, respectively.
Superficial velocity u j of phase j is determined by Darcy's law: where k is the absolute permeability of the porous medium, μ w , μ o , and μ g are the water, oil, and gas viscosities, |∇p| is the absolute magnitude of pressure gradient, and k rj is the relative permeability of phase j. The ratio of (k rj /μ j ) is referred to as the relative mobility of phase j.
Substituting equation (4) into equation (3) yields where superscript f in k f rg indicates the presence of foam. We assume that presence of foam does not alter the relative-permeability function for water or oil but only that for gas. This is supported by a number of experimental observations and greatly simplifies the physics and modeling of foam flow in porous media (Rossen, 1996;Schramm, 1994).
In the STARS foam model (see Appendix A for specific algorithms), k f rg with foam is defined as foam-free gas relative permeability, k rg reduced by a mobility-reduction factor FM (Computer Modeling Group (Calgary, Alberta, Canada), 2015). The scaling factor FM, given in equation (A2) in Appendix A, is a function of a series of physical factors affecting foam stability and degree of gas-mobility reduction by foam. Two major factors in FM concerning water and oil saturations are included in this study, meaning that FM is a function only of saturations.
The commonly used models for k rj include Stone I and II (Stone, 1970(Stone, , 1973 and the Corey model. In a case of three-phase flow, k rj of the intermediate-wetting phase in Stone I and II is related to that of the other two phases based on channel-flow theory, rather than a unique function of its own saturation. Foam alters gas mobility so drastically that Stone's models can give unphysical results with foam, for example, positive relative permeability at zero oil saturation. For simplicity, therefore, we use a Corey-type relative-permeability model, where k rj of phase j without foam is a function of S j alone:

Water Resources Research
TANG ET AL. 10,323 where k 0 rj is the endpoint relative permeability, n j is the Corey exponent that reflects wettability, and S wc , S or, and S gr represent the residual saturations of each phase. Equation (A2) in Appendix A for FM, with F 2 and F 3 included, combined with equation (6) for k rj , suggest that fractional flow f j in equation (5) is a function only of saturations.
To simplify equations (1) and (2), dimensionless position, x D , and time, t D are introduced: where L is the 1-D reservoir length and t D is the number of movable pore volumes injected, with movable pore volume scaled by (1-S wc -S gr -S or ). Equations (1) and (2) then are simplified to where S and F are both vectors, that is, Fractional-flow theory interprets a displacement process in terms of wave propagation (Avraam & Payatakes, 1995;Buckley & Leverett, 1942;Charbeneau, 1988;Lake et al., 2014;Reynolds & Krevor, 2015;Wooding & Morel-Seytoux, 1976). In principle, all saturations between J and I along a composition path exist at the origin at t D = 0. Upon injection, water, oil, and gas propagate starting from the origin with given wave speeds as a function of saturations. Solving for the saturation profile and history is then equivalent to solving for waves speeds, η(S): Substituting η(S) into equation (9) and using chain rule of differentiation, the system is rearranged to (Lake et al., 2014): where I * denotes the 2-by-2 identity matrix, with an asterisk to distinguish it from the initial state, I. J(S) is the Jacobian matrix: where J ij (S) = ∂f i /∂S j with i and j = w, o, and g denoting water, oil, and gas, respectively.
Mathematically, η(S) to be solved is the eigenvalue of the Jacobian matrix J(S). A physical problem of solving for wave speeds, η(S), is eventually converted to a mathematical problem concerning eigenvalues of the Jacobian matrix in equation (12) (Castañeda, 2018;Lax, 1957;Liu, 1974).

Wave-Curve Method
The WCM implemented in the RPn solves the system described by equation (9) for two major outputs: a composition path from J to I that provides S along the path and the associated wave speed for each value of S. The two types of solutions together define the structure of a displacement. From the saturations and wave speeds arising from x D = t D = 0, the saturation profile can be determined at any time and the saturation history at any location. A number of studies of multiphase flow in porous media, concerning a hyperbolic system of conservation laws, provide a detailed description of the method (Azevedo et al., 2010;Castañeda, 2018;. This section briefly explains the general principles of the WCM applied in the RPn program (Liu, 1974;Smoller, 2012), especially to distinguish the Riemann problem that the WCM solves from that of physical displacements.
Generally, the composition path for J displacing I is determined by constructing two families of wave curves via the WCM, as illustrated in Figure 4: a forward slow-wave curve and a backward fast-wave curve. The forward slow wave curve starts from J; the smaller eigenvalue of J (S) in equation (12) at each saturation is assigned in the WCM as the characteristic speed at saturation S. The corresponding eigenvector gives the direction of saturation change from the current saturation to the next saturation. At this saturation, the eigenvalues of J(S) in equation (12) are again calculated and the process continues until the entire slow path is determined.
Saturations within a shock, if it occurs, are unphysical and do not appear on a saturation profile. That leads to a discontinuity in saturations, meaning that the shock speed cannot be resolved through J(S) comprising derivatives that are a function of saturations. The shock and its speed σ that reach I are then determined by constructing a Rankine-Hugoniot (RH) locus via the RH condition (Azevedo et al., 2010;: where S I represent phase saturations at I. The RH condition is derived from a mass balance across a shock. This condition as in equation (13) defines all possible states that may reach I via a shock, seen from the dashed line in Figure 4, referred to as the backward fast-wave curve.
The two families of wave curves as illustrated in Figure 4 cross at some point denoted IJ, at which the slow speed is less than the fast wave speed. The spreading wave velocities of saturations along the forward slow-wave curve increase in the forward direction starting from J. The shock velocities of points, along the backward fast-wave curve, decrease in the backward direction initiating from I to IJ. The whole path, from J to IJ and then to I, is thus guaranteed to meet the compatibility of monotonically increasing wave velocity. This result hinges on the nonlinearity of relative-permeability functions. In the simplified case of linear relative-permeability functions, shocks (equation (13)) occur along the same curve as spreading waves (equation (11)), but this is not the case in general (Lake et al., 2014;Namdar Zanganeh et al., 2011).
To avoid confusion in mathematical and engineering terminologies, we clarify that all the states that connect J to I (solved through the two families of wave curves), as shown in Figure 4 and subsequent illustrations, are referred to as a composition path. Saturation points on a spreading-wave curve are physical, which can be observed in a displacement. The two endpoints of a shock wave (i.e., intersection state IJ and initial state I in Figure 4) are physical and can also be observed, whereas the other saturation points between these two states within the shock are unphysical, which cannot be observed.
Note that the WCM for three-phase flow solves a problem with different constrains than a coreflood with specified injection rates of phases. In a coreflood, at t D = 0, initial state I is defined as follows: Upon injection, initial state I, present for 0 < x D ≤ 1, is displaced forward by injection state J starting at x D = 0 with specified fractional flows of phases, as schematically illustrated in Figure 5. Strictly, the WCM does not specifically solve for a physical displacement of I for 0 < x D ≤ 1 by J at x D = 0. Instead, as shown in Figure 6, it solves for an initial state with I present for x D > 0 and J for x D ≤ 0: where S is a vector of saturations in a case of three-phase flow. Starting at time t D = 0, injection of state J (from x D < < 0) begins, and the state evolves. The analysis below shows that this distinction is crucial to determining the correct displacing state J among multiple possible injection states. In particular, only a composition path with all positive wave velocities can represent a physical displacement by J at x D = 0.

Foam Model
A key to the success of foam-related applications is the physical stability and strength of foam with oil. We here apply the widely used IT foam model in the STARS simulator (Computer Modeling Group (Calgary, Alberta, Canada), 2015; Cheng et al., 2000). It includes two algorithms describing the effect of oil on foam properties: the wet-foam representation and the dry-out representation. Appendix A gives the wet-foam representation that describes impacts of a series of physical factors on foam properties, that is, surfactant concentration, water saturation, oil saturation, salinity, shear-thinning rheology, and capillary number. Two major functions, in equations (A3) and (A4), are included in this study to quantify the effect of S w and S o on foam properties. The dry-out representation (equation (A3)) can also predict multiple steady states (Tang, Ansari, & Rossen, 2019), but we do not employ it here.
The model parameters we used are based on a detailed fit of coreflood data for foam without oil (Alvarez et al., 2001;Cheng et al., 2000) and are consistent with coreflood data on foam with oil (Tang, Vincent-Bonnieu, & Rossen, 2019).

Foam Representation on Ternary Diagram
In corefloods, the mobility reduction in gas by foam is measured through pressure drop or pressure gradient. In modeling, this is evaluated via a mobility-reduction factor FM (equation (A2) in Appendix A) as a function of S w and S o . Through the factor FM, foam modifies gas mobility by reducing gas relative  permeability. A larger value of (1/FM) indicates greater reduction to gas mobility and thus stronger foam. The model parameters used as given in Table A1 in Appendix A indicate a very strong foam, with an abrupt collapse of foam at S w near the limiting water saturation, fmdry (physically denoted as S * w ). Figure 7 plots (1/FM) in equation (A2) in ternary saturation space as defined in the wet-foam representation, where function F 2 in equation (A3), for the effect of S w , and F 3 in equation (A4), for the effect of S o , are considered. Generally, the factor (1/FM) splits the ternary saturation space into three regions: a full-strength foam region, a partially destabilized foam region, and a no-foam region. The region with foam either at full strength or partially destabilized has (1/FM) > 1; this region resides at the lower-left corner of the ternary diagram. The remaining, white region, with (1/FM)~1, indicates absence of foam.
The absence of foam indicated by the white region in Figure 7 arises mainly from two reasons related to foam stability and strength. Foam stability is controlled by the limiting capillary pressure, which corresponds to a limiting water saturation S w * (Khatib et al., 1988;Rossen & Zhou, 1995).
Along a direction parallel to the G-O binary, for S w lower than S w * , foam collapses. The abruptness of collapse depends on an adjustable parameter epdry in equation (A3), the value of which is large in this study as justified in coreflood measurements (Kim et al., 2005;Boeije & Rossen, 2015;Tang et al., 2018). Specifically, for S w just a bit less than S w * , F 2 in equation (A3) is nearly zero, setting the inverse of factor FM in equation (A1) close to unity. There is nearly no reduction in k f rg due to absence of foam in this range of S w . As S w rises, a transition zone in the interval around S w * , seen in Figure 7, is marked by an abrupt increase in (1/FM). This zone corresponds to an abrupt drop in k f rg around S w * in Figure A1 in Appendix A. For S w greater than S w * , foam is at full strength in the absence of oil.
In the other direction, parallel to the G-W binary, increasing S o weakens foam via equation (A4), leading to a larger value of FM and thus less gas-mobility reduction. The effect of S o in equation (A4) is bounded by two limits, an upper-limiting oil saturation fmoil and lower-limiting oil saturation floil. Specifically, for S o < floil, oil-saturation-dependent function F 3 in equation (A4) equals unity, meaning that oil has no detrimental impact on foam. For S o between floil and fmoil, oil shows a nonlinear effect, indicated by the color gradient in Figure 7. When S o > fmoil, F 3 is zero, with (1/FM) = 1; foam is killed completely.
This corresponds to the sudden rise of k f rg around S w , roughly 0.45 along the bottom axis of Figure A1 in Appendix A. The full-strength foam results from combined effects of S w and S o , specifically for S w > S w* and S o < floil.

Multiple Steady States in Foam Model
Experimental observations show that most oils destabilize foam (Farajzadeh et al., 2012;Rossen, 1996). The widely used IT foam model shown in Appendix A suggests a cusp catastrophe as in Figures 1 and 2, leading to multiple steady states. Tang, Ansari, and Rossen (2019) provide a detailed analysis of the occurrence of the multiple steady states. These phenomena are briefly described here for the purpose of our study.
The foam properties as shown in Figure 7 are usually evaluated through foam apparent viscosity μ app . A larger value of μ app indicates stronger foam. Treating multiphase flow in foam flow with oil as a pseudo-single phase and applying Darcy's law gives the definition of apparent viscosity, μ app : where u represents the total superficial velocity. With u = u w +u o +u g , solving for ∇p for three-phase flow (equation (4)) and substituting into equation (16) returns μ app as a function of (S w , S o ):  Table A1 in Appendix A.

Water Resources Research
TANG ET AL. 10,327 where the denominator (k rw =μ w þ k ro =μ o þ k f rg =μ g ) is the total relative mobility. Since fractional flows f j for each phase defined in equation (5) are also functions only of (S w , S o ), μ app is then correlated to f g , f o , and f w through (S w , S o ). Figure 8 plots μ app as a function of f g , f o , and f w , using the wet-foam representation in Appendix A. On a ternary diagram of fractional flows, the shape of foam apparent viscosity μ app forms a surface. The curves on the surface are plotted for μ app at various fixed ratios (f o /f w ), which is equivalent to fixing (u o /u w ). The blue curve, along the f g -f w binary, shows the changing trend of μ app without oil. As a group, the curves illustrate the shape of the surface.
A key feature illustrated in Figure 8 is that the surface folds back toward f g = 1 in the middle, due to the destabilizing effect of oil on foam, as explained by Tang, Ansari, and Rossen (2019). In the folded region, a set of fractional flows f g , f o , and f w corresponds to three possible steady states: a strong foam state with large μ app on the upper surface, a collapsed foam state with very small μ app at the bottom, nearly flat surface, and the intermediate state on the folding surface. The projection of the folded region onto the ternary diagram gives the dashed line, beyond which foam does not exist, either due to low S w or high S o . As briefly mentioned in the section 1, the key issue we answer in the following section is which of the multiple possible steady states is the displacing state to oil for injection with fractional flows in the folded region. The displacing state controls the success of a foam process, since different states have very different mobilities.

Displacing State Among Multiple Steady States
For a number of engineering applications described by a cusp catastrophe, the final state achieved in the folded region depends on the initial state (Golubitsky & Keyfitz, 1980;Guckenheimer, 1986;Holmes & Rand, 1976;Kravtsov & Orlov, 1990). For a given set of injected fractional flows, we define two types of initial conditions I, that is, unfavorable (I 1 ) or favorable (I 2 ) to displacement by stable foam. Figure 9 gives an example of a set of multiple steady states J 1 , J 2 , and J 3 , all fitting the same injected fractional flows, and two types 10,328 of initial states I 1 and I 2 . There are thus six composition paths, from each of J 1 , J 2 , and J 3 , to each of I 1 , and I 2 , respectively. The question is which path represents the solution for the given injected fractional flows and each initial condition. Table 1 summarizes the saturations of each possible injection state, J 1 , J 2 , and J 3 , and of the two types of initial conditions, I 1 and I 2 in Figure 9. All the parameters needed in the construction of composition paths are taken from Table A1 in Appendix A.
Note that I 1 and I 2 in Figure 9 were selected outside and inside the foam region, to illustrate the impact of the initial state on the choice of the displacing state among J 1 and J 2 and J 3 . However, the fundamental division of I in ternary saturation space, for the choice of the displacing state among the given set of J 1 , J 2 , and J 3 , is not necessarily the boundary of the foam region. The dependence of the displacing state on the initial state is captured through a boundary curve shown below.

Initial State Outside Foam Region
The two major outputs in Riemann solutions solved by the WCMcomposition paths and wave speeds-provide physical insights that unravel the occurrence and features of a displacement, that is, wave type, configuration, and propagation. We analyze Scenario 1 in Table 1 first, with the initial state I 1 outside the foam region, based on the Riemann solutions of cases 1, 2, and 3 with J 1 , J 2 , and J 3 as the injection state, respectively. Figure 10 shows the three composition paths from each of J 1 , J 2 , and J 3 to I 1 . Though J 2 might be intrinsically unstable as suggested by catastrophe theory, one could still theoretically construct a composition path from J 2 to I 1 . Tracks of the three paths from J to I on Figure 10 follow, respectively: Each path comprises two groups of wave curves, that is, the forward slow-wave group from J and backward fast-wave group from I as schematically illustrated in Figure 4. These two groups of wave curves in a path are distinguished by the intersection state IJ, that corresponds to the point C 1 in J 1 path, C 2 in J 2 path, and C 3 in J 3 path, respectively. Along each path in Figure 10, the portion from J to IJ represents the slow-wave group, and the rest portion from IJ to I corresponds to the fast-wave group.
The specific structure of the three paths in Figure 10 is indicated by dashed and solid curves that represent the shock and spreading waves, respectively. The labelled letters A k , B k , C k , and D k (k = 1, 2, and 3) indicate the transitions along each path from a shock to a spreading wave or the reverse. The paths from J 1 and J 2 to I 1 exhibit similar wave type and configuration. Both paths start from the injection state with a shock wave (from J 1 to A 1 along the path from J 1 or J 2 to A 2 along the path from J 2 ) followed by a spreading wave (from A 1 to B 1 or A 2 to B 2 ). Thereafter, there occurs a second shock (from B 1 to C 1 or B 2 to C 2 ), that is connected to a second spreading wave (from C 1 to D 1 or C 2 to D 2 ) which eventually reaches I 1 with a shock, that is, from D 1 Figure 9. Illustration of a displacement with multiple possible displacing states represented by J 1 (upper state-low mobility), J 2 (middle state), and J 3 (lower state-high mobility) corresponding to the same injected fractional flows. I 1 and I 2 represent two initial conditions, unfavorable and favorable for displacement by stable foam (upper state), respectively. The colored patch indicates the foam region. Saturations are given in Table 1.

Water Resources Research
to I 1 or D 2 to I 1 . Both paths feature a sharp inflection across the foam boundary at S o = fmoil. This arises from foam collapse that yields an abrupt rise in gas mobility. Capturing this sharp inflection in calculations requires very fine resolution in saturation steps. The path from J 3 (with very small gas fractional flow as given in Table 1), at which foam is fully collapsed, starts with a shock from J 3 to C 3 . A nearly invisible spreading wave in Figure 10 connects C 3 to D 3 , which shocks to I 1 . Figures 11a-11c illustrate trajectories of the associated wave speeds as a function of saturations along each path from J 1 , J 2 , and J 3 to I 1 in dimensionless x D versus t D space. Fractional-flow theory states that, as implied by equation (10) (Lake et al., 2014), the slope of each trajectory line in Figure 11 is the constant wave speed of a given feature (shock, characteristic, etc.). A red trajectory corresponds to a shock wave, as denoted by a dashed curve in Figure 10. Black lines in Figure 11 represent characteristics within spreading waves. I, J, and IJ mark constant-state regions, with IJ being the intersection state at C 1 , C 2, or C 3 as illustrated in Figure 4. The existence of the constant-state region at IJ is a major difference between two-and three-phase flows in porous media. Figures 11a-11c shows that the paths from each of J 1 and J 2 to I 1 pass through the bottom quadrant, reflecting negative wave speeds. In a physical displacement, the injected fractional flows F(J) are maintained fixed at x D = 0 as shown in Figure 5. However, for both paths starting from J 1 and J 2 as in Figures 11a and 11b, the fractional flows at x D = 0 with t D > 0 deviate from those at J 1 or J 2 , due to the negative wave velocities. The appearance of negative wave velocities then rules out a foam state being the displacing state in a physically acceptable displacement. Figure 10. Three composition paths in ternary saturation space constructed from each of J 1 , J 2 , and J 3 , respectively, to I 1 , which resides outside the foam region. The respective tracks of the three paths are denoted by:

Comparing the wave velocities in
The three saturation points at J 1 , J 2 , and J 3 represent three foam states, all fitting the same (f w , f o , f g ). A solid curve represents a spreading wave and a dashed curve represents a shock. Figure 11. Trajectories of the associated wave speeds as a function of saturations along each path from J 1 , J 2 , and J 3 to I 1 in Figure 10 in dimensionless x D versu t D space. A red trajectory marks a shock, with black lines representing characteristics within spreading waves. Only the path from J 3 gives wave speeds that are all positive.
Only the path from J 3 in Figure 11c yields wave speeds that are all positive; state J 3 therefore corresponds to the physical injection condition in Figure 5. J 3 , among the multiple possible steady states, is thus the displacing state in Scenario 1 that displaces initial state I 1 outside the foam region.
To illustrate the physical meaning of negative velocities in transport, we illustrate the Riemann solutions in terms of saturation profiles using the path from J 1 in Figure 10. Figure 12 displays the associated wave speeds from Figure 11a as a function of saturations (on the top axis) and also saturations as a function of position at fixed time (t D = 0.4) on the bottom axis. Upon injection, those saturations with positive velocities in Figure 12 move forward from x D = 0, whereas those with negative velocities propagate backward. For t D > 0, the state at x D = 0 is not J 1 ; thus, the Riemann solution for a displacement with negative wave velocity is physically unacceptable.

Initial State Inside Foam Region
A similar stability analysis is performed for Scenario 2 for the same set of foam states, J 1 , J 2 , and J 3 but with a different initial state I 2 now inside the foam region. Figure 13 illustrates the three composition paths from each of J 1 , J 2 , and J 3 to I 2 , which are solved again based on the definition of equation (15). The structures of the three paths are indicated by their associated tracks: J 1 -M 1 -I 2 ; J 2 -M 2 -I 2 ; J 3 -M 3 -I 2 . Specifically, the intersection state IJ that splits the slow-and fast-wave groups in the paths corresponds to the points M 1 , M 2 , and M 3 on Figure 13; these three points are located very close to each other, as seen in the left expanded view of Figure 13, and all belong to the Hugoniot locus of I 2 . The paths from J 1 and J 3 to I 2 follow a similar structure, both starting with a slow shock (from J 1 to M 1 or J 3 to M 3 ) and then reaching I 2 with a fast shock (from M 1 to I 2 or M 3 to I 2 ). The path from J 2 follows first a spreading wave from J 2 to M 2 and thereafter a shock from M 2 to I 2 . Figures 14a-14c show trajectories of the associated wave speeds as a function of saturations along each path from J 1 , J 2 , and J 3 to I 2 in dimensionless x D versus t D space. Only the path from J 1 to I 2 in the three paths has exclusively trajectories residing in upper quadrant, that is, only positive wave speeds. Therefore, only the strong-foam state, J 1 , among the multiple possible injection states, is the physically true displacing state to the initial state I 2 .
J 2 is not the displacing state in either scenario. The choice of the displacing state, among the multiple possible injection states, shows a dependence on initial state. In section 4.3, we show the nature of the dependence of the occurrence of J 1 or J 3 on initial state.  . Three composition paths in ternary saturation space constructed from each of J 1 , J 2 , and J 3 , respectively, to I 2 located inside the foam region. The track of each path is: The phase saturations at J 1 , J 2 , and J 3 are the same as those in Figure 10, fitting the same (f w , f o , f g ). A solid curve represents a spreading wave and a dashed curve marks a shock. Figure 14. Trajectories of the associated wave speed as a function of saturations along each path from J 1 , J 2 , and J 3 to I 2 in Figure 13 in dimensionless x D versus t D space. A red trajectory marks a shock, with black lines representing characteristics within spreading waves. Only the path from J 1 gives only positive wave speeds.

Boundary Curve for the Dependence of the Nature of the Displacement on I
To determine which of J 1 or J 3 is the displacing state to any initial state I in ternary saturation space, it is tedious and impractical to go through the calculations as in Figures 11 and 14 for every choice of I. It is then necessary to capture the universal dependence of the choice of the displacing state on initial state in whole ternary saturation space. This is especially crucial to improving the prediction and control of foam displacement with a given initial state in a reservoir, in particular with J corresponding to multiple possible injection states.
Intuitively, for an initial state I 2 inside the foam region, with low S o favorable for stable foam, it makes sense that strong-foam state J 1 , with large μ app and ∇p, would make the displacement. By the same logic, an initial state I 1 outside the foam region, where high S o kills foam, would be displaced by the collapsed-foam state J 3 with low μ app and ∇p.
However, the location of I inside or outside the foam region is not decisive for the displacement by J 1 or J 3 . Instead, ternary saturation space is divided by a "boundary curve" (defined through the intermediate state J 2 ) that determines which of J 1 or J 3 makes the physically acceptable displacement. Mathematically, the boundary curve is developed in terms of a forward fast-wave curve starting from J 2 . Specifically, in the cases presented in this study, its definition is given by the following two conditions: where σ∈R þ 0 , S is a vector of saturations, and subscripts f and s denote the fast and slow eigenvalues of Jacobian matrix in equation (12), respectively. Equation (18) is the expanded form of the Rankine-Hugoniot condition in equation (13). This condition gives all states S along the forward Rankine-Hugoniot locus starting from J 2 . Equation (19), as stated in Lax theory (Lax, 1957), defines the 2-Lax shocks starting from J 2 . λ f (J 2 ) in equation (19) is positive, guaranteeing that all admissible 2-Lax shocks from J 2 have positive wave speeds σ. These two conditions therefore in equations (18) and (19) represent a collection of admissible local 2-Lax shocks starting from J 2 , illustrated as the boundary curve in Figure 15.
Ternary saturation space as shown in Figure 15 is split into two regions: the region above the boundary curve, where state J 1 resides, and the region below the boundary curve, where state J 3 is located. For any I in the  (18) and (19) that captures the dependence of the nature of the displacement on initial state I.

Water Resources Research
TANG ET AL. 10,333 upper-left region, the strong foam state J 1 is the displacing state that gives a path with wave speeds all positive. In contrast, for any I in the lower-right region, the collapsed-foam state J 3 makes the physically acceptable displacement.
For a state I exactly on the boundary curve in Figure 15, one could not distinguish which of J 1 , J 2 , or J 3 makes the displacement, since all paths from J 1 , J 2 , and J 3 are physically admissible in theory. The path from J 2 is connected to I by an admissible 2-Lax shock. Given that J 1 , J 2 , and J 3 all fit the same injected fractional flows, the path from J 1 or J 3 would jump first to J 2 with zero velocity and then follow the same track as the path from J 2 to I . However, any perturbation that takes I off the boundary curve would always result in negative velocities along any sort of path from J 2 , which is physically unacceptable. A physically acceptable displacement would then shift to either J 1 or J 3 , depending on the location of I with respect to the boundary curve.

Significance for Field Applications
Foam for enhanced oil recovery is of course never co-injected with oil. Nevertheless, geological formations and fluids have many complexities (Lake et al., 2014): geological heterogeneity, fractures, unfavourable mobility ratios of displacing to displaced fluids, gravity segregation, etc. Part of the oil may remain in place due to limited displacement and sweep efficiency after initial gas injection. Since foam enhances greatly the sweep and displacement efficiency of gas injection, it is likely that injected foam flows with oil starting some distance from the well. Direct applications of foam to removing nonaqueous phase liquid contaminants in aquifers and soils would also involve cocurrent foam flow with fluids that may affect foam stability. All these situations may have multiple possible steady states fitting the same injected fractional flows for most process designs. If a displacing phase is a collapsed-foam state, gas-mobility control is much less effective than desired with foam. Negative results include viscous instability and slow displacement and production of oil. The theory and findings presented in this study may assist in predicting the displacing state for given initial conditions. One can then optimize the design of foam processes to ensure their success in engineering applications.
In a simplified one-dimensional model for foam processes in a homogeneous medium, the injected foam lies along the water-gas binary of the phase diagram and the initial state below and to the right. If gas has not previously been injected, the initial state is on the water-oil binary. The foam performance with the particular model parameters we implemented is not sufficiently effective for such a process. Only an initial state I in the small upper-left region in Figure 15 would be displaced by the strong foam state J 1 . The poor performance arises in part from the value of the oil parameter floil selected, close to S or . This implies that most physically realizable oil saturations destabilize foam somewhat. In addition, in the cases examined foam collapses completely at a modest oil saturation above fmoil.
In contrast, experimental observations and field pilots demonstrate that foam can show very good performance in displacing oil Tang et al., 2018). Thus, the parameters used here for illustration do not by any means represent all foam processes with all oils. By selecting surfactants that enhance foam tolerance to oil, that is, with floil well above S or and a large value of fmoil, a foam process more resistant to oil could be modeled by techniques similar to those here. Further efforts are needed to explore the displacement behavior with greater tolerance of foam to oil, so that one can represent cases of oil displacement by stronger foam and maximize the benefit of foam applications.

Experimental Verification of Multiple Steady States
Here we suggest two ways to confirm the multiple steady states predicted for foam flow with oil using laboratory corefloods. First, one could do a displacement, under different initial conditions, with an injection condition (fixed fluxes) that corresponds to multiple steady states. The core would be initialized either at a high S o such as I 2 or a low S o like I 1 . Foam apparent viscosity or pressure gradient (∇p) achieved at steady state may well depend on the initial conditions. Different foam states achieved for the same injection condition would reflect the existence of multiple steady states, that is, J 1 and J 3 . Second, one could fix ∇p and f g across the core. A series of measurements at steady state with increasing ∇p might give a folding curve similar to Figure 2. Both approaches face challenges. Either of the two approaches needs to start with accurate foam-model parameters and relative-permeability curves that fit a specific situation. This is time consuming and difficult in the presence of oil, which may require imaging techniques such as X-ray CT to obtain phase-saturation information. Core-scale artifacts (entrance region for foam generation, capillary end effect) can distort coreflood measurements, even at steady state. Moreover, as noted in the wave speeds reported above, the time for foam propagation through a core can be extremely long, especially in the cases of displacement by a weak or collapsed foam.

Conclusions
The implicit-texture foam model discussed here predicts multiple steady states for foam flow with oil: an upper strong-foam state J 1 , intermediate state J 2 , and lower collapsed-foam state J 3 , with different apparent viscosities but all fitting same set of fractional flows. A displacement could then correspond to more than one possible injection state for the same injected fractional flows. Our study shows how to determine which state makes the displacement and the dependence of the nature of the displacement on initial condition, using three-phase fractional-flow theory and the wave-curve method.
The WCM for three-phase flow in porous media has a different problem definition than a coreflood with specified injection rates. In a coreflood, injection state J is fixed at, x D = 0 displacing forward initial state I, which is present for 0 < x D ≤ 1. In the WCM, J is injected from x D < < 0, with initial state I present for x D > 0 and J for x D ≤ 0.
The distinction between the problem definitions in the WCM and a coreflood makes the WCM capable of identifying the unique displacing state among multiple possible injection states J 1 , J 2 , and J 3 . A composition path, from each of J 1 , J 2 , and J 3 to initial state I, can be constructed in the WCM, satisfying the conservation equations. Only the path with only positive velocities, in the solutions solved by the WCM, corresponds to a physical injection by J maintained at x D = 0. Any composition path featuring negative velocities does not correspond to this physical injection condition.
A stable displacement could be made by either the upper strong-foam state J 1 or the lower collapsed-foam state J 3 , but never the intermediate, intrinsically unstable state J 2 . The choice of the displacing state shows a dependence on initial state. We define a boundary curve through the unstable intermediate state J 2 in ternary saturation space that captures the dependence of the nature of the displacement on initial state I. For any I in the region where the upper state J 1 resides, the strong foam state J 1 makes the displacement. For any I located in the region where the lower state J 3 lies, it corresponds to a stable displacement by the collapsedfoam state J 3 .
The significance and implications of the findings for field applications are discussed. We also give suggestions for the verification of multiple steady states in foam flow with oil through corefloods in laboratory. Further investigations would account for complexities such as capillary diffusion, formation heterogeneity, and gravity separation in a foam EOR process.
Foam in the STARS model scales gas mobility by modifying k rg in equation (6), the foam-free gas relative permeability, through a mobility-reduction factor FM: where superscript f denotes the presence of foam. k f rg with foam is usually referred to as effective gas relative permeability. The mobility scaling factor FM accounts for the effects of a series of physical factors on foam described by functions F n (n = 1, 2, 3 …), for example, surfactant concentration, water saturation, oil saturation, capillary number, shear-thinning, or salinity: where fmmob is the reference mobility-reduction factor, representing the maximum attainable gas mobility reduction.
This study considers two major factors dominating local-equilibrium foam flow behavior, water (S w ) and oil (S o ) saturations. The effect of S w is captured through the function F 2 : F 2 ¼ 0:5 þ arctan epdry S w −fmdry ð Þ ½ π ; (A3) Figure A1. Gas relative permeability k f rg with foam, as represented by the wet-foam model. Oil is introduced by fixing various ratios (u o /u w ), which is equivalent to fixing (f o /f w ). S o increases with S w along each curve at a fixed (u o /u w ), and the fixed value of S gr corresponds to different values of S w on the curves, which depend on the values of (u o /u w ). k rw and k rg are shown for a comparison with k f rg . All model parameters used for the illustration are in Table A1. where fmdry is the limiting water saturation, physically denoted as S * w (Rossen & Zhou, 1995), the water saturation around which foam collapses. Parameter epdry determines the abruptness of foam collapse in a narrow range of S w around fmdry. In the wet-foam model, fmdry is a fixed constant.
Oil here does not have an impact on foam stability (i.e., on fmdry), which dominates the high-quality regime. It affects foam only in the low-quality regime, through its effect on fmmob: where S wc , S or , and S gr represent the residual saturations of water, oil, and gas, respectively. Equation (A4) is a piecewise smooth function depicting the effect of S o through three factors fmoil, floil, and epoil. The factors fmoil and floil denote the upper-and lower-limiting values of S o for stable foam. For S o > fmoil, foam is killed completely (F 3 = 0). For S o less than floil, oil has no effect on foam (F 3 = 1). For floil < S o < fmoil, oil destabilizes foam in a nonlinear way. Table A1 summarizes the model parameters implemented throughout the study. Note that all of our illustrations in ternary diagrams, that is, Figures 7, 9, 10, 13, and 15, are plotted using normalized saturations rescaled as follows: where S j,a with subscript a denotes the absolute saturation of phase j and S jr represents the residual saturation of phase j. Residual saturations of phase j, that is, S wc , S or , and S gr , are not shown in the illustrations. Figure A1 illustrates gas relative-permeability curves on semi-log scale using parameters in Table A1, as described by the wet-foam representation above. Without foam, k rg (green curve) is much greater than k rw (dashed blue curve) for a wide range of S w , indicating a large gas mobility that is usually unfavourable for oil displacements. With foam but without oil, k f rg (cyan curve) drops down suddenly and drastically, by a factor of 10 4 , within a narrow range of S w around fmdry. The abruptness of the drop depends on the value of epdry in equation (A3).
In Figure A1, the effect of oil is introduced by fixing the ratio (u o /u w ), equivalent to fixing (f o /f w ). Note that parameter fmdry, around which k f rg drops abruptly and suddenly, remains independent of the presence of oil and of (u o /u w ). The unchanging value of fmdry reflects the fact that the wet-foam representation describes the effect of oil only on the low-quality regime. For S w > fmdry, increasing the ratio (u o /u w ) at fixed S w raises S o , as seen from equations (4) and (6). That means increasingly weaker foam (smaller value of F 3 in equation (A4)), causing k f rg in Figure A1 to increase nonlinearly. As S o increases enough to kill foam, k f rg rises largely, because gas mobility is no longer restricted by foam. k f rg thereafter drops drastically down to zero as S g approaches its residual saturation.