Universality of Riemann solutions in porous media

Universality, a desirable feature in any system. For decades, elusive measurements of three-phase flows have yielded countless permeability models that describe them. However, the equations governing the solution of water and gas co-injection has a robust structure. This universal structure stands for Riemann problems in green oil reservoirs. In the past we established a large class of three phase flow models including convex Corey permeability, Stone I and Brooks–Corey models. These models share the property that characteristic speeds become equal at a state somewhere in the interior of the saturation triangle. Here we construct a three-phase flow model with unequal characteristic speeds in the interior of the saturation triangle, equality occurring only at a point of the boundary of the saturation triangle. Yet the solution for this model still displays the same universal structure, which favors the two possible embedded two-phase flows of water-oil or gas-oil. We focus on showing this structure under the minimum conditions that a permeability model must meet. This finding is a guide to seeking a purely three-phase flow solution maximizing oil recovery.


Introduction
An effective Enhanced Oil Recovery method is the so-called WAG injection, which consists of injecting water and gas alternately. A modern solution to this problem is the co-injection of water and gas on displacing oil. Both methodologies the WAG and co-injection are indistinguishable, see e.g., [6,36]. In doing so, oil extraction is improved. The solutions of the resulting mathematical model comprise complex series of nonlinear waves (rarefactions and shocks), for which we have identified a persistent structure within a large class of models (see e.g., [6,7,13,14]). In this paper we discuss properties so that a model presents this structure.
In petroleum science there is a fundamental result found by Buckley and Leverett, [9]. The universal behaviour of this solution is valid for the injection of a single fluid (water or gas) for the displacement of oil in a reservoir. Our aim is to contribute to finding persistent solutions for three-phase flow. A persistent solution was found in [6], here we discuss conditions for the persistent solution to exists.
First, we recall the solution structure for WAG injection which we call universal and discuss the sense in which we say that a Riemann solution is robust. The Riemann problems we focus on take the right state as pure oil (resident fluid) and the left state as a co-injected mixture of water and gas. In the saturation triangle (state space), the right state is a vertex and the left is any state on the opposite edge; this is the setting corresponding to the co-injection of water and gas or the so-called sWAG injection (simultaneous water and gas injection in opposition to wateralternating-gas). In [6], it is proved, that because of fluid mixing within the porous rock, the average proportion in a WAG injection is attained in the interior of the porous rock, see e.g. [36]. For this reason, we do not make a distinction in this study between WAG and sWAG injection. For these problems, the solution has, generically, the following structure: it comprises a 1-wave group across which all saturations change preceded by a 2-wave group possessing states along one of the remaining edges (water-oil or gas-oil) of the saturation triangle; the latter wave is thus a two-phase solution (cf. [9]); namely, the solution of a Riemann problem with right state given as a green reservoir and left state consisting of a mixture of either water and oil or gas and oil. We claim that this behaviour is robust in the sense that left states in the water-gas edge in the saturation triangle near the pure water vertex W give rise to solutions with fast waves involving only water and oil, whereas left states near the vertex G give rise to solutions with fast waves involving only gas and oil. Thus somewhere between W and G there exists a special state with the property that nearby left states give rise to Riemann solutions with either water-oil or gas-oil fastest waves, depending on their position in state space relative to such special state. This special state defines the left Riemann state of what we call the separatrix solution. This is the structure of Riemann solutions which we call ''robust''.
We are interested in injection problems leading to flow in porous media of three phases. Such problems are modeled by a system of two conservation laws; a survey of the mathematical theory may be found in [6,10,36] and references therein. For simplicity in the calculations ahead, we assume that the three phases are incompressible, they do not exchange mass, that gravitational segregation and capillary pressure effects are negligible, and that there is no mass transfer among the phases. It seemed important in previous works (cf. [6,7,[13][14][15][16]) the study of immiscible displacement modeled by phase permeabilities which vanish superlinearly as the corresponding phase saturation disappears. Here we generalize these works, [6], by relaxing hypothesis in which the Riemann solution is constructed for data representing sWAG injection into a reservoir initially containing oil and residual water.
The first class of viscosity dependent models is studied in [13], where we extend the earlier work [6] by establishing that the same solution structure for the WAG injection problem in a green reservoir also holds for the broader class of Corey models with convex permeabilities, the second class. (Related preliminary mathematical work appeared in [14][15][16].) The present work relates to a third class of models in [13], which includes Stone I and Brooks-Corey models, [10,23,33]. Other models in this class arise from the experiments reported in [27,54]; for these models, preliminary results confirm a similar solution behaviour. Other cases where the universal solution behaviour is observed are reported in [1,24,36,44,48,55] as well as for the injection of CO 2 and CH 4 in NOTT-400 nanotubes, [25]. Moreover, preliminary results show that the universal solution behaviour persists for SAG injection (surfactant alternating gas) in the presence and absence of foam, see [35,46,53].
The focus in this work is on constructing a permeability model which encompasses the minimal conditions that preserve the universal structure for WAG injection. Thus, we identify in a peculiar model, the minimal conditions for which the universality is preserved.

The universal Riemann solution
In this section we review the solution structure which we dubbed ''universal structure within the WAG injection setting''; for brevity we refer to this structure as U-structure. To justify our claim of ''universality'' (for WAG injection) we notice that, to the best of our knowledge, only the model discussed in [29] does not exhibit this behaviour. (An extensive study on three-phase models can be found in [5,48].) In Sec. 2.1 we discuss the solutions for the first class of models in [13], namely the quadratic Corey model used in [6] to illustrate the U-structure. In Sec. 2.2 we discuss how this solution structure holds for the second and third classes of models in [13] as well, and we exhibit a model for which the U-structure is not observed.
Some of the concepts introduced in a terse way in Sec. 2.1 are rigorously defined only in Sec. 3.1, where the specific construction of certain relevant manifolds and points for the new model given in Sec. 3 are studied.

The U-structure for the quadratic Corey model
Let s w ðx; tÞ, s g ðx; tÞ and s o ðx; tÞ denote the water, gas and oil saturations at position x along the cylinder at time t. Because s w þ s g þ s o ¼ 1 and 0 s w ; s g ; s o , the state space is the saturation triangle; see Fig. 1. In our analysis, it is convenient to choose s g and s o as the two independent variables, and write S :¼ ðs g ; s o Þ to represent states in the saturation triangle. The vertices of this triangle are W ¼ ð0; 0Þ, G ¼ ð1; 0Þ and O ¼ ð0; 1Þ, representing pure water, pure gas, and pure oil. Three-phase flow in one spatial dimension at constant total flow rate is governed by the nondimensionalized system representing conservation of gas and oil (as well as of water since the saturations add to 1). Any of these equations can be replaced by an equation for water conservation with flux f w defined implicitly by Riemann problem with pure oil O ¼ ð0; 1Þ as the right state and left state L ¼ ðs g ; 0Þ on the edge WG. The concept of wave group is crucial for our analysis, see e.g. [31,34]. A particular Riemann problem arises for the left state B defined by the ratio of water and gas viscosities B ¼ ðl g =l wg ; 0Þ with l wg ¼ l w þ l g . The state B denotes one end point of the separatrix OB in Fig. 1.
Definition 1 Consider a (backward wave) curve in the saturation triangle connecting state O to a state B on the edge WG consisting of left states of wave groups with O on the right. This curve is called a separatrix if the solution of the Riemann problem with right state O and left state LðsÞ ¼ B þ ðs; 0Þ has a gas-oil displacing wave for 0\s\e and a water-oil displacing wave for Àe\s\0, for some e [ 0.

Remark 1
The left state B has the peculiarity of also being the optimal mixture in WAG injection, with this the recovery of crude oil is maximized (cf. [6,7,13]).
The sense in which the separatrix is a wave curve has to be seen under the gaze of a few previous investigations. Typically, it is a composition of waves involving both characteristic families. The complete construction of the path in the saturation triangle aids to characterize a scalar version of this curve; for technical reasons it is better to start from O with a backward wave curve. We define, as in [12], the effective flow function f ðsÞ ¼ f o ðs g ðsÞ; sÞ with S ¼ ðs g ðsÞ; sÞ going through the values of the saturation of gas and oil on the separatrix. This effective flux function works as a lifting, so the scalar equation satisfies Oleǐnik principles, cf. [12], guaranteeing among others the admissibility of shocks. In [13] we showed that something analogous occurs for models in the first, second and third classes. Care must be taken, because the shock from B Ã to O, in these cases, is the limit between an overcompressive and a Lax 2-shock, that is, the speed of the discontinuity is equivalent to the slow characteristic speed at B Ã . In [16] there is an analysis on the composition of these shocks as the concatenation of a ''limiting'' Lax 1-shock (to any of the borders) followed by a Lax 2-shock with the same speed, in this way, the three versions are equivalent and it is better to think of the direct shock from B Ã to O. For the general admissibility of shocks we take an superposition of Lax's and Liu's admissibility criteria, see e.g. [13] for a more detailed description.
In this case, B is located when the ratio s g ¼ l g =ðl w þ l g Þ is satisfied. It can be shown that the line connecting B to O is, in this case, invariant under the dynamics of system (1), see e.g. [7,12]. The separatrix is given by the concatenation of admissible shock waves along part of the Hugoniot locus HðOÞ and the slow-family rarefaction curve for the left state B on WG.
In Fig. 1, we illustrate the U-structure depicting Riemann solutions corresponding to left states I Ã , B and J Ã on either side of the separatrix. The solutions are constructed in a similar way for each state on WB or on BG. For example, for a state L on the open segment BG we have the following wave structure: 1. -Slow wave: 2. -Fast wave: which is a condensed scheme for the sequence of waves in the wave groups for the solution of a Riemann problem. From left Riemann data L to right Riemann data O we have that the rarefaction of the slow family (rar1) connects the left state L to a state S on the slow-family extension of edge GO (denoted as I s GO ), and from there a left-characteristic (or sonic) Lax 1-shock (''limiting'' Lax1) connects S to a state M on the edge GO. (This shock has speed equal to the slowfamily characteristic speed at S, [16].) From such state M the solution proceeds along GO to O. In this framework, the fastest wave is this second wave, known as the displacing wave; thus, the breakthrough is done with a two-phase wave.
As those solutions need this second wave, notice that a gas injection is comprised by the fast rarefaction wave from G until G Ã , a state satisfying k s ðG Ã Þ ¼ rðG Ã ; OÞ and well-known on Buckley-Leverett solution, which typically is known also as the Bethe-Wendroff point. From this state the 2-Lax shock to reach O is left characteristic.
The structure of the second wave depends on the position of M relative to the Bethe-Wendroff point G Ã on the edge GO. There are two cases: it is either a fast rarefaction with an adjacent (left-characteristic) Lax 2-shock (''limiting'' Lax2) or a single Lax 2-shock (Lax2). The change of structure occurs at M ¼ G Ã , so it is easy to track the path back to find J Ã as follows: first, S J is found as the state in HðG Ã Þ which comprises a ''limiting'' Lax 1-shock from S J to G Ã and, second, from S J the slow rarefaction is computed the in backward direction, i.e. with decreasing characteristic speed until J Ã ; thus, all states along J Ã and S J belong to a spreading wave. Thus the solutions for injection states L belonging to J Ã G possess in their second wave a fast rarefaction (rar2) with an adjacent Lax 2-shock wave, whereas for injection states L belonging to BJ Ã the Riemann solution possesses a single Lax 2-shock. The solutions for states L in the interval WB are constructed in a similar way with J Ã replaced by I Ã , and G Ã by W Ã .
As pointed out earlier, the solution for the left state B is important because B is the starting point of the separatrix. The line BO is an invariant manifold for system (1), see [6,7]. In particular, the solution for L ¼ B consists of a single wave group; it is the so-called Buckley-Leverett solution.

The solutions within the U-structure
The structure of Riemann solutions described in Sec. 2.1 holds for a wider class of models, see e.g. [1,25,27,36,44,48,54]. The U-structure is observed when, upon choosing the left state, the displacing wave is a two-phase wave which changes type: either water-oil or gas-oil, depending on the location of the left state. In contradistinction, the solution along the separatrix typically is genuinely a threephase flow; it comprises only waves involving all phases. It is clear that the U-structure, corresponding to co-injection of water and gas, requires a specific behaviour. For example, the invariance of the edges is part of the universal structure, see (i) in Def. 2. This is a natural assumption since in the absence of a fluid, the evolution of the system should not produce such a fluid, so that it is absent throughout the process.
From a historical point of view, in the late 70s a condition for permeabilities (or mobilities) to vanish superlinearly was imposed to permeability models, see e.g. [19,28,52] and references therein. Previously, in 1941, Leverett and Lewis, [32], performed the first recorded experiments for three-phase flows permeabilities, indicating a superlinear vanishing permeabilities. In 1927, the pioneering work of Kozeny, [30], and later the review in 1937 of Carman, [11], led to the known Carman-Kozeny equation for a single fluid in a bed of solids formed by parallel capillary tubes. This equation is, to the best of our knowledge, the first analytical explanation of the superlinear vanishing condition. The modification of Kozeny-Carman equation to multiphase flow was first estimated by Purcell, [43], and continue e.g. with the efforts of Alpak et al., [4]. In our review, the culminating research was done in the 90s, with e.g. the analyses carried out as in [26] provided a set of mild assumptions to be a reasonable set of conditions on permeabilities, see also [28]. In all cases (cf. [19,26,28]), the set of conditions guarantees the invariance of the edges with the fast-family characteristic direction parallel to them. This set of conditions causes the slow-family characteristic speed to be zero at the boundaries and to increase into the interior of the state space. In such circumstances, an interior branch of the Hugoniot locus emanates from each of the respective vertices. Suggesting the formation of the U-structure.
Typically, experiments manage measurements for two-phase flow only, thus three-phase flow are extrapolations from such data inheriting the superlinear vanishing of relative permeabilities. Nonetheless, experimental data with superlinear vanishing have a vast literature, see for example [17,20,38,47,51] for light oils or the recent PhD thesis [2] for heavy oil (cf. [3]). Other studies are carried out in [1,48] and references therein.
In the past few years, the superlinear vanishing condition has been inquired for gas or water permeabilities. The denial of such property may lead to a strictly hyperbolic system in the interior of saturation triangle. A strictly hyperbolic system typically admits a global curvilinear system of coordinates defining both families of integral curves, which cannot possesses universality [31,34]. On one hand, Shearer and Trangenstein, [49,50], established that violation of strict hyperbolicity occurs inside the saturation triangle for immiscible three-phase flow. On the other, Medeiros, [37], proved that a hyperbolic singularity is a consequence of two-phase behavior near the saturation triangle for immiscible three-phase flow. As Falls and Schulte studied, [21,22], the role that umbilic points play are as organizing centers for bifurcations, determining some aspects of the wave structure in Riemann solutions. In this work, we aim for a better understanding of the solution behaviour for WAG injection. We show that this vanishing condition is not needed, and that a loss of strict hyperbolicity is required at least on the boundaries of state space for the occurrence of the U-structure.
A direct conclusion from Def. 2 is that the family of Riemann solutions has the U-structure when there are two distinct solutions along both of the edges of saturation triangle for the system of conservation laws and there exists a single separatrix curve. Nonetheless, properties guaranteeing the existence and uniqueness of a separatrix are elusive. In Sec. 2.1 we showed that the Hugoniot locus for right Riemann datum O has three branches: the edges WO and GO as well as the inner branch connecting B with O, the latter is part of the separatrix. Once the separatrix is given it seems natural to follow the steps shown in Sec. 2, however there are some particularities that can be misleading.
The U-structure has been observed for models which fail to be strictly hyperbolic in the interior of the saturation triangle, viz. convex Corey, Stone I and Brooks-Corey models (see [6,7,13,24,33]). On the other hand, in [29] an example was exhibited without the U-structure when the interior umbilic point coincides with one of the umbilic points at the vertices. In the mean time, we seek a strictly hyperbolic model (at least in the interior of state space) that displays the U-structure.
The model in the current work only fails to be strictly hyperbolic at an umbilic point located on one edge of the saturation triangle (and at two vertices), it is strictly hyperbolic in the interior of state space. This model will aid in identifying the minimal conditions for the U-structure to exist. In the next section we construct a special case for a permeability model possessing the U-structure, that is strictly hyperbolic on the interior of the saturation triangle and for which the Bethe-Wendroff point B Ã belongs to a ''nonlocal'' branch of HðOÞ.
We construct a three-phase flow model in which the phase mobilities violate the assumptions in [26]; yet, it still possesses the U-structure. Our goal is to construct a strictly hyperbolic system inside the state space by displacing the umbilic point to one of the edges of the saturation triangle.
To keep the model as simple as possible, we take the gas and oil mobilities as in quadratic Corey models. For the sake of clarity, we recall them here, For the water mobility we introduce where and for a constant j yet to be defined. Our motivation lies in modifying as minimum as possible a model that allows analytical computations, the Corey quadratic model, [7]. From a phenomenological point of view, in the so-called oil-wet reservoirs, water surrounds the gas and is wrapped by oil within a porous medium, thus revealing a promising candidate for permeability alterations. Note that the water mobility does not vanish superlinearly as s w ! 0 since om w ðSÞ=os w 6 ¼ 0 for all states satisfying s w ¼ 0. Indeed, which guarantees linear vanishing if j 2 ðÀ1=l g ; 1=l o Þ.
Again, for simplicity we take j ¼ 1=l o , which yields linear vanishing except at s o ¼ 0. We also choose l w ¼ 1, l g ¼ 1=2 and l o ¼ 2 to produce the figures shown in this work. More realistic viscosity values yield similar results.
There are several distinct features in the present model (3)-(4) in comparison to Corey models with quadratic mobilities. First, the vertex O is no longer an umbilic point (vertices W and G still are). Second, the model is strictly hyperbolic in the interior of the saturation triangle, and it loses strict hyperbolicity at the umbilic point U in the edge GO, see Fig. 2. The slow family characteristic speed increases away WG and WO edges with transversal corresponding eigenvector. (The eigenvector is transversal in GU but the characteristic speed decreases away this segment.) We will prove these assertions together with invariance of the edges in the next section.

Fundamental waves
Equations (1) have solutions that propagate as nonlinear waves. Because of selfsimilarity of the data and the PDE, the solutions of a Riemann problem depend on x/ t and consist of centered rarefaction waves, shock waves and sectors of constant states (the so-called intermediate states), see e.g. [18,31,34,39]. The characteristic speeds are the two eigenvalues of the Jacobian matrix where for S ¼ ðs g ; s o Þ. (Here, a prime 0 denotes total differentiation and o a is partial differentiation with respect to s a , a ¼ w; g; o, as the water mobility depends on all saturations.) We refer to the smaller and larger eigenvalues as the slow-family and fast-family characteristic speed, respectively. System (1) has two families of continuous solutions called slow-and fast-family rarefaction waves. They arise by solving an ODE, namely, dS=dn ¼ r k ðSÞ, k ¼ slow or fast , defined by the eigenvectors r k ðSÞ of the Jacobian matrix JðSÞ, with fJðSÞ À n Igr k ðSÞ ¼ 0. For n ¼ x=t, the profile SðnÞ defines the forward rarefaction, provided n is monotone increasing; it is called backward rarefaction for n monotone decreasing. Some rarefaction curves, which parameterize the rarefaction waves appearing in the solution of our Riemann problems, are calculated with the ELI package (cf. [42]) and plotted in Fig. 2. System (1) also admits solutions in the form of moving jump discontinuities. To respect conservation of gas and oil (and hence water), the fluxes in and out of the moving discontinuity must balance. In terms of the state S À ¼ ðs À g ; s À o Þ on the left of the discontinuity, the state S ¼ ðs g ; s o Þ on the right, and propagation speed r, this balance is expressed as FðSÞ À rS ¼ FðS À Þ À rS À , or Equations (7) are the Rankine-Hugoniot (RH) conditions. The Hugoniot locus of a state S À , denoted by HðS À Þ, consists of all states S for each of which there exists a value of r ¼ rðS À ; SÞ such that the RH conditions (7) are satisfied. Thus a shock wave ðS À ; SÞ with speed r is classified with the Lax admissibility criterion as slowand fast-shock for Lax 1-and 2-shock respectively (cf. [31]). An analogous construction follows for S þ on the right of the discontinuity. As the mobilities m g , m o vanish superlinearly as s g or s o tend to zero, respectively, it is easy to prove that the Jacobian in (6)  is satisfied, so det JðS; s w ¼ 0Þ is strictly positive. Along this edge, the eigenvalues are distinct, except at U ¼ ðl g =l go ; l o =l go Þ 2 GO where the Jacobian is 2 times the identity matrix, and both eigenvalues are equal so that U is a umbilic point. As in convex Corey models, here the umbilic point represents the maximum pressure gradient or the minimum of total mobility; the proportion of saturation at which each of the three fluids hinders maximally the flow of the other two (cf. [15].) This further ensures that the umbilic point is isolated even when extending the domain outside the saturation triangle. Examples of these extensions can be found in [8]. The first work in our knowledge and of which we are aware of the existence of quase-ubilic points appear in [44,45]. A direct inspection of (6) shows that v s ¼ ð1; holds for S 2 GU (and the opposite inequality for S 2 UO), we conclude that v is associated with the fast-family characteristic speed for the former and with the slowfamily characteristic speed for the latter. The other eigenvector is transversal to GO, except at U, which is an umbilic point. This analysis confirms our numerical plots in Fig. 2.
Remark 2 Some relevant facts are (refer to Fig. 2): (1) the slow-family integral curves are parallel to the edge GO and point to U on the segment UO, and they are transversal to the segment GU;

Relevant manifolds and points in state space
Now we analyze important features of the model under study. Hence the manifolds and points in Fig. 3. First, we calculate points on edge WO and their relevant connections to points on WG. We also locate two important manifolds. We compute and denote the respective state coordinates for our parameter election, see them after equation (5).
Þ The edge GO is invariant, waves along this edge can be treated by a single conservation law with flux restricted to the sonic shock connecting G Ã with O has speed k f ðG Ã Þ; in other words, G Ã is a Bethe-Wendroff point. A simple calculation shows that the analogous point on the edge s g ¼ 0 does not exist since the flux along that segment (with s o as variable) is convex, with zero derivative only at W, i.e., at s o ¼ 0.
From G Ã we construct the backward slow-rarefaction curve which intersects the edge WG at a state denoted by J Ã .
Secondary bifurcation point C ¼ ð1=3; 2=3Þ The RH condition (7) for S þ ¼ O shows that the edges s g ¼ 0 and s w ¼ 0 are part of the Hugoniot locus HðOÞ. Equating and eliminating the speed r in equations (7) with S À ¼ ð0; 1Þ yields to relation s g m w ¼ s w m g , which defines the ''nonlocal'' branch of HðOÞ. This branch is a line segment, the thin line in Fig. 3, with end points The shock ðC; OÞ is admissible since C is closer to O than G Ã . From C we construct the backward slow-family rarefaction curve to edge WG reaching a state denoted as B.
Umbilic  (7), rðS; OÞ ¼ f g ðSÞ=s g for any given S 2 HðOÞ, so that the following shock speeds are obtained using expression (4). Thus ðG Ã ; OÞ is a sonic Lax 2-shock and ðC; OÞ a limiting overcompressive/Lax 2-shock; actually C belongs also to the slow-family extension of WO boundary.

WO
Analytical construction of extension boundaries is given in [16]. We remark that both C and W belong to I s WO . Moreover, this manifold is beneath the slow inflection showed in Fig. 2 and any S state on it satisfies k s ðSÞ ¼ rðS; MÞ for M 2 HðSÞ \ WO, see the dotted curve in Fig. 3.
Notice in Fig. 3 the relative ordering in edge s w ¼ 0 of states G, G Ã , C, U, and O. The depicted configuration maintains the relative positions between C and U as the inequality 2ð1 þ jl g Þ [ 1 þ 2jl g holds for any choice of parameters. However, the configuration may have G Ã and C reversed. Considering, for example, j ¼ 1=l o and l o [ 4=3, it is enough that l g \l o is satisfied to have the realative positions of the Fig. 3. It is important to note that G Ã cannot change positions with U, as the ratio l g =l go is always less than one. The change of positions between G Ã and C alters the wave groups conforming the separatrix, but the results remain similar; we use the chosen configuration since pretends to be more realistic.

Solution construction
We will solve the Riemann problem with left state L on segment WG and right state O. The third fact in Remark 2 about both slow-and fast-family characteristic speeds near O shows that, as right state, O can be reached by a rarefaction only along the edge WO; thus it is important to understand its Hugoniot locus, HðOÞ, see Sec. 3.1.1.
Since the edges GO and WO are invariant, a Buckley-Leverett solution-like may be constructed. (Notice that the flux on WO is convex; thus the solution along WO is not the classical BL solution: it is a single rarefaction wave.) A Bethe-Wendroff point G Ã exists in GO, however in WO there is not an analogous W Ã ; all states in the latter reach O by a single fast-family rarefaction wave. The analogous point to B Ã in Sec. 2 now belongs to the edge GO; actually, it coincides with the secondary bifurcation point C. System (1) reduces to a scalar conservation law along the edge GO and every state on segment G Ã O is the left state of an admissible Lax 2-shock to the right state O. There are no other states reaching O with admissible shocks. Special attention for states on segment DC is needed. A possible discontinuity with O as right state and S on DC as left state has shock speed equal to the slow-family characteristic speed at S, which has ''limiting'' Lax 2-shock classification. However, numerical simulations show that such discontinuities are not a viscous profile solution, therefore are inadmissible shocks. Thus a solution profile emanating from WG must reach a state on one of the edges before proceeding to O.
An analogous curve to those invariant edges is constructed as the separatrix. Notice we can use the admissible shocks in segment CO and continue a curve by means of the backward slow-family rarefaction to the edge WG at B. Thus, for state O as the initial reservoir state, the separatrix can be regarded as the set of injection states which give rise to an injection solution of (1)-(2) comprising states over this curve, see Fig. 3a, 3c.
We can even construct an effective flux function over the separatrix as in [13] (see also [12]) with the oil saturation as the wave group parameter. This effective flux function comprises the slow rarefaction from B to C, then the limiting overcompressive/Lax 2-shock to O, and all states on CO (see Fig. 4 ahead). Indeed, C belongs to I s WO and the shock speed satisfies rðC; OÞ ¼ k s ðCÞ [ k f ðOÞ. In the same way, we denote by J Ã 2 WG the starting point of a slow-family rarefaction reaching G Ã .
The other slow wave groups follow directly and in Fig. 3 the solutions for left states in WG are given by 1.-Slow wave: for L 2 BG 8 > < > : : 2.-Fast wave: Clearly, the solution structure bifurcates at B and J Ã as left states, see Fig. 3. The class of solutions are regular in the sense that the L 1 loc -continuity holds as the left state passes from segment WB to segment BG. For the former, say with L ¼ I 2 WB (see Fig. 3 where the characteristic shock to M 2 WO has speed k s ðSÞ, M is an intermediate state and the 2-wave comprises the fast-family rarefaction from M to O. As L approaches B, S and M approach C and O respectively, thus the shock speed rðS; MÞ approaches k s ðCÞ, and ðC; OÞ is a ''limiting'' Lax 1-shock which speed is characteristic at the left of the discontinuity and is larger than the impinging characteristics with speed k f ðOÞ at the right of it. For the latter, with L 2 BJ Ã , the solution comprises the slow-family rarefaction to an intermediate state M 2 G Ã C and a Lax 2-shock to the right state O. As L approaches B, M approaches C and the shock speed approaches k s ðCÞ; a characteristic shock which turns out to be a limiting overcompressive/Lax 2-shock. Thus the L 1 loc -continuity holds for L near B. A simpler analysis follows for L near J Ã : the fast-family rarefaction and adjacent to the ''limiting'' Lax 2-shock collapses into a single Lax 2-shock as M approaches G Ã , with characteristic speed k f ðG Ã Þ. Remarkably, the L 1 loc -continuity holds for the solutions with L near B with a simpler analysis than the triple-shock rule used in [13] (cf. [16]) and it also holds for L near J Ã . Thus, we conclude that the L 1 loc -continuity of solutions holds among the left states on WG.

The effective flux function
It is clear that the presence of a separatrix curve is a crucial ingredient for the occurrence of the U-structure. In all the models previously analyzed the separatrix comprises the left states of admissible shocks with left state as an interior point of the saturation triangle and the vertex O as right state. In [6,7,13,14] the construction of such separatrix is straightforward with the use of the ''internal Hugoniot'', e.g. in [13] we take the interior part of the Hugoniot locus HðOÞ and look for its Bethe-Wendroff point, the state from which we construct a backward slow-family rarefaction curve that reaches the edge WG. Here the concept of ''interior part of the Hugoniot locus HðOÞ'' must be revisited.
In all previous models studied by us, see [13,14], the Rankine-Hugoniot relation (7) when S þ ¼ O gives three branches that emanate from vertex O; it is an umbilic point. Now O is ''regular'' (i.e., where strict hyperbolicity holds) and its RH locus has only two local branches. However, the Hugoniot locus of O has a secondary bifurcation point at C. Thus, C belongs to both the ''nonlocal'' branch and a local one of the Hugoniot locus HðOÞ. A simple analysis shows that shocks with left state on the half-open segment DC are not admissible for solutions since Lax admissibility conditions are fulfilled but they do not possess a viscous profile. This is not the case for C, which also belongs to the set of Lax 2-shocks with right state O.
It is well known (cf. [7]) that for the quadratic Corey model, the slow-family integral curves form a regular foliation of the saturation triangle everywhere far from the umbilic points, both interior and the vertices. (For the present model, the foliation is regular in the interior of the saturation triangle; the umbilic points are at the boundary.) Since k s ðCÞ ¼ rðC; OÞ holds but the effective flux function along the separatrix has distinct left and right derivatives at C, we can regard C as a generalization of a Bethe-Wendroff point (with respect to the slow-family), and the backward slow rarefaction curve starting at C reaches the edge WG at B, the starting point of the separatrix. Actually, from left states on BG the slow rarefaction curves foliate the state space, reaching segment GC (see left panel in Fig. 2). Thus a solution with left Riemann data B satisfies O: Since C belongs to I s WO , rðC; OÞ ¼ k s ðCÞ holds, thus the shock ðC; OÞ is also viewed as a left-characteristic overcompressive shock.
We show that such a concatenation of curves is actually a separatrix. In Fig. 4 we draw the flux functions f o restricted to the edges GO and WO; we also draw the lifting of the separatrix solution (cf. [12]) which acts as a flux function along the separatrix. Notice in Fig. 4 that C belongs to two curves, the slope of their tangents represent the slow-and fast-family characteristic speeds at C. The shock speed rðC; OÞ is k s ðCÞ, thus the secant line connecting C with ð1; 1Þ in Fig. 4 is tangent to the continuous curve at C. In the same way, the shock speed rðG Ã ; OÞ is k f ðG Ã Þ, thus the secant line connecting G Ã with ð1; 1Þ is tangent to the stripped curve at G Ã .
An important feature of Fig. 4 lies to show equal-speed shocks that are seen as a single shock in the physical space. This follows from the triple shock rule (cf. [7,13]). First of all, notice that the point ð1; 1Þ represents the right state O; second, notice that points over the striped or dashed curves in Fig. 4 represent left states along HðOÞ, and third, following Oleǐnik construction, the slope of a secant passing through such points and O is the shock wave speed. Thus, if a secant slashes more than one flux function, each intersection represents a left state on an invariant edge or on the separatrix. These states belong to a single Hugoniot locus with equal shock speed. In addition, the tangents to the flux functions show the location of the Bethe-Wendroff points G Ã and C.
4 Discussion on the criteria for models with universal structure In Sec. 2.1 we discussed the universal structure for a specific permeability model. The same structure arises in more complex models as we have pointed out in Sec. 2.2. These models lose strict hyperbolicity somewhere in the interior of the saturation triangle. The model introduced in Sec. 3 is strictly hyperbolic everywhere in the interior of the saturation triangle with strict hyperbolicity failing at two vertices and at one point on an edge. We have empirically observed that failure of strict hyperbolicity in the closed saturation triangle is a necessary condition for a model to posses the U-structure. We further conjecture the requirements for a sufficient condition as follows: Conjecture 1 Consider the 2 Â 2 system of conservation laws given by (1), (2) in the saturation triangle. The Riemann problem with data L 2 WG and R ¼ O possesses the U-structure (Def. 2) if the following two conditions hold: (i) The equation s w m g ðSÞ À s g m w ðSÞ ¼ 0 has a unique solution on the open segment WG. That is to say, one of the branches of HðOÞ intersects the edge WG precisely at a stateB ¼ ðb g ; 0Þ with b g 2 ð0; 1Þ. (ii) There exists a state B on the open segment WG such that the slow-family rarefaction emanating from B reaches a state B Ã on the branch of the Hugoniot locus HðOÞ that containsB and such that k s ðB Ã Þ ¼ rðB Ã ; OÞ.
The first condition in Conjecture 1 is easy to evaluate and it also provides reasons for considering the second condition. From the RH conditions (7), a state S in the interior branch of HðOÞ satisfies s g m o ðSÞ ¼ s o m g ðSÞ, and a simple manipulation leads to condition (i) for S ¼B, once we recall that m ¼ m w þ m g þ m o . (Another way to view this is to use the triple-shock rule, cf. [15].) Thus, condition (i) guarantees thatB belongs to HðOÞ and condition (ii) allows constructing the separatrix curve.
Conjecture 1 guarantees the existence of a separatrix. Regardless of the origin of the permeabilities, or for given flow functions, we have a simple test that models in petroleum engineering must fullfill to have a U-structure. However, it does not guarantee that it is unique. In addition, we must prove that all states along this separatrix form a set of invariant states in the sense that their respective admissible solutions comprise only states over this separatrix. These two prerequisites are crucial for the U-structure (see Def. 2). Regarding the admissibility of the Riemann solution, the construction of an effective flux function (as in Sec. 3.3) appears to be a promising path. Uniqueness is tougher, since a topological analysis is missing. We need to guarantee that any Riemann solution with right state as pure oil O has: (1) for left state on WB, its displacing wave along WO as well as (2) for left state on BG, its displacing wave along GO.
We remark that it is also possible that B Ã does not exist; for example, when all states on the internal branch of HðOÞ are related to admissible shocks to O. In such a circumstance condition (ii) must be read as the shock ðB; OÞ being admissible.

Concluding remarks
In Sec. 2.2 we outlined sufficient conditions for the U-structure to hold: for the right state at one vertex, namely O, it was pointed out that fast characteristic vectors must be parallel to the respective edges, with decreasing speed in the oil direction and, of course, the construction of the separatrix is needed. The model discussed in Sec. 3 had mild conditions which inspired us to coin Conjecture 1. Notice also, that this model displays the universal behavior despite failing to globally satisfy the Ustructure, which is made evident for solutions with left Riemann datum in WO and right Riemann datum at G. Indeed, condition (i) of Conjecture 1 fails for any left state in WO such that the slow wave reaches WG with a slow-family characteristic shock wave; the separatrix solution does not exist. On the other hand, with a left Riemann datum in GO and right datum at W the U-structure follows trivially: condition (i) is satisfied, moreover the segment UW is invariant and works as the separatrix solution, following the construction in [6,13]. Conjecture 1 is satisfied by a wide class of permeability models, as was shown in [13], in addition to the new model given in Sec. 3. We also remark that physical derivations of permeability functions usually lead to these conditions, see e.g. those analyzed in [27,54] which will be tackled in future work. Notice that the set of reasonable physical conditions given by the two-phase flow over boundaries WO and GO and the immiscibility of the fluids which leads to a S-shaped two-phase flow for the water/gas mixture imply that condition (i) of the Conjecture 1 is fulfilled. Therefore, models as those presented in [24,33] also lead the to U-structure and, the study presented in [21] shows that condition (i) of Conjecture 1 is satisfied, promoting the robustness of the U-structure. Moreover, it is noteworthy that the U-structure was also reported for a model including gravity effects, [44], but only for the cases in which condition (i) is satisfied. Hence the strength of this Conjecture.
Since WAG injection converges to a sWAG injection with the corresponding proportion of fluids injected (cf. [6,36]), the separatrix existence is encouraging. On the one hand, the optimal state of injection is located at the end of this curve on WG, and is the one that maximizes the recovery of oil (cf. [13]), on the other, controlling how to reach this situation is due to increasing the injection of the absent fluid to that produced in the leading wave.