Stone-Marchesin Model Equations of Three-Phase Flow in Oil Reservoir - - PDF document

Stone-Marchesin Model Equations of Three-Phase Flow in Oil Reservoir Simulation

Fumioki ASAKURA, Osaka Electro-Communication Univ. JAPAN asakura@isc.osakac.ac.jp

1

• 1. Introduction

WAG Injector Oil Producer Oil Water Gas Water-Gas Shale

Figure 1: WAG Enhanced Oil Recovery (schematic picture) 2

Primary Recovery By the underground pressure, usually about 20% of the oil in an oil reservoir can be extracted. Secondary recovery By injecting water and gas (air or CO2), generally 25% to 35% of the oil can be extracted. Water-Alternating-Gas (WAG) Enhanced Oil Re- covery: Water Injection Good sweep efficiency, but 40 to 60% of the original oil on-site is left behind. Gas Injection Good displacement efficiency, but an expen- sive operation. WAG Injection More efficient than injection of water or gas alone.

3

Overview: Marchesin D. & Plohr B. (2001), Theory of Three- Phase Flow Applied to Water-Alternating-Gas Enhanced Oil Recovery, Proceedings of the 8th International Conference in Magdeburg, Vol.II, Birkh¨ auser Verlag, 693–702. Plan of this Talk:

• Model equations
• Hyperbolicity, elliptic region
• Geometry of characteristic field, 2-phase like flow curves
• Compressive, under and overcompressive shock waves
• Geometry of Hugoniot curves
• Entropy functions

4

Stone’s Model [3]: Neglect the gravity. Assume: the medium is homogeneous and the flow is incompressible and immiscible. water gas

• il

Volume Fractions: sw = u sg = v so = 1 − u − v Permeability Functions: kw kg ko Fluid Viscosity: µw µg µo Fluid Velocity: vw vg vo Pressure: pw pg po Capillary pressure: pc = pnon wetting phase − pwetting phase. Water-oil interface: water is the wetting phase. Gas-oil interface: oil is the wetting phase. pow = po − pw, pgo = pg − po.

5

Leverett’s assumption: pow is a decreasing function of u = sw, and pgo is an increasing function of v = sg. Darcy’s Law vi = −ki µi ∇pi, i = w, g, o. Mass Conservation Laws: ∂si ∂t + ∇ · vi = 0, i = w, g, o. By eliminating ∂po

∂x and denoting

D =

• j=w, g, o

kj µj

6

       ∂sw ∂t + ∂ ∂x kw µwD

• =

∂ ∂x kw µw

• 1 − kw

µwD ∂pwo ∂x − kg µgD ∂pgo ∂x

• ,

∂sg ∂t + ∂ ∂x kg µgD

• =

∂ ∂x kg µg

• − kw

µwD ∂pwo ∂x +

• 1 − kg

µgD ∂pgo ∂x

• Stone’s assumption: The water and gas permeability

functions depend only on the water and gas volume fraction kw = kw(u), kg = kg(v). If the capillary pressure is negligible, by using relative perme- ability functions f(u) = kw(u) µw , g(v) = kg(v) µg h(u, v) = ko(u, v) µo

7

Water: ∂u ∂t + ∂ ∂x

• f(u)

f(u) + g(v) + h(u, v)

• = 0,

(1) Gas: ∂v ∂t + ∂ ∂x

• g(v)

f(u) + g(v) + h(u, v)

• = 0

(2) for (u, v) ∈ Ω : 0 < u + v < 1, u, v > 0. If w = 0, 2-phase flow is governed by the Buckley-Leverett equation ∂u ∂t + ∂ ∂x

• f(u)

f(u) + g(1 − u)

• = 0.

N.B. 1 An example of a single conservation law with a non-convex flux function.

8

• 2. Geometry of Characteristic Vector Field

Hyperbolicity: Denote D = f(u) + g(u) + h(u, v), F(U) = t

f D, g D

• .

Hyperbolic F ′(U) has real eigenvalues λ1(U), λ2(U) for any U ∈ Ω. Strictly Hyperbolic Eigenvalues are distinct: λ1(U) < λ2(U). Umbilical Point U ∗ λ1(U ∗) = λ2(U ∗) and F ′(U) is diago- nalizable, hence a scalar matrix. The eigenvalue λj(U) is called the jth characteristic speed and corresponding right eigenvector Rj(U) is called the jth characteristic vector field.

9

2-Phase Like Flow Curves (Medeiros [16]): F ′(U) = 1 D2 f ′(g + h) − fhu −f(g′ + hv) −g(f ′ + hu) g′(f + h) − ghv

• Ru =

1

• , Rv =

1

• , Rw =
• 1

−1

• , respectively, are

characteristic vectors at v = 0, u = 0, w = 0, respectively. 2-phase like flow curves L1, L2, L3 are defined by L1 : (f ′ − g′)h + (f + g)(hu − hv) = 0, L2 : g′ + hv = 0, L3 : f ′ + hu = 0. where Rw, Rv, Ru, respectively are characteristic vectors on L1, L2, L3, respectively.

10

U*

u + v = 1 v = 0 u = 0

L2 L1 L3

R R R R R Rw Rv Ru

Figure 2: 2-Phase Like Curves (Quadratic Marchesin’s Model) 11

Introduce the 2-phase like variable: ξ = (f ′−g′)h−(f +g)(hu−hv), η = g′+hv, ζ = f ′+hu. Lemma 1 The discriminant has the form: Dchar = 1 D4

• {f ′(g + h) − g′(f + h) − fhu + ghv}2

+ 4fg(f ′ + hu)(g′ + hv)] = 1 D4

• (ξ − fη + gζ)2 + 4fgηζ
• =

1 D4

• ξ2 − 2ξ(fη − gζ) + (fη + gζ)2

.

12

Theorem 1 1. The system is hyperbolic in the following 3 regions: (1) ηζ > 0, (2) ξ > 0, η > 0, ζ < 0, (3) ξ < 0, η < 0, ζ > 0

• 2. Elliptic regions appear in the following 2 regions:

(1) ξ < 0, η > 0, ζ < 0, (2) ξ > 0, η < 0, ζ > 0

• 3. If ξ = η = ζ = 0 at U ∗, then U ∗ is an umbilical point.

13

L1 L2 L3 elliptic region ξ < 0 ξ > 0 ζ > 0 ζ < 0 η < 0 η > 0

Figure 3: Existence of Elliptic Region (ξ > 0, η < 0, ζ > 0) 14

Integral Curves of Characteristic Vector Fields: Let A be a 2 × 2 matrix. Note that x : an eigenvectors of A ⇔

tx⊥Ax = 0,

(tx⊥x = 0). The integral curve of the characteristic vector fields: solutions to the differential equation

t ˙

U ⊥F ′(U) ˙ U = 0. The equation of the trajectory: gζdu2 + (ξ − fη + gζ)dudv − fηdv2 = 0. (3) Equivalently dv du = 1 2fη

• ξ − fη + gζ ±

√ ∆

• r

du dv = − 1 2gζ

• ξ − fη + gζ ±

√ ∆

• 15

where ∆ =

• ξ2 − 2ξ(fη − gζ) + (fη + gζ)2

. Notice that |ξ − fη + gζ| < √ ∆, if ηζ > 0, |ξ − fη + gζ| > √ ∆, if ηζ < 0. Marchesin’s model: Marchesin, Paes-Leme, Schaeffer & Shearer [17]. Theorem 2 (Existence of Umbilical Point) Assume: h(u, v) = h(1 − u − v), f(0) = g(0) = h(0) = 0, f ′′(u), g′′(v), h′′(w) > 0. Then the system of equations (1), (2) is hyperbolic and has a unique umbilical point in Ω.

16

The integral curves of the characteristic vector fields sketched:

fast integral curves slow integral curves

Figure 4: Integral Curves

N.B. 2 Impossible to construct globally in Ω fast or slow characteristic fields.

17

• 3. Quadratic Marchesin Model

Quadratic Relative Permeability Functions: f(u) = αu2, g(v) = βv2, h(u, v) = γ(1−u−v)2, α, β, γ > 0 The model equations: ∂u ∂t + ∂ ∂x

• αu2

αu2 + βv2 + γ(1 − u − v)2

• = 0,

(4) ∂v ∂t + ∂ ∂x

• βv2

αu2 + βv2 + γ(1 − u − v)2

• = 0.

(5) A unique umbilical point and the coincident characteristic speed: U ∗ = γ βγ + γα + αβ β α

• ,

λ∗ = 2αβγ(βγ + γα + αβ) αβ2γ2 + βγ2α2 + γα2β2.

18

2-phase like curves: L1 : αu − βv = 0 λ = 2γαwu D2 , R = β α

• , λ⊥ = 2αu

D , R⊥ =

• 1

−1

• ,

L2 : (β + γ)v − γ(1 − u) = 0 λ = 2αβuv D2 , R = β + γ −γ

• , λ⊥ = 2βv

D , R⊥ = 1

• ,

L3 : (α + γ)u − γ(1 − v) = 0 λ = 2αβuv D2 , R = −γ α + β

• , λ⊥ = 2αu

D , R⊥ = 1

• ,

Theorem 3 Each 2-phase like flow curve Lj, j = 1, 2, 3 is a line and the characteristic vector field R is parallel to the direction of the line.

19

Useful lemmas: Along the 2-phase like flow curve L1 : u = βτ, v = ατ, D = (α + β)(βγ + γα + αβ)τ 2 − 2γ(α + β)τ + γ. Hence Dτ = 0, if and only if (βγ + γα + αβ)τ = γ. Lemma 2 The quantity D attains its minimum at the um- bilical point. Confine our attention to the 2-phase like flow curve L1. At u = v = 0, we have λ = λ⊥ = 0 and at u + v = 1 (w = 0), λ = 0, λ⊥ > 0. Lemma 3 The characteristic speed λ attains its maximum in the interior of each 2-phase like flow curve.

20

Lemma 4 ∂λ ∂w

• U=U∗

= 2αβγ(αβ − βγ − γα) (α + β)(βγ + γα + αβ)D

• n

L1 ∂λ ∂v

• U=U∗

= 2αβ(βγ − γα − αβ) (βγ + γα + αβ)D

• n

L2 ∂λ ∂u

• U=U∗

= 2αβ(γα − αβ − βγ) (βγ + γα + αβ)D

• n

L3 Study L1 as a rarefaction or shock curve:

• 1. If αβ > βγ+γα, then L1 consists of slow rarefaction curve

and fast shock curve. The slow rarefaction curve ends in the interior of the 2-phase like flow curve.

• 2. If αβ = βγ + γα, then L1 consists of slow and fast shock

curves.

21

• 3. If αβ < βγ +γα, then L1 consists of slow shock curve and

fast rarefaction curve. The fast rarefaction curve ends in the interior of the 2-phase like flow curve. N.B. 3 Note that: if αβ ≥ βγ + γα, then βγ < γα + αβ, γα < αβ + βγ.

22

αβ > βγ + γα αβ < βγ + γα 2−phase like flow curve λ U∗ 2−phase like flow curve λ U∗ λ λ

Figure 5: Characteristic Speeds 23

U*

slow shock fast rarefaction

Case I

slow shock slow shock fast rarefaction fast rarefaction

U*

slow rarefaction fast rarefaction

Case II (αβ > βγ + γα)

slow shock slow shock fast rarefaction fast shock

Theorem 4 (Case I) Suppose that the following (all) three inequalities hold αβ < βγ + γα, βγ < γα + αβ, γα < αβ + βγ. (6) Then each 2-phase like flow curve in a neighbourhood of the umbilical point consists of slow shock curve and fast rarefaction curve. The fast rarefaction curves end in the interior of 2-phase like flow curves.

24

Theorem 5 (Case II) Suppose that one of the following three inequalities hold αβ > βγ + γα, βγ > γα + αβ, γα > αβ + βγ. (7) Then then one of the three 2-phase like flow curve in a neighbourhood of the umbilical point consists of slow rar- efaction curve and fast shock curve. Each of the other two 2-phase like flow curves has the property of Case I. The fast and slow rarefaction curves end in the interior of 2-phase like flow curves.

25

Schaeffer-Shearer’s Classification [17]: May assume that U ∗ = O and F(O) = O. Taylor expansion of the flux function F(U) near U = O: F(U) = λ∗U + Q(U) + O(1)|U|3 where λ∗ = λ1(U ∗) = λ2(U ∗) and Q : R2 → R2 is a homoge- neous quadratic mapping. Schaeffer-Shearer [17] shows that every hyperbolic quadratic mapping Q(U) with an isolated umbilical point U = O is equivalent to Q(U) = 1 2 au2 + 2buv + v2 bu2 + 2uv

• = 1

2∇C(U), (8) C(U) = 1 3au3 + bu2v + uv2. (9)

26

where a and b are two real parameters satisfying a = 1 + b2. In their classification: Case I : a < 3 4b2

• r

Case II : 3 4b2 < a < 1 + b2.

27

• 4. Undercompressive and Overcompressive Shock

Waves Rankine-Hugoniot condition: A jump discontinuity: U(x, t) = UL for x < st, UR for x > st, (10) a piecewise constant weak solution to the the conservation laws⇐ ⇒ the Rankine-Hugoniot condition: s(UR − UL) = F(UR) − F(UL). (11) The weak solution (10) satisfying (11) is often called a shock wave of speed s joining the state UL, on the left, to the state UR, on the right.

28

In the quadratic model, each 2-phase like flow curve Lj, j = 1, 2, 3 is a line and the characteristic vector field R is parallel to the direction of the line. Lemma 5 (B. Temple) Suppose that an integral curve of a characteristic vector field constitutes a line. Then the line is also a Hugoniot curve.

29

Compressive shock wave: The shock wave is said to be a j-compressive (j = 1, 2) if the speed satisfies the Lax entropy conditions: λj(UR) < s < λj(UL), λj−1(UL) < s < λj+1(UR) Here we adopt the convention λ0 = −∞ and λ3 = ∞.

1-compressive

x = st λ2

+

λ2

−

λ1

−

λ1

+

2-compressive

x = st λ2

+

λ1

−

λ2

−

λ1

+

Figure 6: Compressive Shock waves 30

Undercompressive shock wave: Undercompressive if s satisfies λ1(UR) < s < λ2(UR), λ1(UL) < s < λ2(UL)

Undercompressive x = st λ2

+

λ2

−

λ1

−

λ1

+

Figure 7: Undercompressive Shock wave 31

Overcompressive shock wave: Overcompressive if s satisfies λ1(UR) < s < λ1(UL), λ2(UR) < s < λ2(UL) Overcompressive shock waves appear only in Case II.

Overcompressive x = st λ2

+

λ2

−

λ1

−

λ1

+

Figure 8: Overcompressive Shock wave 32

Stability and Admissibility of Shock Waves: It is generally believed

• Compressive shock waves are generally stable and admis-

sibility is independent of diffusion matrices in a generic class.

• Undercompressive shock waves are stable with additional

(kinetic) condition and admissibility depends on diffusion matrices.

• Overcompressive shock waves are generally unstable.

33

Undercompressive shock waves appear in Case I and II. Theorem 6 (Case I) Suppose that the (all) three inequal- ities αβ < βγ + γα, βγ < γα + αβ, γα < αβ + βγ. hold Then on each 2-phase like flow curve, there exist un- dercompressive shock waves connecting two states on the curve.

34

Overcompressive shock waves appear only in Case II. Theorem 7 (Case II) Suppose that one of the three in- equalities αβ > βγ + γα, βγ > γα + αβ, γα > αβ + βγ.

• holds. Then on one of the three 2-phase like flow curve,

there exist overcompressive shock waves connecting two states

• n the curve.

Each of the other two 2-phase like flow curves has the property of Case I.

35

• 5. Hugoniot Loci

The Hugoniot locus of U0, denoted by H(U0): the projection

• f the set of (U, s) satisfying the Rankine-Hugoniot condition

HU0(U, s) := −s(U − U0) + F(U) − F(U0) = O (12)

• n to the U plane.

Equivalently, the set of the states U satisfying F(U)−F(U0) ∝ U − U0. For the Stone-Marchesin model: the plane cubic curve de- fined by (v−v0)(D0f(u)−Df(u0)) = (u−u0)(D0g(u)−Dg(v0)) (13) where D0 = f(u0) + g(v0) + h(w0).

36

Since the above cubic curve has a singularity at U0 = t(u0, v0), it is a rational curve. U0 is the primary bifurcation point. By introducing a parameter ξ as v − v0 = ξ(u − u0).

Lemma 6 The Hugoniot locus has a rational parametrization u − u0 = 2[αu0{h0ξ+(1+ξ)g0}−βv0ξ{h0+(1+ξ)f0}+γw0(ξf0−g0)(1+ξ)]

βξ2{h0+(1+ξ)f0}−α{h0ξ+(1+ξ)g0}+γ(ξf0−g0)(1+ξ)2

. and v − v0 = ξ(u − u0), where f0 = αu2

0, g0 = βv2 0 and h0 = γw2 0.

N.B. 4 A shock curve is a C1 piece of the Hugoniot locus.

37

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

Figure 9: α = 0.9, β = 0.8, γ = 1.0, UL = (0.20, 0.35)

Study degenerate cases. U0 on the boundary: The boundaries u = 0, v = 0 and w = 0 are themselves shock curves. Let w0 = 0 and set

38

u0 = z ≥ 0, v0 = 1 − z ≥ 0. Then u − u0 = 2αβz(1 − z) {1 − z(1 + ξ)} α(β + γ)z2ξ2 + γ {αz2 − β(1 − z)2} ξ − β(γ + α)(1 − z)2 and v − v0 = ξ(u − u0) showing that the Hugoniot locus is composed of a hyperbola plus the boundary w = 0. There are no secondary bifurcation points.

39

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

Figure 10: Hugoniot Locus of U0 : w0 = 0. (left: 0 < z < 1

2, right: 1 2 < z < 1)

U0 on the 2-phase like curves: The 2-phase like curves Lj = 0 (j = 1, 2, 3)are themselves shock curves. Let U0 ∈ L1 and set u0 = βz α + β, v0 = αz α + β (0 ≤ z ≤ 1).

40

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

Figure 11: Hugoniot Locus of U0 ∈ L1 and Secondary Bifurcation Point UII

Then u − u0 = 2αβz2(1 + ξ)

• γ(1 − z) − αβ

α+βz

• αβγz2(1 + ξ)2 + αβz2(1 + ξ)(βξ + α) + γ(α + β)(1 − z)2ξ

41

and v − v0 = ξ(u − u0) showing that the Hugoniot locus is a part of a hyperbola plus the 2-phase like curve L1. By setting ξ = α β, we have uII − u0 = 2βz2 γ(1 − z) − αβ

α+βz

• γ(α + β)z2 + 2αβz2 + γ(α + β)(1 − z)2.

Here UII = t(uII, vII) represents the secondary bifurcation point. Notice that the umbilical point is expressed as u∗ − u0 = β

• γ(1 − z) − αβ

α+βz

• βγ + γα + αβ

. By direct computation, we find UII = U ∗ ⇔ z = 1 2

42

and uII < u∗ for z = 0, 1. Thus, by setting U ∗

0 = U|z=1

2 =

t β 2(α+β), α 2(α+β)

• , we have

u∗

0 − u∗ =

β {αβ − (βγ + γα)} 2(α + β)(βγ + γα + αβ). Theorem 8 If αβ < βγ + γ then u∗

0 < u∗

and uII < u∗ if U0 / ∈ U ∗

0U ∗,

uII ≥ u∗ if U0 ∈ U ∗

0U ∗.

43

N.B. 5 The three 2-phase like curves are Hugoniot loci of the umbilical point. However these curves are not parallel to asymptotes of generic Hugoniot curves.

• 3
• 2
• 1

1 2 3

• 3
• 2
• 1

1 2 3

Figure 12: α = 0.9, β = 0.8, γ = 1.0, UL = (0.20, 0.35) 44

• 6. Invariant Region

The boundary of Ω consists of 3 segments each of which is a rarefaction curve and at the same time a shock curve. Hoff [8] says that Ω is an invariant region. Lax-Friedrichs and Godunov approximations: Lax- Friedrichs and Godunov approximations are based on taking averages and integral averages,respectively, of solution to Rie- mann problems. Since Ω is convex, we have Theorem 9 Suppose that all initial all states are suffi- ciently close and Courant-Friedrichs-Lewy condition holds. Then Ω is an invariant region for Lax-Friedrichs and Go- dunov approximations.

45

• 7. Entropy Functions

Entropy Equation: The entropy-entropy flux pair (H, Q) satisfies the compatibility condition ∇Q = F ′(U)t∇H. The integrability condition for Q yields the second order linear partial differential equation fηHuu + (ξ − fη + gζ)Huv − gζHvv = 0. (14) Proposition 1 (Lax [11]) The linear partial differential equation (14) is strictly hyperbolic if and only if (u, v) be- longs to the strictly hyperbolic region of the original equa- tions (1), (2)

46

The characteristic condition is fηφ2

u + (ξ − fη + gζ)φuφv − gζφ2 v = 0

(15) and solutions to the differential equation fη dv du 2 − (ξ − fη + gζ) dv du

• − gζ = 0

are called bicharacteristic curves. Proposition 2 (Lax [11]) The bicharacteristic curves are the integral curves of the characteristic vector fields and the Riemann invariants φ = w, z satisfy the characteristic condition (15). Let λ, µ (λ < µ) denote the characteristic speeds in a hy- perbolic region. By using the Riemann (characteristic) coor-

47

dinates w, z, the entropy equation (14) has the form Hwz + 1 λ − µ {λzHw − µwHz} = 0. (16) Suppose that the rectangle [w0, w1] × [z0, z1] is contained in a strictly hyperbolic region. Entropy functions are constructed by assigning the boundary condition H(w, z0) = Φ(w), H(w0, z) = Ψ(z), (17) which is the Goursat problem. Theorem 10 Let Φ(w), Ψ(z) are Lipschitz continuous func- tion defined on the interval [w0, w1], [z0, z1], respectively. Then there exists a unique entropy function H(w, z) de- fined on the rectangle [w0, w1]×[z0, z1] satisfying the bound- ary condition (17). Full proof is found in Sobolev [18].

48

Covering Regions: It is impossible to construct globally in Ω the fast and slow characteristic fields. Let Ωǫ denote Ω with removing small neighbourhoods of three vertices and the umbilical point. Theorem 11 There exist a double covering of Ωǫ and a triple covering of corresponding region in wz-plane such that the fast and slow characteristic fields are globally well

• defined. The entropy equation (14) and its characteristic

form (16) are well-deifined on those covering region.

49

w = −1 A D E F G H I J L B C K w = −1 w = −1 A B C F H J K L D E G I w = 1 w = 1 w = 1 w = 1 w = 1 w = 1 w = −1 w = −1 w = −1 : w = const. z = 1 A D E F G H I J L B C K z = −1 z = 1 A B C F H J K L D E G I z = −1 z = −1 z = 1 z = 1 z = −1 z = 1 z = −1 z = 1 z = −1 z = 0 z = 0 z = 0 z = 0 z = 0 z = 0 : z = const.

Figure 13: Fast and Slow Integral Curves 50

z = 1 z = −1 w = 1 w = −1 G E F H J G H I z = 1 z = −1 w = 1 w = −1 K I J L A K L B z = 1 z = −1 w = 1 w = −1 C B A D F C D E

Figure 14: Covering of wz-Plane 51

References

[1] Asakura F. & Yamazaki M. (2005) Geometry of Hugoniot curves in 2×2 systems of hyperbolic conservation laws with quadratic flux functions, IMA J. Appl. Math., 70, no. 6, 700–722. [2] Asakura F. & Yamazaki M. (2009) Viscous Shock Profiles for 2 × 2 Systems of hyperbolic conservation laws with an umbilical point, J. Hyperbolic Differential Equations, Vol. 6, 483–524. [3] Bell J.B., Trangenstein J. A. & Shubin G. R. (1986) Con- servation laws of mixed type describing three-phase flow in porous media, SIAM J. Appl. Math. 46, No. 6, 1000–1017 [4] Courant R., Friedrichs K. O. & Lewy H. (1928) ¨ Uber

52

die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann. , 100, 32–74 (English Translation: IBM Journal, 11, 215-234, 1967). [5] Friedrichs K. O. (1954) Symmetric hyperbolic linear dif- ferential equations, Comm. Pure. Appl. Math., 7, 345–392 [6] Gomes M. E. S. (1989) Riemann problems requiring a vis- cous profile entropy condition, Adv. Appl. Math., 10, 285– 323. [7] Harten A., Lax P. D. & Van Leer. (1983) On upstream dif- ferencing and Godunov-type schemes for hyperbolic con- servation laws, SIAM Review, 25, 35–61. [8] Hoff D. (1985) Invariant regions for systems of conserva- tion laws, Trans. Amer. Math. Soc., 289, 591–610.

53

[9] Holden L., Risebro H. & Tveito A. (1995) Maximum prin- ciple for a class of conservation laws, SIAM J. Appl. Math., 55, 651–661. [10] H¨

• rmander L. (1997) Lectures on Nonlinear Hyperbolic

Differential Equations, Math´ ematiques & Applications 26, Springer. [11] Lax P. D. (1971) Shock waves and entropy. E. Zaran- tonello (ed.), Contributions to nonlinear Functional Analysis, Academic Press, New York, 603–634. [12] John F. (1952) On integration of parabolic equations by difference methods, Comm. Pure. Appl. Math., 5, 155–211 [13] John F. (1978) Partial Differential Equations (3rd edi- tion), Applied Mathematical Sciences 1, Springer

54

[14] Majda A. & Pego R. (1985) Stable Viscosity Matrices fo Systems of Conservation Laws, J. Differential Equations, 56, 229–262. [15] Marchesin D. & Plohr B. (2001) Theory of Three-Phase Flow Applied to Water-Alternating-Gas Enhanced Oil Re- covery, Hyperbolic Problems; Theory, Numerics, Appli- cations, Vol.II, Birkh¨ auser Verlag, 693–702. [16] Medeiros H. B. (1992) Stable Hyperbolic Singularities for Three-Phase Flow Models in Oil Reservoir Simulation, Acta Applicandae Mathematicae, 28, 135–159. [17] Schaeffer D. & Shearer M. (1987) The classification of 2×2 systems of non-strictly hyperbolic conservation laws, with applications to oil recovery, Comm. Pure Appl. Math., 40, 141–178.

55

[18] Sobolev S. L. (1964) Partial Differential Equations of Matheamtical Physics, Pergamon Press (1989, Dover edi- tion) [19] Stone H. (1970) Probability model for estimating three- phase relative permeability, J. Petr. Tech., 22, 214–218.

56