Capillary instability in models for three-phase flow.
Standard models for immiscible three-phase flow in porous media exhibit unusual behavior associated with loss of strict hyperbolicity. Anomalies were at one time thought to be confined to the region of nonhyperbolicity, where the purely convective form of the model is ill-posed. However, recent abstract results have revealed that diffusion terms, which are usually neglected, can have a significant effect. The delicate interplay between convection and diffusion determines a larger region of diffusive linear instability. For artificial and numerical diffusion, these two regions usually coincide, but in general they do not.Accordingly, in this paper, we investigate models of immiscible three-phase flow that account for the physical diffusive effects caused by capillary pressure differences among the phases. Our results indicate that, indeed, the locus of instability is enlarged by the effects of capillarity, which therefore entails complicated behavior even in the region of strict hyperbolicity. More precisely, we demonstrate the following results. (1) For general immiscible three-phase flow models, if there is stability near the boundary of the saturation triangle, then there exists a Dumortier-Roussarie-Sotomayor (DRS) bifurcation point within the region of strict hyperbolicity. Such a point lies on the boundary of the diffusive linear instability region. Moreover, as we have shown in previous works, existence of a DRS point (satisfying certain nondegeracy conditions) implies nonuniqueness of Riemann solutions, with corresponding nontrivial asymptotic dynamics at the diffusive level and ill-posedness for the purely convective form of the equations. (2) Models employing the interpolation formula of Stone (1970) to define the relative permeabilities can be linearly unstable near a corner of the saturation triangle. We illustrate this instability with an example in which the two-phase permeabilities are quadratic.Results (1) and (2) are obtained as consequences of more general theory concerning Majda-Pego stability and existence of DRS points, developed for any two-component system and applied to three-phase flow. These results establish the need for properly modelling capillary diffusion terms, for they have a significant influence on the well-posedness of the initial-value problem. They also suggest that generic immiscible three-phase flow models, such as those employing Stone permeabilities, are inadequate for describing three-phase flow.
Introduction
The standard model for three-phase
ow in petroleum reservoirs is based on the two-component
system of partial dierential equations
@
@t
@x
@x
@
@x
U
Date: November, 1999.
This work was supported in part by: FAP-DF under Grant 0821 193 431/95; CAPES under Grant
BEX0012/97-1; FEMAT under Grant 990003; CNPq under Grant CNPq/NSF 910029/95-4; CNPq under
Grant 520725/95-6; CNPq under Grant 300204/83-3; MCT under Grant PCI 650009/97-5; FINEP under
Grant 77970315-00; CNPq under Grant 301411/95-6; NSF under Grant DMS-9732876; DOE under Grant
DE-FG02-90ER25084; ONR under Grant N00014-94-1-0456; and NSF under Grants DMS-9107990 and
DMS-9706842.
denote the saturations of two of the
three-phases (viz., water, gas, and oil), f 1 and f 2 are their respective fractional
ow functions,
and B encodes the eects of capillary pressure. The variables u 1 and u 2 take values on the
saturation triangle  := f g. A physical derivation of this
model is given in Sec. 4.1.
It is standard practice to neglect the capillary terms and consider Eq. (1.1) from the
point of view of hyperbolic conservation laws. A well-known diculty with this approach is
that, generically, three-phase
ow models give rise to elliptic regions, i.e., regions where the
eigenvalues of the Jacobian F 0 (U) are not real [5, 37, 38, 41, 20, 21, 30]. In these regions, the
initial-value problem for the equations with B(U)  0 are ill-posed, and numerical solutions
of them display rapid oscillations with large-amplitudes on ner and ner scales as the mesh
spacing tends to zero. However, at least at one time, there seems to have been a general
feeling that the conservation laws should yield a satisfactory hyperbolic theory outside of
the elliptic region. See, e.g., Refs. [18, 19, 34, 14, 13], in which the structure of Riemann
solutions is studied for certain Stone models in the absence of capillarity eects.
The role of elliptic regions in three-phase
ow models has been investigated by several
authors. Bell, Shubin, and Trangenstein studied numerically the Riemann solution of a
version of Stone's model [5]. The model they considered has a long, thin elliptic region,
which appears to be unrelated to the elliptic region that occurs in the present paper [25].
Keytz [24, 25, 26] studied rigorously the wave structure for a model having a strip as an
elliptic region. She explained many of the wave features observed near the boundary of
the elliptic region of Ref. [5]. Riemann solutions for other systems of conservation laws
having innite strips as elliptic regions have been studied extensively by Slemrod [39] and
Shearer [36]. In all these works, simplied or numerical diusion terms were utilized to
determine admissibility conditions for shock waves.
More recently, results in Refs. [27, 8, 2, 3, 17] have emphasized a more primary cause of
anomalous behavior: linearized instability of the complete diusive system, Eq. (1.1), rather
than ill-posedness of its purely convective form. These two concepts coincide in the case of
articial diusion B  I, or more generally when convective and diusive eects commute;
however, as shown by Majda and Pego, the role of dissipation structure is quite important
when other than such special diusion terms are considered. In Ref. [29], they developed
a useful sucient condition in terms of B(U) and F 0 (U) for the linearized instability of
Eq. (1.1). The corresponding \Majda-Pego instability region" contains, but is typically
larger than, the elliptic region. Nonuniqueness [1, 3] or nonexistence [6, 7, 33] of solutions
for Riemann problems, the latter manifested as highly oscillatory waves (which are measure-valued
solutions) [17, 16, 15], can occur in the Majda-Pego instability region, even in zones
of strict hyperbolicity.
Accordingly, we study here a model with a physically correct diusion matrix, taking full
account of capillarity. The resulting diusion matrix is not a multiple of the identity matrix
and gives rise to the eects mentioned above. Notice that numerical diusion typically
commutes with convective terms (indeed, for standard schemes the diusion matrix B(U) is
a polynomial in F 0 (U)); hence these eects (having to do with noncommutativity) do not
arise or are not captured by standard hyperbolic dierence schemes. Of course, numerical
simulations which fully resolve the parabolic equation (1.1) do capture the eects discussed
here.
More precisely, we study a model with (a) permeabilities obtained using Stone's interpolation
formula from quadratic two-phase permeabilities and (b) Leverett's capillary pressure
functions. Stone's model is widely used in petroleum engineering. The purely convective
form of this model generically features an open elliptic region in the interior of the saturation
triangle  along with three isolated points of nonstrict hyperbolicity (umbilic points) at the
corners of . For articial or numerical diusion, the linear instability region consists of the
elliptic region and the three corners. We nd that, when physical (viz., Leverett) capillarity
is taken into account, each of these regions of linear instability expands into the region of
strict hyperbolicity.
In particular, we identify model parameters for which there are wedges at the corners
inside (respectively, outside) of which the model is linearly unstable (resp., stable). This
phenomenon is rather unexpected, since the instability occurs near the boundary, where
modeling should be relatively accurate and the system might be expected to be well-behaved.
On the other hand, if a particular Stone's model is well-behaved near the boundary, or
the corner instability regions are suciently small, we show that a Dumortier-Roussarie-
Sotomayor (DRS) bifurcation point must exist in the interior of . Such a point lies on the
boundary of the Majda-Pego instability region (see the denition in Sec. 2) and generically
lies in the strictly hyperbolic region (where the eigenvalues of F 0 (U) are real and distinct).
As shown in Ref. [3], a nondegenerate DRS point gives rise, in its vicinity, to multiple
solutions of Riemann problems for the purely convective form of the governing equations.
Both linear instability and nonuniqueness are outside the range of \good" behavior for
hyperbolic conservation laws.
The paper divides into two parts. In the rst part, consisting of Secs. 2 and 3, we investigate
existence of DRS points in general two-component models following a degree-theoretic
approach. Our main result is to show that generically there exists a DRS point within any
simple closed curve on which: (i) all points are strictly hyperbolic and Majda-Pego stable;
and (ii) the eigendirections of F 0 (U) rotate by an odd multiple of  as is traversed. This
greatly generalizes an existence theorem proved in Ref. [3] for special, quadratic
ux models
by direct calculation. The proof relies on a detailed decomposition of the Majda-Pego
instability region.
In the second part, consisting of Secs. 4{6, we consider the three-phase
ow model described
above. We determine the Majda-Pego stability of edge and corner points and calculate
the rotation of the eigendirections of F 0 (U) around a contour near @, thereby
verifying the hypotheses for existence of DRS points. At the same time, we give a detailed
discussion of capillarity and its relationship to well-posedness of Eq. (1.1), which we hope is
of general use to the reader.
2. Linear Stability and DRS Bifurcation
Consider a general two-component system of conservation laws of the form in Eq. (1.1),
where F is C 3 and B is C 2 in a particular open subset U of state space; further, assume
that B is strictly parabolic in the sense that, for each U 2 U , the eigenvalues of B(U) have
strictly positive real parts. In this section, we explore a connection between linear stability of
a constant solution and DRS bifurcation via the stability condition of Majda and Pego [29].
2.1. Majda-Pego stability condition. We begin by recalling, and slightly extending, the
results of Ref. [3] on bifurcation of admissible shock waves. An admissible shock wave is a
traveling wave solution U(x; t) of system (1.1). If  denotes its speed and U(x;
x !1, then such a traveling wave corresponds to an orbit of the system of two ordinary
dierential equations
_
which we regard as a dynamical system parameterized by (U
In Ref. [3], the bifurcation of nonclassical (overcompressive and transitional) shock waves
from a constant solution was investigated. It was shown that this bifurcation corresponds
to a codimension-three bifurcation, studied under certain nondegeneracy conditions by Du-
mortier, Roussarie and Sotomayor [12, 11], that occurs at a DRS point, viz., a point (U  ;
U  R satisfying:
(D1)   is an eigenvalue of F 0 (U  ), with associated right and left eigenvectors r  and '  ;
Associated with DRS bifurcation are the phenomena of nonuniqueness of Riemann solutions
and nontrivial asymptotic behavior, as described in the introduction.
Remark. We have substituted the condition ( e
for the original condition
One consequence of our analysis below is that, given condition (D1), conditions (D3) and ( e
are equivalent except if U  is an umbilic point, i.e., F 0 (U
will prove
to be more convenient in calculations.
One of the main observations of Ref. [3] is a link between DRS bifurcation and linear
instability of the constant solution U(x; t)  U  of system (1.1). This link, however, was
explored only when U  is a point of strict hyperbolicity. Our rst task, therefore, is to establish
a relationship between DRS bifurcation and linear instability when strict hyperbolicity
fails at U  .
By Fourier analysis, L 2 linearized stability of the constant solution U(x; t)  U  of system
(1.1) requires the following condition:
for all real k, each eigenvalue of ikF 0 (U  nonpositive real part. (2.2)
Standard matrix perturbation theory shows that the eigenvalues
s (k) and
f (k) of the
have the following expansion around
where a
s and a
f denote the eigenvalues of F 0 (U  ). Moreover, if the eigenvalues a
are
distinct, and '
denote left and right eigenvectors associated to a
normalized so that
Therefore necessary conditions for linear stability are hyperbolicity at U  and, in case U  is
a point of strict hyperbolicity,
Denition 2.1. We will refer to the inequality
as the Majda-Pego stability condition for family j, and we let the Majda-Pego stable region
MP be the set of strictly hyperbolic points U  for which the Majda-Pego stability condition
holds for both families.
Evidently, condition (D3), in the context of condition (D1), represents neutral failure of the
Majda-Pego stability condition (2.6). As mentioned above, condition (D3) can be replaced
by condition ( e
within the region of strict hyperbolicity [3]. Conditions (D1) and ( e
imply that the trace and determinant of B(U  Following Ref. [8],
we dene the Bogdanov-Takens locus (or BT locus) as
Bogdanov-Takens bifurcation in the structure of weak traveling waves occurs at points of BT
satisfying certain nondegeneracy conditions [25, 8, 9]. Similarly, we dene the coincidence
locus, where the characteristic speeds coincide, to be
It was observed in Refs. [8, 3] that if U is a point of strict hyperbolicity and U 2 @MP ,
then U 2 BT . We will show in this section that the assumption of strict hyperbolicity is
super
uous; in other words, the boundary of the Majda-Pego stable region consists entirely
of points satisfying conditions (D1) and ( e
D3). To prove this statement, let us introduce the
following real-valued functions dened for U
The signicance of E and BT is that solves the equation tr[F 0 (U)
and solves the equation tr[B(U) 1
Lemma 2.2. Let U 2 U . If a s and a f denote the eigenvalues of F 0 (U), then
Moreover,
is 1=4 times the discriminant of F 0 (U), and if F 0 (U) is diagonalizable (in particular, if
right eigenvectors of F 0 (U) associated to a j , normalized so
that
6 AZEVEDO, MARCHESIN, PLOHR, AND ZUMBRUN
Proof. Equations (2.15) and (2.16) follow immediately from the denitions of g E and g BT .
Substituting Eq. (2.15), we obtain Eq. (2.17). Similarly, as B(U)
is a 2  2 matrix,
which yields Eq. (2.18) when substituted in Eq. (2.16).
Equations (2.9), (2.15), and (2.16) imply the following result.
Proposition 2.3. The functions g E and g BT are related by
Another consequence of Lemma 2.2 concerns E \ BT .
Proposition 2.4. If U 2 E is not an umbilic point, then the following are equivalent:
(a) U 2 BT ;
(c) the unique eigendirection of F 0 (U) is an eigendirection of B(U).
Proof. Statements (a) and (b) are equivalent by Eq. (2.21).
As the matrix F 0 (U) E (U)I is nilpotent; and because U is not an umbilic
point, this matrix is nonzero. Let r 6= 0 span its kernel. Since the square of this matrix is
zero, r also spans its range. In other words, F 0 (U)
If statement (a) holds, then with  denoting the common value of E (U) and BT (U ),
nilpotent and nonzero, and r belongs to its kernel. In fact,
since the square of this matrix is zero, r also spans its range. Let s be a vector such that
cr for some c, by nilpotency, and therefore
statement (c) is true.
Conversely, if statement (c) holds, trfB(U) 1 [F 0 (U) Hence
which implies statement (b).
Proposition 2.5. We have that
In particular, @MP  g.
Proof. By the strict parabolicity of B, 0 < tr Therefore
Eq. (2.18) shows that, given strict hyperbolicity at U , g BT (U) < 0 if and only if U 2 MP.
By Prop. 2.3, strict hyperbolicity (g E (U) < 0) follows from the condition g BT (U) < 0.
This result was proved for strictly hyperbolic points in Prop. 2.1 in Ref. [3]. Therefore
the new content concerns points of E . In this regard, it is useful to determine the behavior
of eigenvectors of F 0 (U) as U approaches a point U  2 E from within the strictly hyperbolic
region.
Lemma 2.6. Suppose that U  2 U has an open neighborhood V 0 containing no umbilic
points. Also assume that g E (U  )  0, and that DgE (U
has an open neighborhood V  V 0 such that the following statements hold.
(a) The sets VH := f U of hyperbolic
and strictly hyperbolic points in V are connected.
(b) In the hyperbolic region VH , the eigenvalues of F 0 are continuous functions a s  a f .
(c) Corresponding to the eigenvalues a j , for there exist: continuous right eigenvector
elds r j on VH , normalized so that jr 1, such that det(r s ; r f ) has a xed sign
continuous left eigenvector elds ' j in
normalized so that ' j r
(d) Suppose that U  2 E \V. If U ! U  from within V SH , then ' s (U)=j' s
Proof. By the hypothesis concerning g E and the Implicit Function Theorem, U  has an open
neighborhood statement (a) holds. We can also assume that if U
Statement (b) holds because F is C 1 . For each U 2 VH and the
matrix F 0 (U) a j (U)I has rank one (because there are no umbilic points in V 0 ); let R j (U)
denote its kernel. Then by statement (b), R j is a continuous line eld on VH . Moreover,
Let r
vectors lying in R j (U  ); if U  2 E , take r
f . By choosing
smaller, if necessary, we can assume that the angle between R j (U) and R j (U  ) is less
than =2 for all U 2 VH . Then for each U 2 VH be the unit vector
that has angle less than =2 relative to r
. Note that det(r s ; r f ) does not vanish in V SH ,
and hence has xed sign by connectivity. Therefore we can dene ' s (U) and ' f (U) to be
the rows of the inverse of the matrix (r s (U); r f (U)). The vectors r j and ' j so dened satisfy
statement (c). Statement (d) is a consequence of the normalization ' j r and the fact
that r
Remark. In the context of this lemma, consider a point U  2 E \ V, and suppose that
2 BT . By Prop. 2.4, the unique eigendirection of F 0 (U  ) is not an eigendirection of
. By
Prop. 2.3, g BT (U  ) > 0, so that g BT > 0 throughout an open neighborhood W  V of
U  . Hence, by formulae (2.17) and (2.18), ' s Br s and ' f Br f have xed and opposite signs
throughout W SH := W \V SH . In fact, by statement (d) of Lemma 2.6, these two quantities
have the signs of s  (r  ) ? B(U  )r  , respectively. Thus we see that, if s  (r  ) ? B(U  )r  is
positive (respectively, negative), W SH consists of points at which the Majda-Pego stability
condition is violated for the slow (resp., fast) family. Moreover, with
. Referring to Eq. (2.4), we see that this result re
ects
the transition from the second-order (in instability in the strictly hyperbolic region to the
rst-order instability in the strictly elliptic region.
2.2. Decomposition of the Majda-Pego unstable region. Let us dene
to be the Majda-Pego unstable region. Also
let
denote the elliptic region, slow-family instability region, and fast-family instability region,
respectively. Here, for a strictly hyperbolic point U , ' j and r j denote any left and right
eigenvectors associated to the eigenvalue a j of F 0 (U ), normalized so that ' j r
Proposition 2.7. The
sets
s ,
and
f form a partition
of
. In particular, the elliptic
region is contained in the Majda-Pego unstable region. Moreover, the boundaries
@
@
s ,
and
@
f are contained in E [ BT .
Proof.
and
by strict
parabolicity,
s
2.3,
If, on the other hand, g E (U) < 0,
then, by Eq. (2.18), g BT (U)  0 if and only if one of l s B(U)r s and ' f B(U)r f is nonpositive.
Thus
s
f
. The statement about boundaries is evident from the denitions of
the sets involved.
We now make the following nondegeneracy hypotheses on the model parameters of system
(1.1) for the domain U of interest:
(H1) at each point U 2 U such that g E
(H2) at each point U 2 U such that g BT
By the Implicit Function Theorem, the loci are
smooth curves. Also, E is the same as the elliptic boundary
@
and BT is the same as the
Majda-Pego boundary
Remark. A simple calculation shows that DgE is an umbilic point, i.e., F 0 (U)
is a multiple of the identity matrix. Therefore hypothesis (H1) excludes umbilic points.
In particular, we can restate one of the conclusions of Prop. 2.4 as: if g E
only if the unique eigendirection of F 0 (U) is an eigendirection of B(U ).
Let U be such that g E
also shows that D(gE g BT In
other words, E and BT are tangent where they meet. In view of this result, we are led to
adopt a higher-order nondegeneracy assumption:
(H3) at each point U 2 U such that g E
for all vectors V 6= 0 such that DgE (U)
Under this hypothesis, the points where E and BT intersect are isolated.
Proposition 2.8. Assume hypotheses (H1){(H3). Away from the discrete set E \ BT , the
boundaries
@
s and
@
f are smooth curves coinciding locally with either E or BT . Moreover,
@
s \
@
Proof. Consider a point U  2
@
s that does not lie in E \ BT . (The proof if U  2
@
f is
analogous.) Then, by Prop. 2.3 and Eq. (2.18), either (1) g E (U
Case (1). By assumption (H1), there exists an open neighborhood V of U  such that
g. As g BT (U  ) > 0, we can assume
CAPILLARY INSTABILITY IN MODELS FOR THREE-PHASE FLOW 9
that (ii) g BT > 0 throughout V. By Lemma 2.6, we can further assume that properties (a){
(c) hold as well. Properties (ii) and (a){(c), together with Eqs. (2.17) and (2.18), imply that
Therefore
so that
(@
Thus
@
s coincides with E in V.
Notice also that U  cannot lie in the
closure
because V is a neighborhood of U  that
does not
intersect
denition
comprises points of V SH for which ' f B r f  0,
Case (2). By assumption (H2), there exists an open neighborhood V of U  such that
g. As
we can assume that (ii) g E < 0 throughout V. By Lemma 2.6, we can further assume
that properties (a){(c) hold as well. Because B is strictly parabolic and ' s Br s  0 at
U  , we must have that ' f Br f > 0 at U  ; therefore we can assume that (iii) ' f Br f > 0
throughout V. Properties (ii), (iii), and (a){(c), together with Eqs. (2.17) and (2.18), imply
that
so that
(@
@
s coincides with
BT in V.
Finally, consider a point U  2
@
s \
@
f . If g E (U  ) were negative, then both '
s B(U  )r
s
and '
f B(U  )r
f would vanish, which is impossible because of the strict parabolicity of B.
Hence so that g BT (U  )  0 by Prop. 2.3. But it was just demonstrated that
cannot belong to
both
s
and
f . Therefore
Remark. As we traverse
@
the family for which the Majda-Pego condition (2.6)
fails can switch only at points common
to
s
and
f , which necessarily belong to the elliptic
boundary
@
To see how this switch can occur, let U  2 E \ BT and consider
traversing instead the elliptic boundary (which is tangent to BT ).
s be as in Lemma 2.6. By taking V smaller if necessary, we may
assume that E \ V and BT \ V are smooth curves meeting only at U  , where they have
rst-order, but not second-order, contact. Let U  2 E \ V. According to the discussion
in the remark following Lemma 2.6, if U
belongs to
@
s
(respectively,
@
according as s  (r  ) ? B(U  )r  is positive (resp., negative). Therefore
the family for which the Majda-Pego stability condition (2.6) fails switches if and only if
changes sign as U  passes U  . This change could be assured, for example, by
another nondegeneracy hypothesis.
These considerations show that, from the original perspective of traversing BT , the family
for which stability fails typically switches at elliptic points of BT . An example of this
phenomenon is illustrated in Fig. 3.1 below.
3. Existence of a DRS Point
We are now ready to establish a general result concerning existence of DRS points. We seek
a point at which conditions (D1){( e
are satised. By Prop. 2.5, it is sucient to identify
a point U  on the boundary of the Majda-Pego stable region MP for which condition (D2)
holds, i.e.,
where r  and '  are the eigenvectors of F 0 (U  ) determined by conditions (D1) and (D3).
Our strategy is to show that the directions of '  and r  rotate by an odd multiple of  as
U  traverses @MP . Since the left-hand side of Eq. (3.1) is homogeneous of odd degree in '
and r  , we can conclude that it vanishes at some point on @MP .
3.1. Degree. In the hyperbolic region f U g, the eigenvalues of the
are continuous functions a s (U)  a f (U ). Under hypothesis (H1), the matrix
F 0 (U) a j (U)I has rank one throughout the hyperbolic region, as there are no umbilic points
(see the remark after condition (H2)). Therefore the left and right kernels of F 0 (U) a j (U)I
are continuous left and right line elds L j and R j . Notice that for U
and R s
Let be a continuous, oriented closed curve in R 2 , and let R be a continuous line eld
on . Then the degree of R with respect to , denoted by deg(R; ), is 1= times the angle
through which R(U) rotates as U traverses once around (see Ref. [22, Chapter III, Sec. 1]).
For our purposes, we shall require a slightly more general notion of degree, dened for
boundaries of suitably nice sets. Toward this end, dene deg(R; ) as above even if is
not closed. (Of course, this number is not necessarily an integer.) If is decomposed as a
succession of continuous, oriented curves 1
(see Ref. [10, Chap. XVI, Sec. 1]).
Denote by C the class of bounded, nonempty sets W  R 2 such that @W consists of a
nite union of piecewise smooth curves, intersecting at a discrete set of points, that can be
oriented unambiguously by the outward-normal convention (i.e., at each point of @W that is
not one of the intersection points, W is locally dieomorphic to a half-plane). In particular,
by Props. 2.7 and 2.8, MP,
s
and
s
f belong to C.
If W 2 C and R is a continuous line eld on @W, we can decompose @W as a succession
of continuous, oriented curves 1
deg(R; @W) :=
This number is independent of the decomposition. Moreover, it is an integer because any
decomposition of @W into continuous, oriented curves can be grouped into a union of closed
curves by forming the boundaries of the connected components of the set obtained from W
by removing the nitely many intersection points.
The following useful properties of the degree hold.
(P1) The degree deg(R; ) is invariant under homotopies of both and R
disjoint members of C, and suppose that W :=
belongs to C. If R is a continuous line eld on @W 1 , @W 2 and @W, then deg(R;
R is a continuous line eld on W , then deg(R;
For property (P1), see, e.g., , Ref. [42, Sec. 34]. Property (P2) holds because @W 1 and
consist of @W together with a nite number of pairs of coinciding curves with opposite
orientations. Property (P3) is proved for a simply connected region W by shrinking @W to a
point and invoking homotopy invariance; it is proved in general by subdividing the compact
set W into simply connected regions.
3.2. Main result. We can now state the main result of this section.
Theorem 3.1. Let hypotheses (H1){(H3) hold on a bounded region U , and suppose that
is a (nontrivial) piecewise smooth simple closed curve lying within MP. If deg(R
, then there exists a DRS point strictly inside .
We will prove this theorem using a series of lemmas. Let  denote the closed bounded
region with boundary . Then  belongs to the class C. Without loss of generality, we may
assume that the Majda-Pego unstable
region
contained in the interior of
. Notice
that
is nonempty; otherwise  would be contained in the strictly hyperbolic
region and R k would be a continuous line eld on , so that deg(R
would vanish by property (P3).
Lemma 3.2. Under the hypotheses of Theorem 3.1, deg(R
and the degrees deg(R s ; @
, and deg(L f ; @
all equal this common
value.
Proof. Since, for each U in the strictly hyperbolic region, R s (U) and R f (U) are linearly
independent directions in the plane, R s and R f have equal degrees around . Likewise, by
orthogonality of left and right eigenvectors, deg(L
For the line eld R j is continuous throughout
because this set is contained
in the (non-strictly) hyperbolic region. Further,
belongs to class C by the assumed
smoothness of
@
and its separation from = @. Therefore
Similarly deg(L
@On
s
f , dene the new line elds
e
R s
in
in
and e
in
in
The line eld e
R is well-dened and continuous
on
because, by Prop. 2.8, R s and
on
s
Similarly, e
L is well-dened and continuous. Moreover, for
@
e
Lemma 3.3. Under the hypotheses of Theorem 3.1, deg( e
and deg( e
Proof. By properties (P2) and (P3) of the degree,
@
@
@
@
R. But by Eq. (3.6), deg(R;
@
independent of the choice R
R f , or e
R. Therefore Lemma 3.2 implies that deg( e
Similarly, deg( e
@ equals this value.
Proof of Theorem 3.1. By hypothesis (H1), conditions (D1) and ( e
are satised on
@
R and '  2 e
L. It remains only to show that for some U  2 @
On the hyperbolic region, in particular on BT , the continuous line elds e
R and e
locally continuous, nonzero vector elds r  and '  ; indeed, each line eld locally generates
precisely two branches. Choose a connected component of
@
around which e
R and e
L have
odd degree, and choose a particular branch for r  and '  at some initial point on this
component. Following this branch while traversing once around this component, we nd by
Lemma 3.3 that the quantity '  F 00 (U  )(r  ; r  ) must change sign, since both e
L and e
R rotate
by an odd multiple of . Thus it vanishes somewhere on @
Remark. Arguments with a similar topological
avor have been used to establish the existence
of points at which strict hyperbolicity fails [37, 38]. In [20, 21, 41, 30], strict hyperbolicity
failure was established with dierent arguments.
Remark. The DRS point generically lies in the region of strict hyperbolicity. If it were to
occur on E , Prop. 2.4 indicates that r  would have to be an eigendirection of B(U  ).
Example 3.4. For generic two component models with quadratic
ux F and constant matrix
B, the elliptic boundary and the BT locus intersect in precisely two points, unless B is
a multiple of the identity matrix, in which case they coincide. Indeed, one can show that
the eigenvector r  undergoes a rotation of angle  as U  traverses the elliptic boundary
(it suces to compute this rotation for homogeneous quadratic models). Thus the direction
of r  coincides precisely once with each of the eigendirections of B. If B has more than
two eigendirections, then B is a multiple of the identity matrix, in which case E and BT
coincide. Otherwise B has only two eigendirections, and there are precisely two points in
may verify that the family in which Majda-Pego stability fails switches at
these two points, illustrating the conclusions of the remark at the end of Sec. 2.2.

Figure

3.1 shows MP and the decomposition
of
s
f for a model such
that
is compact. The points in E \ BT are indicated by two small open circles, the three
lled circles represent DRS points, and the rectangles mark line eld singularities [31, 23].
The curves through the rectangles constitute the in
ection locus I, along which genuine
nonlinearity fails in the family f or s indicated by whether the curve is solid or dashed.
4. Three-Phase Flow
In this section, we introduce a class of systems of the type given by Eq. (1.1), which are
used in Petroleum Engineering to model three-phase
ow in porous media. We discuss the
properties of the model that have a role in determining features such as Majda-Pego stability
regions and DRS points.
One of the main consequences of the analysis of important classes of three-phase
ow
models in Sec. 5, under the generic assumptions that such models satisfy (H1){(H3), follows
from Theorem 3.1. It is the following.
Theorem 4.1. Either the Majda-Pego instability
region
accumulates at the boundary of
the saturation triangle @ or else there exists a DRS point interior to .
MP
s
f
I

Figure

3.1. The Majda-Pego stable region MP , the decomposition of the
Majda-Pego unstable
region
s
f , and the coincidence, Bogdanov-
Takens, and in
ection loci E , BT , and I for a quadratic model.
4.1. The basic equations. We consider one-dimensional, horizontal
ow of three immis-
cible
uid phases in a porous medium [32]. For concreteness, we consider a
uid composed
of gas, oil and water, mixed at macroscopic level. The dierences among these phases lie
in some
ow properties. We assume that the whole pore space is occupied by the
uid and
that there are no sources or sinks. Compressibility, thermal and gravitational eects are
considered to be negligible.
The equations expressing conservation of mass of water, gas, and oil are
@
@t
@x
respectively, where  denotes the porosity of the porous medium. For the phase i, s i denote
the saturation,  i the density and v i is the seepage velocity (the product of the saturation
by the particle velocity of the phase i). Since the
uid occupies the whole pore space, the
saturations satisfy
As a consequence, any pair of saturations in the saturation triangle  may be chosen to
describe the state of the
uid.
The theory of multiphase
ow in porous media is based on the following form of Darcy's
law of force [32, 35, 4]:
@
@x
where K denotes absolute permeability of the porous medium,  i  0 is the mobility of
phase i, and p i is the pressure of phase i. The mobility is usually expressed as
the ratio of the relative permeability k i and the viscosity  i of phase i.
14 AZEVEDO, MARCHESIN, PLOHR, AND ZUMBRUN
The porosity  and absolute permeability K are associated to the rock; we take them to be
constant. Neglecting thermal eects and compressibility,  i and  i are constant too, and we
can rewrite Eqs. (4.1) without  i . Each relative permeability k i depends on the saturations.
Experimentally, k i increases when s i increases, and the relative permeabilities never vanish
simultaneously.
Let us denote the dierence between the pressures in phases i and j (i 6= j) by
This pressure dierence, called the capillary pressure, is measured experimentally as a function
of the saturations. Dene the total mobility  and the fractional
ow functions f i
by
Of course,
1. Introducing the total seepage velocity
using algebraic
manipulation, we can write that
@
@x
Adding and subtracting v i in the last equation and noting that @
we see that
@
@x
Therefore, by Eqs. (4.1) and (4.7), the equations governing the
ow are
@
@x Kw [f g
@
@
wo
@
@x K g [f w
@
@
@
@t
@x
@
ow
@
Summing Eqs. (4.1), we nd @
so that v is a function of t alone. (This simplication
occurs only for
ow in one spatial dimension. In general, we obtain an elliptic equation for the
pressure, with coecients depending on the saturations.) Assuming that v never vanishes, it
is possible to change the variable t so that v is constant. We do not consider the case v(t)  0
because in this case the system (4.8) does not contain the terms related to transport. As v
is nonzero, we can set removing v, , and K from
system (4.8). For simplicity of notation, we drop the tildes.
Of course any one of the equations in the three-component system (4.8) is redundant and
the system can be reduced to a two-component system. As a result of the redundancy in
the system (4.8), one of its characteristic speeds is 0. The two other characteristic speeds
are the same for subsystem of two equations; see Ref. [30].
We will nd it convenient later to be
exible in our choice of the two saturations entering
this two-component system. Hence we will use u 1 and u 2 to denote two of the saturations
in the reduced system, where the other saturation u 3 is replaced by 1
Also it is useful to replace p 12 by . Thus we obtain the two-component system for
three-phase
ow
@
@t
@
@
@t
@
This system can be written in compact form as Eq. (1.1) with
and
and
The quantities U , F and B(U) are called the state vector,
ux function and diusion matrix,
respectively. We refer to Q(U) and P 0 (U) as balance matrix and capillary pressure Jacobian.
There are two major assumptions that we make concerning the mobilities. The rst is
that
(More generally,  i vanishes when s i is below the irreducible saturation for phase i, which,
for simplicity, we take to be zero.) Because  the edge
for system (4.8). Indeed, system (4.9) reduces to the scalar conservation law
@
@t
on the edge and a similar reduction occurs on the other edges. In other words,
three-phase
ow equations reduce to the Buckley-Leverett equation for two-phase
ow.
There are two characteristic speeds for system (4.9), one of which reduces to the Buckley-
characteristic speed on the edges. Our second assumption is that the other characteristic
speed is positive and strictly smaller than the Buckley-Leverett characteristic speed
near the edges. For example, on the edge
Thus we see that the characteristic speeds are f 1;1 and f 2;2 , that f 1;1 is the Buckley-Leverett
speed, and that the second assumption can be written as
4.2. Assumptions on the fractional
ow functions. In this section we present assumptions
used in the rest of the paper.
For brevity we use a subscript \; j" to denote the partial derivative with respect to u j for
2. The Jacobian matrix of F is therefore
According to Eq. (4.5),
2: (4.17)
We shall make the following assumptions concerning the fractional
ow functions: f 1 , f 2 ,
and f 3 are continuously dierentiable functions on the closed saturation triangle  such that
on the open edge
f 1f 2;1 has a continuous extension to the open edge where
in the interior of  near
are bounded away from zero in the interior of  near
Each of the assumptions above actually represents three dierent assumptions, corresponding
to a choice of two saturations u 1 and u 2 among s w , s g and s
The rst part of assumption (4.18) corresponds to the requirement that the permeabilities
be nonnegative. By the second part of this assumption, f on the edge
so that the Jacobian F 0 is upper triangular and its eigenvalues are f 1;1 and f 2;2 . Hence
assumption (4.19) states that the eigenvalues are distinct on each open edge, and that the
fast eigenvector is parallel to the edge. Thus the three-phase
ow reduces to two-phase
ow
along each edge, with the fast eigenvalue corresponding to the Buckley-Leverett speed (see
Eq. (4.13)). In particular, f at the corner re
ecting
the immiscibility of two-phase
ow. Assumption (4.21) implies that the eigenvalues are
distinct in the interior of  near each corner. (Necessarily, the eigenvalues coincide at each
corner. Thus, the model is strictly hyperbolic near the boundary except at the corners.) This
assumption also implies that the saturations u 1 and u 2 vary in opposite directions across slow
rarefaction waves near each corner. Assumptions (4.20) and (4.22) simplify the conditions
for stability at the boundary of . We will show that these assumptions hold for certain
models employed in petroleum reservoir engineering in Sec. 6.
4.3. Assumptions on the capillary pressures. In principle, the capillary pressures p ow
and p og are experimentally measured functions of all saturations. (The third capillary pressure
is obtained as In general, the capillary pressure Jacobian (4.11)
is
We make the following assumptions:
positive denite and diagonally dominant inside ; (4.24)
on the closed edge
Assumption (4.25), which guarantees that the two-phase
ow equation (4.13) for the edge
has a positive diusion coecient, actually represents three dierent assumptions,
corresponding to a choice of two saturations u 1 and u 2 among s w , s g and s
4.4. Well-posedness. Next, we investigate the eects of capillary pressure in system (1.1),
as represented by the diusive term (B(U)U x ) x . We rst verify that this term is strictly
parabolic (i.e., the eigenvalues of B(U) have strictly positive real parts) in the interior of
the saturation triangle .
Lemma 4.2. Under assumption (4.18), the balance matrix Q(U) is symmetric and positive
denite in the interior of the saturation triangle , and det
Proof. Rewriting
it is apparent that Q(U) is symmetric and diagonally dominant, with nonnegative diagonal
entries that are positive in the interior of , by assumption (4.18). One can verify that
det so that det
Lemma 4.3. If P 0 satises assumptions (4.24) and (4.25), then the eigenvalues of the diffusion
matrix B(U) have positive real part in the interior of , whereas on the boundary of
, one eigenvalue of B(U) is zero and the other is positive.
Proof. A real matrix has eigenvalues with positive real part provided that its trace and
determinant are positive. The determinant of both Q and P 0 are non-negative, by diagonal
dominance (4.24), so that det strict inequality in the
interior of . Consulting Eq. (4.26), we nd that
by assumptions (4.24) and (4.25).
Remark. It is easily checked (see the calculations below) that the normal to  is a left
eigenvector of F 0 (U) and a left null vector of B(U) on @. This is a necessary condition in
order that  be invariant under the evolution of system (1.1), see Theorem 1.47, p. 200 of
[40]. However, suciency does not follow from this standard result, due to degeneracy of
B(U) on the boundary; indeed, the question of invariance in this degenerate situation seems
to be an interesting one in its own right.
5. Linear Stability Conditions and Degree
In this section we compute the degree for the fast-family line eld for three-phase
ow
models. We also nd conditions under which the Majda-Pego stability conditions (2.6) hold
near @.
5.1. Eigenvalues and eigenvectors. It proves convenient to work with the scaled Jacobian
matrix
and to introduce the following notation:
Let v s and v f denote the slow and fast eigenvalues of A, respectively. Then
a  d b
c a  d
We may choose the matrix of right eigenvectors to be
This a good choice for calculations near the edge where The
corresponding matrix of left eigenvectors is given by
d)
a
Using the relation d bc, we nd that the determinant of R in Eq. (5.5) 2d[a
On the other hand, for calculations near the corner where
choice for the right eigenvector matrix is
(The reason is that the projections (1; 1) r s and (1; 1) r f are positive.) The corresponding
matrix of left eigenvectors is
d)
Again by the relation d bc, the determinant of R in Eq. (5.7) is 4d[d
5.2. Degree. In this section we show that, for models satisfying the assumptions of Sec. 4.2,
the fast-family line eld has winding number 1 around a suitable contour close to the
boundary of the saturation triangle.
Lemma 5.1. For a model satisfying assumptions (4.18), (4.19), and (4.21), consider the
arc of a circle centered at leading from the edge to the edge
If the radius of the arc is suciently small, the fast-family line eld turns by angle =2
upon traversing the arc.
Proof. The projection of the eigenvector r f , as dened in Eq. (5.7), onto the vector (1; 1) is
assumption (4.21), this quantity is positive along any open arc of suciently
small radius. At the edge assumption (4.19) and assumption (4.18),
so that d = a and r which is a positive multiple of (1; . Similarly, at the
edge so that d = a and r which is a positive
multiple of (0; 1) T . Therefore r f turns by angle =2.
Denition 5.2. Let the contour  be the boundary of the set   of points in  at least
distance  from the boundary.
Proposition 5.3. Consider a model satisfying assumptions (4.18), (4.19), and (4.21). For
suciently small  > 0, the degree of the fast-family line eld around  is 1.
Proof. According to the preceding lemma, the fast-family line eld rotates by =2 at the
corner Taking account the changes of coordinates, the same lemma shows
that the fast-family line eld rotates by =4 at each of the two corners
traversing edges, the fast-family line eld does not rotate because it
is parallel to the edge. Therefore the line eld rotates by  around  .
CAPILLARY INSTABILITY IN MODELS FOR THREE-PHASE FLOW 19
5.3. Stability analysis near edges. In this section we determine a necessary and sucient
condition for a model to be linearly stable near the open edges of the saturation triangle,
in that conditions (2.6) hold with strict inequality except on the edge. Without loss of
generality, we focus on the edge The calculations presented in this section are not
applicable at the corners, so we assume that 0 <
We choose the eigenvector matrices (5.5) and (5.6). If u assumption
(4.18) and a > 0 by assumption (4.19). In particular, d = a and det Therefore,
when evaluated along the edge, the fast family eigenvectors r f and ' f are positive multiples
of (1; 0) T and (1; b=(2a)), respectively. Assumption (4.18) also entails that, when
the only nonzero entry of the balance matrix Q is Q
sequently, by denition (4.10) of the diusion matrix B, ' f B r f is a positive multiple of
along the edge We have proved the following result.
Proposition 5.4. For a model satisfying assumptions (4.18), (4.19), and (4.25), the fast-
family Majda-Pego condition holds in an open neighborhood of the open edge
The right eigenvector for the slow family is r However, we shall
see that ' s hence to verify that ' s B r s  0, we cannot simply set
As ' s is a positive multiple of (c; (a +d)), Eq. (4.11) implies that ' s Q is a positive multiple
of the vector d)). By assumption (4.20),
we may extend  1
Q, to be a continuous function at
positive multiple of ( 1
It proves useful to write this last vector as r T
nd a simplied expression for
First we note that when u
Using these results in the denition (5.9) of  , we nd that
(The quantity  3;1  3;2 is the derivative of  3 with respect to u 1 with u 3 xed.) We can
now formulate a condition equivalent to the Majda-Pego condition for the slow family.
Proposition 5.5. For a model satisfying assumptions (4.18), (4.19), (4.20), and (4.25), let
let  be given by Eq. (5.13). Then
r T
is a necessary and sucient condition for the model to satisfy the slow-family Majda-Pego
condition in an open neighborhood of the open edge
strict inequality when u
5.4. Stability analysis at corners. In this section we determine a necessary and sucient
condition for a model to be linearly stable near the corners of the saturation triangle, in that
conditions (2.6) hold with strict inequality except on the edges. Without loss of generality,
we focus on the corner
Condition (2.5) is equivalent to the nonnegativity of the diagonal elements of LQP 0 R.
For the eigenvector matrices R and L we choose those dened by Eqs. (5.7) and (5.8),
respectively. Using the balance matrix dened in Eq. (4.11), we see that
(det R)' s
(det R)' f
By assumption (4.21), det positive in the interior of  in a neighborhood
of the corner.
Proposition 5.6. For a model satisfying assumptions (4.18), (4.22), and (4.25), the fast-
family Majda-Pego condition holds in the interior of  in a neighborhood of the
corner only if
in the interior of  in a neighborhood of this corner, where a, b, c, and d are given by
Eqs. (5.2) and (5.3).
Proof. The stability condition for the fast family is that (det R)' f QP 0 r f > 0. By assumption
(4.18) (and the corresponding assumption along the edge tend to
zero at the corner. In particular, we may neglect f 1 and f 2 relative to 1 in the expression
for (det R)' f Q. Moreover, by assumptions (4.18) and (4.19), a, b, c, and hence d vanish at
the corner; therefore, by assumption (4.22), f 1 (a d) is negligible relative to
d) is negligible relative to . Thus (det R)' f Q can be approximated
by d) T , the inequality
(det R)' f QP 0 r f > 0 is equivalent to  f > 0 near the corner.
By an analogous argument, we can prove the following result.
Proposition 5.7. For a model satisfying assumptions (4.18), (4.22), and (4.25), the Majda-
Pego condition for the slow family holds in the interior of  in a neighborhood
of the corner only if
d)
in the interior of  in a neighborhood of this corner, where a, b, c, and d are given by
Eqs. (5.2) and (5.3).
5.5.

Summary

. Combining the results of this section, we have the following theorem.
Theorem 5.8. Consider a model satisfying assumptions (4.18){(4.22) and (4.25) along with
conditions (5.14), (5.17), and (5.18). For suciently small  > 0, the contour  of Deni-
tion 5.2 lies within the interior of MP and the degree of the fast-family line eld around
is 1.
6. A Three-Phase Flow Model
In this section, we present a model for the relative permeabilities and the capillary pressures
used in petroleum reservoir engineering and verify the assumptions of Secs. 4.2 and 4.3.
Using the results of Sec. 5, we show that the model is linearly stable near the boundary of
the saturation triangle except in an open set containing two of the corners in its closure.
6.1. Stone's Permeability model. The permeability of each of the three phases is an
experimentally measured function that depends on the saturations of the phases. Leverett
and Lewis [28] found that, for many reservoirs, kw depends only on s w and k g depends only
on s g . We make the same assumption here.
For most reservoirs, the oil permeability depends on two saturations. As there is little
experimental data concerning the three-phase region, this function is usually measured only
along two edges of the saturation triangle and then interpolated to the interior of . Restricted
to the s reduces to the relative oil-water permeability, k ow ; similarly,
edge is the relative oil-gas permeability k og . To satisfy assumption (4.12), we
require that k ow (s w ) and k og (s g ) are non-negative and that k ow (s w
proposed an interpolation
scheme for determining k throughout  from k ow and k og :
ow (s w )k og (s g )
The particular model we examine is the quadratic Stone's model, which uses kw (s w
ow (s w so that
w , and k
In the notation of system (4.9), the mobilities of the quadratic Stone's model are
if we associate u 1 to s w and u 2 to s g , whereas
if we associate u 1 to s w or s g and u 2 to s
6.2. Leverett's Capillary pressure model. Based on the work of Leverett [28], Aziz
and Settari [4] propose that the the capillary pressure functions for three-phase
ow can be
taken to have the form p wo ow ow and P og are certain
monotone decreasing functions. Denoting
, we assume that
Associating u 1 to s w and u 2 to s g in system (4.9), p 13 and p 23 represent wo and p go ,
respectively. Therefore the Leverett capillary pressure Jacobian is
On the other hand, considering the oil to be one of the phases in Eq. (4.9), the Leverett
capillary pressure Jacobian takes a dierent form. Associating u 1 to s w and u 2 to s
22 AZEVEDO, MARCHESIN, PLOHR, AND ZUMBRUN
loss of generality, we have that u 3 is s g , p 13 is p . In this case,
It is easy to see that assumptions (4.24) and (4.25) of Sec. 4.3 hold for the Leverett
capillary pressures under assumption (6.5).
6.3. Verication of assumptions. We now verify the assumptions of Sec. 4.2.
Lemma 6.1. The quadratic Stone's model satises assumptions (4.18){(4.22) if we associate
u 1 to s w or s g and u 2 to s that the oil viscosity is no less than half of the
water and gas viscosities, i.e.,2 w
Proof. The mobilities are given by Eqs. (6.4). Assumption (4.18) is evidently true.
Let us verify assumption (4.19). Since
we see from Eq. (4.17) that
is nonnegative everywhere and reduces to
u2=0
on the edge On the other hand,
so that
reduces to
u2=0
In particular, f at the corner on the open edge
only if
These inequalities hold by virtue of assumption (6.8).
Now we verify assumptions (4.20), (4.21) and (4.22). Since
dividing by  2 gives
so that f 1
f 2;1 has a continuous extension to to the open edge where
In
f 2;1 has a positive limit at (0; 0). Moreover,
so that f 1
bounded away from zero in the interior of  in a neighborhood of
(0; 0). Similarly,
so that f 1;2 > 0 in the interior of  in a neighborhood of (0; 0). In fact, f 1
f 1;2 has a positive
limit at (0; 0).
Lemma 6.2. The quadratic Stone's model satises assumptions (4.18){(4.22) if we associate
u 1 to s w and u 2 to s g .
Proof. The mobilities are given by Eqs. (6.3). Assumption (4.18) is evidently true.
Since
we see by Eq. (6.10) that f 2;2 is nonnegative near is zero on the
edge On the other hand,
so that by Eq. (6.13), f 1;1 on u by Eq. (6.14). Therefore f 1;1 > f 2;2 on
the open edge assumption (4.19) holds.
shows that f 1
f 2;1 has a continuous extension to the edge
and tends to 2 at the corner (0; 0). Similarly, by Eq. (6.20), f 1
tends to 2 at the corner.
Thus assumptions (4.20), (4.21), and (4.22) hold.
6.4. Linear stability near edges. In this section we prove the following result.
Theorem 6.3. Assuming conditions (6.5) and (6.8), the quadratic Stone's model with the
capillary pressures satises the stability conditions (2.5) near each open edge of ,
with strict inequality away from the edges.
Proof. By Props. (5.4) and (5.5) and the results of Secs. 6.2 and 6.3, the proof reduces to
verifying condition (5.14) of Prop. (5.5).
First, we associate u 1 to s w or s g and u 2 to s o . By Eqs. (6.14), (6.11), and (6.21) the
vector r reduces to
on Eqs. (6.12) and (6.12),  1
so that
Using the capillarity matrix P 0 given in Eq. (6.7), we nd
r T

Therefore assumptions (6.5) and (6.8) imply condition (5.14).
Second, we associate u 1 to s w and u 2 to s g . By the calculations in the proof of Lemma 6.2,
and
Using the capillary pressure Jacobian (6.6), we nd that
Thus assumption (6.5) is implies condition (5.14).
6.5. Linear stability and instability at corners. In this section we investigate Majda-
Pego conditions at the corners of the saturation triangle .
Theorem 6.4. Consider a model satisfying assumptions (4.18), (4.22), and (4.25) and using
the Leverett capillary pressures. This model satises the Majda-Pego condition (2.6) for the
fast family in the interior of  in neighborhoods of each corner; it also satises the Majda-
Pego condition for the slow family in the interior of  in neighborhood of the corner s
Proof. First associate u 1 to s w and u 2 to s g . Then Leverett's capillarity matrix P 0 is given
by Eq. (6.6) and inequalities (5.18) and (5.17) become
d)
d) > 0; (6.31)
d)
d) > 0: (6.32)
In the interior of  in a neighborhood of the corner, b and c are positive and d > jaj, so that
the coecients of
and  in  s
wg and  f
wg are positive; thus  s
wg > 0.
Now associate u 1 to s w or s g and u 2 to s
by Eq. (6.7), and the inequality (5.17) becomes
d)
c) > 0:
By the same argument as given above, the condition  f
wo > 0 is satised in the interior of
in a neighborhood of the corner.
However, the Majda-Pego condition for the slow family is violated for the quadratic Stone's
model and Leverett capillary pressures, as we now show.
Theorem 6.5. Consider the quadratic Stone's model with the Leverett capillary pressures.
Assume that the oil viscosity is no less than half of the water and gas viscosities. This model
in an open set containing the corners s
in its closure.
Proof. Without loss of generality, we may associate u 1 to s w and u 2 to s Consider the set
where Eqs. (6.13), (6.10), (6.12), and (6.9),
so that in a neighborhood of the corner, the set a = 0 is a curve through the corner tangent
to the line
When
c ]((
(det R)' s
as seen from Eqs. (5.7) and (5.15). By Eqs. (6.21) and (6.19),
so that
c) as along the curve a = 0. Therefore
near the corner.
7. Conclusion
Linear stability is a desirable property. We have numerical evidence that certain three-phase
models satisfying assumptions (4.18){(4.22) are linearly stable near the boundary of
the saturation triangle and a curve  satisfying the hypotheses of Theorem 4.1 can be
constructed. This is not the case for quadratic Stone's model with Leverett's capillary
pressures: as we can see from Theorem 6.5, there is an instability region near the two
vertices. On the arc with radius 0:04, the instability region corresponds to the negative
part of the lowest curve in Fig. 7.1. However this gure shows that for larger radii, the
arcs avoid the instability region. We have veried numerically that the eigenvectors wind
correctly along these arcs. Based on these facts, we have constructed curves  which also
satisfy the hypotheses of Theorem 4.1. Also for the quadratic Stone's model we found
numerically regions in parameter space for which the hypotheses (H1){(H3) are satised.
Theorem 4.1 establishes the existence of a DRS point inside  , which generically implies
other instabilities. (Indeed, the above discussion shows that a DRS point is to be expected
when the instability regions near the vertices are disconnected from the elliptic region.)
For some nongeneric values of parameters, the elliptic region may collapse into umbilic
points that coincide with the DRS points. The same collapse occurs for certain other three-
phase
ow models. We have no evidence of DRS-instabilities for such degenerate DRS points.
We regard the stability inequalities as providing mathematical guidelines for the construction
of adequate three-phase
ow models. Three-phase
ow uses quantities such as
26 AZEVEDO, MARCHESIN, PLOHR, AND ZUMBRUN2610140
f

Figure

7.1. Plot of ' s B r s =r 3 (multiplied by 10 3 det R) versus polar angle
for 0:08. For this gure, the viscosities are
the values of the capillary pressure derivatives,
evaluated at the corner, are
permeabilities and capillary pressures that can only be measured accurately at the boundary
of the saturation triangle. We may regard the stability conditions as constraints on
interpolation schemes for this data into the interior of state space.

Acknowledgments

We thank Beata Gundelach for careful assistance in preparing the manuscript. We also
thank Somsak Orankitjaroen for thoroughly reading the manuscript and locating errors in a
preliminary version.



--R

Multiple viscous solutions for systems of conservation laws
Nonuniqueness of nonclassical solutions of Riemann problems

Elsevier Applied Science
Conservation laws of mixed type describing three-phase ow in porous media



Cubic Li
Features of three-component

Existence and asymptotic behavior of measure-valued solutions for three-phase ow in porous media

Oscillation waves in Riemann problems inside elliptic regions for conservation laws of mixed type


On the Riemann problem for a prototype of a mixed type conservation law
On the strict hyperbolicity of the Buckley-Leverett equations for three-phase ow in a porous medium

A global formalism for nonlinear waves in conservation laws




Steady ow of gas-oil-water mixtures through unconsolidated sands
Stable viscosity matrices for system of conservation laws
Stable hyperbolic singularities for three-phase ow models in oil reservoir simulation

Fundamentals of numerical reservoir simulation


The physics of ow through porous media
The Riemann problem for a class of conservation laws of mixed type

Loss of real characteristics for models of three-phase ow in a porous medium
Admissibility criteria for propagating phase boundaries in a van der waals uid
Shock waves and reaction-diusion equations

Departamento de Matem
--TR
Conservation laws of mixed type describing three-phase flow in porous media
Loss of strict hyperbolicity of the Buckley-Leverett equations for three phase flow in a porous medium
On the strict hyperbolicity of the Buckley-Leverett equations for three-phase flow in a porous medium
Admissibility conditions for shocks in conservation laws that change type
Oscillation waves in Riemann problems inside elliptic regions for conservation laws of mixed type
Nonuniqueness of solutions of Riemann problems
On the oscillatory solutions in hyperbolic conservation laws
Fundamentals of Numerical Reservoir Simulation
Nonexistence of Riemann solutions for a quadratic model deriving from petroleum engineering

--CTR
E. Abreu , J. Douglas, Jr. , F. Furtado , D. Marchesin , F. Pereira, Three-phase immiscible displacement in heterogeneous petroleum reservoirs, Mathematics and Computers in Simulation, v.73 n.1, p.2-20, 6 November 2006
capillary pressure instability;nonunique Riemann solution;flow in porous media