Maximal oil recovery by simultaneous condensation of alkane and steam

Maximal oil recovery by simultaneous condensation of alkane and steam

J. Bruining*

Dietz Laboratory, Centre of Technical Geoscience, Mijnbouwstraat 120, 2628 RX Delft, The Netherlands

D. Marchesin†

Instituto Nacional de Matemática Pura e Aplicada, Estrada Dona Castorina 110, 22460-320 Rio de Janeiro, Rio de Janeiro, Brazil

共Received 4 November 2005; revised manuscript received 12 October 2006; published 26 March 2007兲

This paper deals with the application of steam to enhance the recovery from petroleum reservoirs. We formulate a mathematical and numerical model that simulates coinjection of volatile oil with steam into a porous rock in a one-dimensional setting. We utilize the mathematical theory of conservation laws to validate the numerical simulations. This combined numerical and analytical approach reveals the detailed mechanism for thermal displacement of oil mixtures discovered in laboratory experiments. We study the structure of the solution, determined by the speeds and amplitudes of the several nonlinear waves involved. Thus we show that the oil recovery depends critically on whether the boiling-point of the volatile oil is around the water boiling temperature, or much below or above it. These boiling-point ranges correspond to three types of wave struc-tures. When the boiling point of the volatile oil is near the boiling point of water, the striking result is that the speed of the evaporation front is equal or somewhat larger than the speed of the steam condensation front. Thus the volatile oil condenses at the location where the steam condenses too, yielding virtually complete oil recovery. Conversely, if the boiling point is too high or too low, there is incomplete recovery. The condensed volatile oil stays at the steam condensation location because the steam condensation front is a physical shock. DOI:10.1103/PhysRevE.75.036312 PACS number共s兲: 44.30.⫹v, 44.35.⫹c

I. INTRODUCTION

Steam drive is an economical way of producing oil and is used worldwide for heavy oil. An overview of the last forty years of steam drive recovery in California is given in Ref. 关1兴. Steam drive is also considered an efficient method to clean polluted sites关2–4兴. During the steam drive, however, a certain amount of oil is left behind in the steam swept zone 关5兴.

In the late 1970s Dietz关5兴 proposed to add small amounts of volatile oil to the steam to reduce the oil left behind. Similar ideas were put forward independently by Farouq-Ali 关6兴. The volatile oil coinjected with the steam in almost in-finitesimal amounts would ideally condense at the same lo-cation where the steam condenses. The condensed volatile oil acts as a solvent for the heavy oil. As such it pushes the oil away from the steam-swept zone leaving no oil behind 共see Fig.1兲. At the time the crucial importance of the boiling temperature of the volatile oil was not suspected. Experi-ments investigating the mechanism are described in Refs. 关5–9兴. Still, there was a discrepancy between the original idea and the experimental observations. At least 5 wt % 共volatile oil/water兲 was required to reduce considerably the saturation of the oil left behind关5兴. However, it is possible that the requirement of this large percentage was caused by transient effects in the experiments. One of the goals of our work is to clarify this point.

In his pioneering experiment, Willman, in 1961, used a large percentage of initially present volatile oil关9兴. His ex-periment led to the belief that any volatile oil component,

initially present in the oil, would lead to virtually complete recovery from the steam-swept zone. Therefore, the virtue of adding volatile oil was criticized at the time. The second goal of our work is to establish the difference between steam-drive recovery with coinjection of volatile oil and recovery of oil already containing a fraction of volatile oil. It can be expected that an efficient condensed volatile oil region is too short for the resolution of standard simulators.

Our approach关10–13兴 is to simplify the model equations in such a way that the essential elements are retained 关14兴 while avoiding the complexities of solving pressure equa-tions and nonlinear compositional equaequa-tions at every grid cell. As such the model is a straightforward extension of a one-dimensional共1D兲 model proposed by Ref. 关15兴, but al-lowing for immiscible three-phase flow in the steam zone 关16,17兴 共see also Ref. 关18兴兲. The simplification is accom-plished by the assumption that the steam drive runs at con-stant pressure as to the thermodynamical behavior; any pres-sure increase causes an immediate production of fluids. Therefore the pressure equation decouples and we can solve the transport equations locally, reaching resolutions that are unattainable in standard simulators.

The solution of these simplified transport equations is ob-tained by following each physical state in space time, using the method of characteristics. If the transport equation were linear with constant coefficients all states would move at the same characteristic speed and the wave profiles would re-main unchanged. In our case, however, the equations are nonlinear, therefore characteristic speeds depend on the state. If characteristic speeds increase in the flow direction, states spread out giving rise to a rarefaction fan共rarefaction wave兲. On the other hand, if characteristic speeds decrease the states collapse on each other giving rise to a discontinuity or a shock wave. It is this nonlinear collapse that both generates and stabilizes shock waves. The mathematical theory of *Electronic address: J.bruining@citg.tudelft.nl

†

these nonlinear waves is very well developed关19,20兴. How-ever, in our problem mass transfer can occur, viz., evapora-tion or condensaevapora-tion, giving rise to evaporaevapora-tion or condensa-tion rarefaccondensa-tions or shocks, which are not so well known. The relative speeds of the waves occurring in the solution are crucial in determining the physical phenomena and the dis-tinctive behavior in terms of the values of the physical prop-erties.

Knowing the solution obtained by the method of charac-teristics has three advantages. Firstly, it is a time-asymptotic solution, which is relevant at the field scale. Secondly, it allows us to validate the numerical solution. Thirdly, it al-lows the study of bifurcation phenomena, i.e., change of structure of solutions under different injection conditions. The bifurcations of this model in the absence of thermal effects are described in Refs.关21–24兴. 共See also the review in the appendix of Ref.关25兴兲.

The model we used carries three important simplifica-tions. Firstly, the diffusion mixing between volatile oil and heavy oil in the liquid phase and between volatile oil vapor and water vapor in the gaseous phase are disregarded. The model is not valid for extremely low injection rates, where capillary diffusion dominates convection, because we ignore capillary effects. Finally, we do not specify a detailed model for the kinetics of the condensation process关13兴. These as-pects determine the internal structure of the shocks, which sometimes affect the structure of the whole Riemann solu-tion, and are subjects for future work关23兴.

The range of validity of these simplifications can be ex-pressed in terms of dimensionless numbers关26兴. For diffu-sion effects to be negligible, both Péclet numbers 关Pe = Luinj_/共␸D

i兲兴, i.e., the one based on molecular diffusion

共D1兲 and the one based on capillary diffusion 共D2兲, must be much larger than one. For field conditions L is the distance between wells. The Péclet number is at least one million, even for capillary-diffusion phenomena. For laboratory con-ditions it is a factor 100 smaller, but still PeⰇ1. Thirdly, the ratio between the rate of mass transfer between phases and convective mass transport, expressed by the Damkohler number关26兴, must be very big so that the thickness of the condensation zone can be disregarded. This aspect is dis-cussed in Ref.关27兴, where it is shown that a practical value of the Damkohler number Da= qbL / uinjwould be of the

or-der of 108_{. Here q}

b关s−1兴 is the rate of steam condensation. In

the same paper it is shown that local equilibrium is obtained when Da⬃104_{. So the condition of local thermodynamic}

equilibrium is definitely satisfied at Da= 108. This shows that to leading order our model, where we use local thermody-namic equilibrium and disregard diffusion effects, is correct. However, for more precise and quantitative statements these effects must be analyzed. This is however, beyond the scope of this paper.

Finally we also use Darcy’s law without inertia correc-tion, which requires that the Reynolds number␳vdp/␮based

on the grain size, is not larger than one. A typical value for field conditions is Re= 0.007. Therefore for field situations these conditions are always satisfied and they were satisfied in most of our laboratory experiments.

Section II describes the physical model and the relevant thermodynamical relations. The flow is described by balance equations in Sec. III. Self-similar waves, i.e., rarefaction and shocks, are analyzed in Sec. IV. An implicit finite difference method requiring the solution of small matrices is described in Sec. V. Section VI summarizes earlier results on the injec-tion of steam displacing heavy oil. Our results concerning the solution structure and the recovery in terms of the boiling point of volatile oil are described and discussed in Sec. VII. We summarize our conclusions in Sec. VIII. Appendix A contains further details. Appendix B describes physical quan-tities, symbols, and values. Some calculations are found in Appendixes C and D.

II. PHYSICAL MODEL A. Flow of fluids

The model is based on conventional models for steam drive关28,29兴. We consider the injection of steam and volatile oil into a linear horizontal porous rock cylinder with constant porosity and absolute permeability共see Fig. 1兲. The tube is completely thermally isolated. The injection temperature is determined by the three-phase equilibrium condition for the given volatile-oil and steam injection ratio. The cylinder is originally filled with oil and water. The oil consists of dead oil, i.e., oil with a negligible vapor pressure, possibly with dissolved volatile oil. The dissolution of volatile or dead oil in water is negligible. Three-phase flow occurs in the high temperature zone, while oil and water flow occurs in the low temperature zone. The fluids are in local thermodynamic equilibrium. Physical quantities are evaluated at a represen-tative constant pressure throughout the cylinder; this is a good approximation if the total pressure variation is small relative to the total pressure. It is certainly valid in laboratory

experiments. Thermal expansion of liquids is disregarded. The liquid volatile oil and dead oil heat capacities are not the experimental heat capacities, but slightly adapted so that the enthalpy of the oleic phase is independent of composition. This minor adjustment leads to a major simplification of the mathematical analysis. All fluids are considered incompress-ible. We assume Darcy’s law for multiphase flow 关30,31兴. The cylinder diameter is sufficiently small so that capillary forces equalize the saturation in the transverse radial direc-tion and temperature is homogeneous radially. As the flow is horizontal we ignore gravity effects.

B. Thermodynamic fundamentals

Our interest is confined to共1兲 three-phase flow, i.e., flow of the aqueous共w兲, oleic 共o兲, and gaseous 共g兲 phases in the steam zone and共2兲 two-phase flow, i.e., flow of the aqueous and oleic phases in the liquid zone. For liquids, we distin-guish between an aqueous 共waterlike兲 phase and an oleic 共oil-like兲 phase because they do not mix. We use the follow-ing convention: the first subscript 共w,o,g兲 refers to the phase, the second subscript共w,v,d兲 refers to the component, i.e., water, volatile oil, and dead oil. Capital subscripts 共W,V,D兲 are used to denote phases consisting of a single component. The densities of the pure liquids are denoted as

␳W,␳V, and␳D. The densities of the pure vapors, i.e., water

and volatile oil are denoted by␳gW,␳gV.

We disregard any heat or volume contraction effects re-sulting from mixing. The concentration关kg/m3_{兴 of volatile} 共dead兲 oil in the oleic phase is denoted as ␳ov 共␳od兲. The

concentration of water vapor 共volatile oil兲 in the gaseous phase is␳gw 共␳gv兲. For ideal fluids we obtain

␳ov ␳V +␳od ␳D = 1, ␳gw ␳gW + ␳gv ␳gV = 1. 共1兲

The densities of the pure liquids␳V,␳D关kg/m3兴 are

consid-ered to be independent of temperature, and the densities of the pure vapors to obey the ideal gas law, i.e.,

␳gW=

MWP

RT , ␳gV= MVP

RT , 共2兲

where MW, MV denote the molar weights of water and

vola-tile oil, respectively. T is the temperature and the gas con-stant is R = 8.31关J/mol K兴. P is not a variable in this prob-lem, but the fixed prevailing pressure value; here we use one atmosphere, because most of the experiments were carried out at atmospheric pressure.

The water vapor pressure Pw is determined by the

Clausius-Clapeyron equation关32兴 Pw共T兲 = Poexp

冋

− MW R ⌳W共Tb w_兲

冉

1 T− 1 Tb w

冊

册

, 共3兲 where⌳W 共Tb

w_{兲 关J/kg兴 is its evaporation heat at its normal}

boiling temperature Tb

w_{共K兲 at P}

o, the atmospheric pressure.

We also use Clausius-Clapeyron for the volatile-oil vapor pressure. In addition, we use Raoult’s law关32兴, which states that the vapor pressure of volatile oil is equal to its pure

vapor pressure times the mole fraction xov of volatile oil in

the oil phase. Therefore we obtain

Pv共T兲 = xovPoexp

冋

− MV R ⌳V共Tb v_兲

冉

1 T− 1 Tb v

冊

册

, 共4兲

where⌳V共Tbv兲 is the evaporation heat of the volatile oil at its

normal boiling temperature T_bv. We assume that the prevail-ing pressure P is the sum of the two vapor pressures. From Eqs.共3兲 and 共4兲, and P= Pw共T兲+ Pv共T兲, we find for the mole

fraction of volatile oil in the liquid-oil phase xov共T兲

=共␳ov/ MV兲/共␳ov/ MV+␳od/ MD兲, xov共T兲 = P − Poexp

冋

− MW R ⌳W共Tb w_兲

冉

1 T− 1 Tb w

冊

册

Poexp

冋

− MV R ⌳V共Tb v_兲

冉

1 T− 1 Tbv

冊

册

. 共5兲

From this we derive an expression for the volatile-oil con-centration in the oleic phase,

␳ov=

xov␳D␳VMV

xov␳DMV+共1 − xov兲␳VMD

. 共6兲

Note that in the gaseous phase no dead oil component is present, whereas in the oleic phase volatile and dead oil are present. Figure2shows the projections of the phase diagram of cyclobutane共left兲, and heptane 共right兲 on the plane of the temperature vs the volatile-oil mole fraction. The special three-phase point “3ph” indicates where pure liquid volatile oil, liquid water, and vapor coexist. Other three-phase points are on⌫, the curve where liquid water, volatile-oil and dead-oil mixtures, and vapor coexist as explained below. For each

T, the mole fraction of volatile oil in the vapor phase on⌫ is

indicated by ygv=关P− Pw共T兲兴/ P. This equation is used for

heptane in Fig.2共right兲 to find the lower branch ⌫ extending from 共ygv, T兲=关0.0,373.15 共K兲兴→共0.551,351.71兲. At the

latter point x_ov= 1. Similarly, the lower branch ⌫ extends from 共y_gv, T兲=关0.0,373.15 共K兲兴→共0.983,285.20兲 for cy-clobutane 关Fig. 2 共left兲兴. For dodecane it extends from 共ygv, T兲=关0.0,373.15 共K兲兴→共0.0356,371.98兲. All its

fea-tures of interest occur near the very small branch at the left, which makes a figure less illustrative. Therefore we do not show it. In Fig.2we assume that the prevailing pressure is atmospheric, i.e., P = Po.

Furthermore, Fig. 2 contains projections of 3D figures with the temperature as the vertical axis, the volatile-oil frac-tion in the vapor phase, ygv as the horizontal axis, and the

composition of the oil phase x_ovas the axis perpendicular to the paper. The projection is made on a surface for which

xov= constant.

Figure 2 共left兲 contains four phase diagrams for xov

= 1 , 0.6, 0.2, 0.1. Consider as an example the phase diagram for x_ov= 0.2, i.e., the behavior of liquid oil with a volatile-oil mole fraction of xov= 0.2. This phase diagram consists of the

curve兵xov= 0.2其 共see the next paragraph兲 and the part of the

intersection point of these two curves, liquid water and liquid-oil phase with composition xov= 0.2 are in equilibrium

with vapor with a volatile-oil fraction y_gv indicated in the horizontal axis. On ⌫, left of the intersection point liquid water is in equilibrium with vapor with a volatile-oil fraction

ygvread from the horizontal axis. Above these curves there is

only vapor. Because in all cases considered here there is liquid water, we only use the left curve, the curve⌫, for the three-phase zone and the region underneath for the two-phase zone.

The procedure to find the branches on the right emanating from⌫ for which x_ov= 0.2, 0.4, 0.6, 1.0 is the following. We pick a value for x_ov, choosing a curve among these branches. This curve is described by the graph of ygv共T兲= Pv共T兲/ P,

where Pv共T兲 is obtained from Eq. 共4兲. Recall that we used

ygv共T兲=关P− Pw共T兲兴/ P to obtain the plot of ⌫. Only at the

intersection point of the curves liquid oil with the chosen xov,

liquid water and vapor are in equilibrium and we have that

y_gv共T兲=关P− Pw共T兲兴/ P= Pv共T兲/ P. On 兵xov= 1其 there is no

dead oil, rather there is volatile-oil vapor besides water vapor and liquid volatile oil. For xov⬍1 all branches to the right of

⌫ contain dead oil too. Therefore, these branches with con-stant xov describe the two-phase oleic-gaseous equilibrium;

there is no liquid water. We can use Eq.共5兲 to obtain expres-sions for the concentrations ␳ov共T兲, ␳od共T兲, ␳gw共T兲, and

␳gv共T兲.

III. BALANCE EQUATIONS

The energy conservation equation in terms of enthalpy is given as关33兴

⳵

⳵t关Hr+␸共HWSw+ HoSo+ HgSg兲兴

+ ⳵

⳵xu共HWfw+ Hofo+ Hgfg兲 = 0, 共7兲

where the enthalpies per unit volume HW, Ho, and Hg are

defined in the table in terms of densities and enthalpies per unit mass hW, hoV, hoD, hgV, hgW. These enthalpies depend on

temperature 共and on the fixed pressure兲. The enthalpy of

steam in the gaseous phase is hgW, and hWis the enthalpy of

water in the aqueous phase, while hgV is the enthalpy of

volatile oil in the gaseous phase. Furthermore, hoV and hoD

are the enthalpies of liquid volatile oil and dead oil. Their values are chosen so that the heat capacity per unit volume

Hois only a function of temperature 共see Appendix B兲. The

rock enthalpy Hr is per unit volume. The saturation of the

oleic, aqueous, and gaseous phases are So, Sw, Sg, while fo,

fw, fg are their fluxes, defined in Eq. 共B16兲. We use u to

denote the total Darcy flow velocity and␸the constant rock porosity. We can write for the mass conservation equations of water, volatile oil, and total oil关34兴,

␸_⳵⳵ t共␳gwSg+␳WSw兲 + ⳵ ⳵xu共␳gwfg+␳Wfw兲 = 0, ␸⳵ ⳵t共␳gvSg+␳ovSo兲 + ⳵ ⳵xu共␳gvfg+␳ovfo兲 = 0, ␸_⳵⳵ t

冉

␳gv ␳V Sg+ So

冊

+ ⳵ ⳵xu

冉

␳gv ␳V fg+ fo

冊

= 0. 共8兲

Equations共8兲 and 共7兲 can be written in condensed form as

⳵ ⳵tGᐉ+

⳵

⳵xuFᐉ= 0 forᐉ = w,v,o,T. 共9兲

We use the subscriptᐉ to denote the components 共w,v,o兲 and the energy 共T兲. Notice that in the three-phase zone, G_ᐉ and F_ᐉ are functions of the variables Sw, Sg, T and the

de-pendent variables of Eq. 共9兲 are Sw, Sg, T, and u. In the

two-phase zone Eqs.共8兲 and 共7兲 simplify by using␳gw=␳gv

= 0. Here fwdepends on Swand T, foon Sw,␳ov, and T. The

dependent variables in the two-phase zone in Eq.共9兲 are Sw,

␳ov, T, and u. Thus a state in the three-phase zone is defined

by the values of Sw, Sg, T, u and in the two-phase zone by the

values of Sw,␳ov, T, u.

IV. ANALYSIS OF ELEMENTARY WAVES

Considering Eq. 共9兲 and the fact that we use constant injection conditions and homogeneous initial data we

serve that solutions must exist that are invariant with respect to scaling x→ax, t→at, where a is any positive constant. Such solutions depend only on the similarity coordinate x / t and are called Riemann solutions. These solutions represent large-time asymptotic solutions for many initial and bound-ary data. Standard theory of conservation laws say that Rie-mann solutions consist of sequences of smooth rarefaction waves, discontinuities or shocks, and constant states. Shock waves satisfy the Rankine-Hugoniot共RH兲 conditions, which express mass conservation. We refer the interested reader to Smoller 关20兴 and Dafermos 关19兴. Excellent engineering in-troductions in this field can be found in the papers by Pope 关35兴, Hirasaki 关36兴, and Dumoré, Hagoort, and Risseeuw 关37兴, and in the book by Lake 关34兴.

The theory of nonlinear conservation laws relates the speed of a shock with its left and right states through the RH conditions. We find explicit formulas for RH conditions for all shocks, including condensation shocks, one of which is the steam condensation front 共SCF兲. We derive the charac-teristic speeds for rarefaction waves. We have also obtained the rarefaction curves, which represent the rarefaction waves, but we omit their lengthy derivations here. We have used these formulas to verify the correctness of every single wave found numerically in Sec. VII. The concatenation of the waves according to speed and the extended Lax entropy con-ditions关20,19兴 were verified as well. As far as the authors are concerned the treatment of the velocity variable is original as there are no time derivatives for this variable in the equa-tions.

A. Shocks

We use the standard notation for the jump of a quantity U across a shock as关U兴=U+_{− U}−_{. The RH conditions for Eqs.} 共8兲 and 共7兲 can be written as follows 关for a shock with speed v and left and right states共⫺兲 and 共⫹兲兴:

−␸v关␳gwSg+␳WSw兴 + 关u共␳gwfg+␳Wfw兲兴 = 0, 共10兲 −␸v关␳_gvSg+␳ovSo兴 + 关u共␳gvfg+␳ovfo兲兴 = 0, 共11兲 −␸v

冋

␳gv ␳V Sg+ So

册

+

冋

u

冉

␳gv ␳V fg+ fo

冊

册

= 0, 共12兲 −v关Hr+␸共HWSw+ HoSo+ HgSg兲兴 +关u共HWfw+ Hofo+ Hgfg兲兴 = 0. 共13兲

Whenv, u+_{, u}−_{solve the equations above, then av, au}+_{, au}− also solve the equations.

We distinguish six kinds of shocks. 共1兲 The volatile-oil evaporation shock共speed vE兲, with three-phase conditions at

the left. Its main feature is that the volatile-oil concentration increases in the downstream共right兲 direction. The tempera-ture, saturations, and velocity change across the shock.共2兲 The steam condensation shock 共speed vSCF兲, with a

three-phase condition at the left. The vapor saturation decreases drastically in the downstream direction. Again all the quan-tities change across the shock.共3兲 The volatile-oil condensa-tion shock 共speed vC兲, with a three-phase condition at the

left. The volatile-oil concentration decreases in the down-stream direction.共4兲 The volatile-oil two-phase composition shock共speed v_␷兲, which is a contact discontinuity. A contact discontinuity represents the moving interface between two fluids in the same phase. In reality such an interface is not infinitesimally thin. In loose mathematical terms, a contact discontinuity is defined as a shock for which the character-istic speeds at the right and left are equal to the shock speed. 共5兲 The saturation shock 共speed vS兲. Only the saturations

change, while temperature, composition, and the velocity are constant, so that Eq. 共13兲 does not play a role. 共6兲 The Buckley-Leverett shock共speed vBL兲, with only the liquid oil

and water phases present. All quantities except the liquid saturations are constant, so that Eq.共13兲 again does not play a role.

B. Characteristic speeds

Using G_ᐉand F_ᐉfrom Eq. 共9兲, we define G_ᐉn=⳵Gᐉ ⳵Vn , F_ᐉn=⳵Fᐉ ⳵Vn , F_ᐉu= ⳵ ⳵u共uFᐉ兲 = Fᐉ, 共14兲

where V =共V1, V2, V3兲=共Sw, Sg, T兲. Note that G_ᐉdoes not

de-pend on u. Without loss of generality u⬎0 and Eq. 共9兲 can be rewritten forᐉ=w,v,o,T as

兺

n=w,g,T

冉

G_ᐉn⳵Vn ⳵t + uFᐉn ⳵Vn ⳵x

冊

+ Fᐉu ⳵ln u ⳵x = 0. 共15兲

Let us consider solutions of Eq. 共9兲 that depend on 共x,t兲 through the similarity coordinate␩= x / t. Then Eq.共15兲 with 共Sw, Sg, T , ln u兲† denoting a column vector becomes

M

冉

Sw,Sg,T, ␩ u

冊

d d␩共Sw,Sg,T,ln u兲 †_{= 0, 共16兲} where M共V,␭兲 with ␭=␩/ u is the 4⫻4 matrix

冢

F11−␭G11 F12−␭G12 F13−␭G13 F1

F21−␭G21 F22−␭G22 F23−␭G23 F2

F₃₁−␭G31 F₃₂−␭G32 F₃₃−␭G33 F₃ F41−␭G41 F42−␭G42 F43−␭G43 F4

冣

. 共17兲

We have replaced the subscripts w ,v , o , T by 1,2,3,4.

Ex-plicit expressions for the flux functions, the accumulation function, and their derivatives can be found in Appendix D.

1. Three-phase flow

冨

⳵fw ⳵Sw −␭␸ ⳵fw ⳵Sg A13共␭兲 A1 ⳵fg ⳵Sw ⳵fg ⳵Sg −␭␸ A33共␭兲 A₃ 0 0 A₂₃0 −␭A231 A₂ 0 0 F43−␭G43 F4

冨

=

冨

⳵fw ⳵Sw −␭␸ ⳵fw ⳵Sg ⳵fg ⳵Sw ⳵fg ⳵Sg −␭␸

冨

冏

A₂₃0 −␭A₂₃1 A₂ F43−␭G43 F4

冏

, 共18兲

where A1, A3, A2, A230, A123,F4, F43, and G43 depend on V, while A13共␭兲, A33共␭兲 are linear expressions in ␭ with coeffi-cients that depend on V. The calculation and expressions are found in Appendix C. Let us define D =

共

⳵fw ⳵Sw− ⳵fg ⳵Sg

兲

2 + 4⳵fw ⳵Sg ⳵fg ⳵Sw and T = ⳵fw ⳵Sw+ ⳵fg ⳵Sg. The

slow and fast characteristic speeds for saturation rarefaction waves are given as

␭S,1=

− T −

冑

D

2␸ , ␭S,2=

− T +

冑

D

2␸ . 共19兲

The corresponding characteristic vectors have constant T , u; only the saturations vary along these waves. Within the satu-ration triangle, spanned by Sw, So, Sg, which add to one, there

is a point where X₁= 0. Here the two characteristic speeds coincide, giving rise to a rich wave structure共see, e.g., Ref. 关24兴兲.

The third and last characteristic speed is

␭c=␭c共Sw,Sg,T兲 = 1 ␸ A2F43− A23 0_F 4 A₂G₄₃+ A₂₃1F₄. 共20兲 This is associated with a condensation rarefaction wave, in which all quantities vary.

2. Two-phase flow

In the absence of the gaseous phase, there are three kinds of rarefaction waves. One is the thermal rarefaction wave, along which Sw,␳ov, T change. Its speed is

␭T=

Ho

⬘

+共Hw

⬘

− Ho

⬘

兲fw

Hr

⬘

+␸关Ho

⬘

+共Hw

⬘

− Ho

⬘

兲Sw兴

. 共21兲

Then we have the Buckley-Leverett rarefaction with speed␭BL, along which the liquid saturation Swchanges.

Fi-nally, we have the composition wave with speed␭C, which is

a contact discontinuity, along which the composition and the liquid-water saturation change. The speeds ␭BL共Sw,␳ov, T兲

and␭C共Sw,␳ov, T兲 are ␭BL= 1 ␸ ⳵fw ⳵Sw , ␭C= 1 ␸ 1 − fw 1 − Sw . 共22兲

We have computed the characteristic vectors for both two-and three-phase flow, but we do not provide the formulas here. They are necessary to compute the rarefaction waves.

V. NUMERICAL SOLUTION OF THE EQUATIONS We will use the notation V =共Sw, Sg, T兲 in the three-phase

region and V =共Sw,vov, T兲 in the two-phase liquid region.

A. Upstream scheme

Consider Eq.共9兲, where the fluxes F_ᐉare functions of V. We can write the upstream implicit finite-difference scheme 共ᐉ=w,v,o,T兲

G_ᐉm共t + ⌬t兲 + um共t + ⌬t兲F_ᐉm共t + ⌬t兲

= G_ᐉm共t兲 + 共⌬t/⌬x兲um−1共t + ⌬t兲F_ᐉm−1共t + ⌬t兲, 共23兲 where m denotes the grid-cell number. The unknowns are

um_{共t+⌬t兲 and the three components of V}m_{共t+⌬t兲, which}

show up in the expressions for G_ᐉm共t+⌬t兲 and F_ᐉm共t+⌬t兲. Let us rewrite Eq.共23兲 and shorten the unknowns as fol-lows: um_{共t+⌬t兲 as u and V}m_{共t+⌬t兲 as V. We obtain the}

non-linear implicit scheme

G_ᐉ共V兲 + 共⌬t/⌬x兲uF_ᐉ共V兲 = R_ᐉm,m−1, 共24兲 where we have introduced the notation R_ᐉm,m−1for the right-hand side of Eq. 共23兲. We assume that F_ᐉm−1共t+⌬t兲 and um−1_{共t+⌬t兲 have been precomputed when solving the}

previ-ous cell m − 1, which may be in a phase condition different from that of cell m. We emphasize that R_ᐉm,m−1does not de-pend on the condition of cell m at the new time t +⌬t.

B. Solution of the nonlinear system

The system共24兲 is solved using Newton-Raphson. Given an approximate solution in the kth iterationVk_{and u}k_{of Eq.}

共24兲, we find a better approximation in the 共k+1兲th iteration. Equation共24兲 becomes

0 = G_ᐉ共Vk+1兲 + 共⌬t/⌬x兲uk+1F_ᐉ共Vk+1兲 − R_ᐉm,m−1. Substituting Vk+1_=Vk_{+ dV,u}k+1_{= u}k_{+ du, and neglecting}

second-order terms we obtain

冉

⳵G_ᐉ ⳵V 共V k_{兲 +} ⌬t ⌬xu k⳵Fᐉ ⳵V共V k_兲

冊

_{dV +}⌬t ⌬xFᐉ共V k_{兲du = − R} ᐉ k ,

where we have defined R_ᐉk as

R_ᐉk

= G_ᐉ共Vk_{兲 + 共⌬t/⌬x兲u}k_F

ᐉ共Vk兲 − Rᐉm,m−1.

This is solvable for共dV,du兲 if u⌬t/⌬x is not a character-istic speed, which can be achieved by taking ⌬t small enough. After division of this equation by uk⌬t/⌬x, we ob-tain the following linear system to be solved at each Newton iteration, written in the notation of Eq.共17兲:

M

冉

V,− ⌬x uk_⌬t

冊

冢

dV1 dV2 dV3 du/uk

冣

=

冢

−R₁k −R₂k −R₃k −R₄k

冣

. 共25兲 C. Numerical implementation

specify the fluxes of all components at the injection bound-ary. Initially, all cells contain a homogeneous distribution of water and an oleic phase at low temperature.

All calculations in the Newton-Raphson scheme depend on the old phase condition of cell m, as well as on available information from cells to the left of cell m. The method of solution depends on the new condition of the cell, two-phase or three-phase.

The iterative procedure is simple for a cell that starts and stays in the same condition. When a cell starts in the two-phase condition, but in the two-two-phase calculation a tempera-ture arises that exceeds the boiling temperatempera-ture of the water and oleic-phase mixture 关see Fig. 2 共left兲兴 then the calcula-tion is replaced by a three-phase calculacalcula-tion.

Simulations use a uniform grid with 2000 blocks. This implicit method is inexpensive as it only involves the solu-tion of many 3⫻3 matrices as opposed to a single big ma-trix. As far as the authors are concerned, this upstream box finite difference method is original in the way the total ve-locity is treated, as there is no time derivative for it in the system.共See Ref. 关38兴 for a related scheme.兲

VI. METHOD OF CHARACTERISTICS FOR STEAM INJECTION

Figure 3 compares the numerical solutions obtained by the current finite-difference scheme共FD兲 and by the method of characteristics共MOC兲 used in Ref. 关13兴 for the saturations 共Sw, So, Sg兲 versus the length along the cylinder for pure

steam injection in a cylinder filled initially with dead oil only. The profiles are shown after the injection of 0.057 P共ore兲 V共olume兲 共cold water equivalent兲. In region I 共the steam zone兲, where Sg⬎0, we observe a saturation

rarefac-tion wave, in which the temperature and Darcy velocity are constant. At the steam condensation shock or front 共SCF兲, where the temperature drops to the initial temperature, the gas saturation drops to zero. Here we use the word shock and

front interchangeably. The water saturation is larger than the initial water saturation共Swc= 0.15兲, both in regions I and II.

Downstream of region II there is a second shock to the initial conditions, i.e., Sw= Swc, So= 1 − Swc, in region III. The total

downstream Darcy velocity divided by the injection velocity is constant spatially at 1.19⫻10−3 _{in the entire liquid zone} 共regions II and III兲, but it shows numerical fluctuations of 20% between time steps. Nevertheless the average is correct. The oil saturation at the SCF is about 0.3. The observed behavior is approximately independent of the number of grid blocks. The reduced temperature is plotted, but the total Darcy velocity is not because they visually coincide.

VII. RESULTS

We distinguish two classes of results. In the first class the volatile oil is initially present in the reservoir, but it is not coinjected. In the second class no volatile oil is initially present, but it is coinjected with the steam. For each of the classes the volatile oil is cyclobutane, heptane, or dodecane. These alkanes were chosen because they have low, medium, and high boiling points, in such a way that each one gives rise to a different type of solution.

The initial conditions for all the calculations are the fol-lowing. The initial temperature is 293 K and the gas satura-tion is zero. The initial water saturasatura-tion is given as Sw= Swc

= 0.15. We consider the cases where the oleic phase consists of dead oil and of a volumetric 50% mixture of dead oil and volatile oil; however, for cyclobutane we use a volumetric 20% mixture of dead oil and volatile oil, as 50% is above the solubility limit. In the former case we inject an alkane and steam vapor with mass fraction 0.2 关alkane/共alkane + steam兲兴. In the latter case we displace with pure steam. The injection temperature is 373 K and the injection pressure is one atmosphere. We use atmospheric pressure because these results are easiest to validate by laboratory experiments. The volumetric injection flux is 9.52⫻10−4_{m / s. From now on,} all figures plot reduced quantities versus the distance. The reduced velocity共u兲 is the total Darcy velocity divided by the injection velocity. The reduced temperature is 共T − To兲/共Tb

w

− To兲 共see Appendix B for terminology兲. The

re-duced concentration共volume fraction兲 is vov=␳ov/␳V.

Long time runs were subdivided in shorter ones to dampen transient behavior faster; in each one the initial data consisted of the previous one, where every other grid data was omitted关39兴. Each run was stopped before breakthrough of the fastest wave, i.e., before it reached the end boundary.

A. Cyclobutane and steam mixture displacing dead oil Figure 4 shows displacement of dead oil by steam and cyclobutane. There are four regions from the injection point to the initial condition. In region I, there is a fast three-phase saturation rarefaction wave with speed given by Eq. 共19兲 共right兲; T and u are constant on this wave; the volatile oil concentration 共vov兲 is a small constant. The SCF separates

region I from region II. In region II the temperature and the total velocity are constant, but much lower than in region I. Regions I and II are three-phase regions, whereas regions III

FIG. 3. Comparison of the MOC with the FD solutions for the case that volatile oil is neither injected nor present in the initial oil. The curves with sharp edges are obtained with MOC, and the smoother ones with FD. The saturation curves are indicated as

S_w, S_g, S_oin obvious notation. The temperature is indicated with a

and IV are two-phase liquid regions. Regions II and III are separated by a cyclobutane evaporation shock, where the temperature jumps to its initial value in the reservoir and the total velocity to its final downstream value. There is a con-stant state in region III. Regions III and IV are separated by a Buckley-Leverett shock. Region IV contains the initial saturations.

The rarefaction in region I starts at the injection state 共Swc, 1 − Swc, Tb

w

, uinj兲 and ends at the left state of the SCF, 共Sw,SCF− , Sg,SCF− , Tb

w

, uinj_{兲. The right state of the SCF is SCF}+ =共S_w,SCF+ , S_g,SCF+ , T_SCF+ , u_SCF+ 兲. Left and right states and vSCF satisfy the RH conditions共10兲–共13兲. The velocity vSCFis the same as the speed of fast three-phase rarefaction关Eq. 共19兲兴 at the end of region I, i.e., the SCF shock is left-characteristic. Region II starts at SCF+_{, which is also the upstream 共left兲} state of the cyclobutane condensation shock with speedvC.

The right state of this shock is C+=共S_w,C+ , S_g,C+ = 0 , To, uC

+_兲. Left and right states and vC satisfy the RH conditions

共10兲–共13兲. Region III is a constant state. Therefore C+_{is the} upstream共left兲 state of the Buckley-Leverett shock. Region IV is a constant state, which is the right state of this shock, with initial reservoir saturation and temperature.

B. Pure steam displacing a dead-oil and cyclobutane mixture As shown in Fig.5, there are four regions again. In region I, there is a fast three-phase saturation rarefaction wave at constant T and u. Separating region I from region II, there is the SCF. Note that the temperature at the right side of the SCF is not the initial reservoir temperature, but an interme-diate temperature. Region II consists of a constant state with temperature and total velocity lower than in region I. Re-gions II and III are separated by a three-phase saturation shock, which does not change the temperature but reduces the water saturation to its initial value. The gas saturation in region III is slightly lower than in region II. Between region III and region IV there is a cyclobutane condensation shock. Region I 共vov= 0兲 starts at the injection state 共Swc, Sg= 1

− Swc, Tb w

, uinj兲 and ends at 共S_w,SCF− , S_g,SCF− , Tb w

, uinj兲, which is

the left state of the SCF. The right state of the SCF is SCF+=共S_w,SCF+ , S_g,SCF+ , T_SCF+ , u_SCF+ 兲. Left and right states and

vSCF satisfy the RH conditions 共10兲–共13兲. The SCF is left characteristic. Region II consists of the constant state SCF+. Region II ends at the three-phase saturation shock with speed

vS. The right state is denoted as 共Sw,S

+ _{, S} g,S + _{, T} SCF + _{, u} SCF + _{兲 and} continues in region III. This constant state ends at the con-densation shock with speedvC. Region IV is a constant state,

which is the right state of this shock, with initial reservoir saturation and temperature.

C. Steam and heptane mixture displacing dead oil As shown in Fig.6, in region I there is again a fast three-phase saturation rarefaction. At the SCF the temperature drops to the initial temperature, and a volatile-oil bank 共re-gion II兲 builds up downstream of the SCF. The volatile-oil bank does not contain any dead oil. Such a pure volatile-oil bank displaces all dead oil. Downstream of the volatile-oil

FIG. 4. Steam and cyclobutane displacement of dead oil. The water saturation S_w, oil saturation S_o, and steam saturation S_gare shown as the dashed, solid, and dashed-dotted curves. The volatile-oil fractionv_ov in the oil phase and the reduced temperature are shown as the dotted and filled-square curves.

FIG. 5. Steam displacement of a dead-oil or cyclobutane mix-ture. The water saturation, oil saturation, and steam saturation are shown as the dashed, solid, and dashed-dotted curves. The volatile oil fraction in the oil phase and the reduced temperature are shown as the dotted and filled-square curves.

bank there is a contact wave, which marks the boundary between regions II and III. The contact wave is smooth in the simulation due to numerical diffusion. Downstream there is only dead oil. Region III consists of a constant state. Regions III and IV are separated by a Buckley-Leverett shock.

Region I starts at the injection state 共Swc, Sg= 1

− Swc, Tinj, uinj兲, and it ends at the left state of the SCF, viz.,

共Sw,SCF

−

, S_g,SCF− , Tinj, uinj兲. Again the SCF is left characteristic. The right state of the SCF is SCF+=共S_w,SCF+ , S_g,SCF+ ,vov

= 1 , To, u+兲 and continues as a constant state in region II.

Between regions II and III there is a volatile-oil contact wave with right state SCF+ _{and velocity v}

␷. Downstream of the contact wave the constant state 共S_w,C+ , S_g,C+ ,vov= 0 , To, u+兲

spans region III. The Buckley-Leverett shock separates re-gion III from rere-gion IV. This solution agrees with the obser-vation made previously by Bruining and collaborators in the laboratory experiments关5,7兴. This is the case when analyti-cal, numerianalyti-cal, and experimental results are all available. They all agree.

D. Pure steam displacing a dead-oil and heptane mixture As shown in Fig.7, there are only three regions. In region I, there is the usual rarefaction wave with constant T and u. Separating region I from region II there is the SCF. In region II the temperature is equal to the initial temperature and the total velocity attains its constant downstream value. At the SCF there is a remarkable spike of volatile oil. The volatile oil concentration vov vanishes at the left of the SCF, it

reaches almost one at the SCF, and then it declines to its initial value. Region II consists of a constant state. Regions II and III are separated by a Buckley-Leverett shock. Region III contains the initial saturations.

Region I starts at the injection state 共Swc, Sg= 1

− Swc, Tb w

, uinj兲 and ends at the left state of the SCF, viz., 共Sw,SCF− , Sg,SCF− , Tb

w

, uinj_{兲. The right state is SCF}+ =共S_w,SCF+ , S_g,SCF+ = 0 , To, u+兲. Left and right states and vSCF sat-isfy the RH conditions共10兲–共13兲. Again the SCF is left char-acteristic. Region II consists of the constant state SCF+. A

Buckley-Leverett shock with velocityvBLseparates region II

from region III.

Let us discuss the evolution of the volatile-oil bank. Ini-tially there is no volatile-oil bank. It starts to be formed after injection. It grows as long as the volatile-oil–dead-oil mix-ture reaches the steam zone. The growth of the condensed volatile-oil bank stops as soon as no more volatile oil, carried by the liquid oil, can reach the steam zone.

Now we compare the case in Sec. VII C of heptane coin-jection to the case in Sec. VII D where a mixture of dead oil and heptane is displaced by pure steam. For the coinjected 共20%兲 case a large oil bank is built up. On the other hand, with as much as 50% oil in the initial oil mixture the volatile-oil bank is very small. This is so because only the volatile oil that is stripped from the dead oil enters the steam zone and contributes to the building up of the volatile-oil bank for case共b兲. In particular, when the initial volatile-oil fraction is small it can take a long while before such a bank is built up, whereas such a building up is much faster with coinjected volatile oil.

E. Steam and dodecane mixture displacing dead oil As shown in Fig.8, besides the fast three-phase saturation rarefaction in region I, there is a volatile-oil condensation wave in region II with velocity共20兲; in this wave both T,u vary. At the state C joining rarefactions in regions I and II the two characteristic speeds coincide. Region II is separated from region III by the SCF. Region III consists of a constant state, which is separated from region IV by a Buckley-Leverett shock.

Region I starts at the injection state 共Swc, Sg= 1

− Swc, Tinj, uinj兲 and it ends at the coincidence point C

=共Sw,C, Sg,C, Tinj, uinj兲, the left state of the condensation

rar-efaction. The three-phase rarefaction wave in region I is con-tinued in region II as a condensation rarefaction wave, which is connected to the SCF. The left state共S_w,SCF− , S_g,SCF− , T_bw, u−_兲, the right state共S_w,SCF+ , Sg= 0 , To, u+兲, and vSCFsatisfy the RH

conditions共10兲–共13兲. The SCF is left characteristic.

FIG. 7. Steam displacement of a dead-oil and heptane mixture. The water saturation, oil saturation, and steam saturation are shown as the dashed, solid, and dashed-dotted curves. The volatile-oil frac-tion in the oil phase and the reduced temperature are shown as the dotted and filled-square curves.

F. Pure steam displacing dead-oil and dodecane mixture As shown in Fig. 9, there are five regions. In region I, there is a three-phase saturation rarefaction. Near the injec-tion point and at some other points very small transient ef-fects are observed in the simulation. Region II consists of a constant state starting approximately at distance 0.1 m. Sepa-rating region II from region III there is a dodecane evapora-tion shock, followed by a fast composievapora-tion rarefacevapora-tion wave. The evaporation shock speedvE coincides with the speed of

the left part of the composition rarefaction wave. In region III the temperature and the total velocity are lower than in region I. Regions III and IV are separated by the SCF, which is left characteristic. The temperature drops to its initial res-ervoir value and the total velocity to its final downstream value. There is only a constant state in region IV. Regions IV and V are separated by a Buckley-Leverett shock with speed

vBL. Region V contains the initial saturations.

Region I starts at the injection state 共Sw, Sg, T , u兲

=共Swc, 1 − Swc, Tb w

, uinj_{兲 and finishes at 共S} w,E − _{, S} g,E − _{, T} b w , uinj_兲.

This state also represents the left side of the dodecane evapo-ration shock as region II is a constant state. The dodecane evaporation shock has speed vE and the right state E+

=共S_w,E+ , S_g,E+ , TE+, uE+兲. Left and right states and vE satisfy the

RH Eqs.共10兲–共13兲. The evaporation shock is left character-istic too. Region III starts at E+_{with a composition} rarefac-tion, which ends at共S_w,SCF− , S_g,SCF− , TE

− , uE

−_{兲, the left state of the} SCF.

The right state of the SCF is SCF+_=共S

w,SCF

+ _{, S}

g,SCF

+ = 0 , To, uSCF

+ _{兲. Left and right states and vSCF satisfy the RH} conditions. The SCF is left characteristic. Region IV consists of the constant state SCF+and ends with a Buckley-Leverett shock.

G. Comparison of cases

We can distinguish three important mass transfer waves, viz., the SCF, the evaporation wave upstream of the SCF, and the condensation wave downstream of the SCF. Inspection of the results reveals the crucial role of the speeds of these

waves, i.e., the speeds of the evaporation shock, the steam condensation shock, and the condensation wave 共shock or rarefaction兲. For the medium boiling-point alkane 共heptane兲 the three waves merge into a single wave, leading to high recovery. Both high recovery 关5兴 and the existence of the volatile oil bank had been observed in the coinjection experi-ments关7兴. Similar observations can be made for Willman’s experiments关9兴, where medium boiling temperature oil was present initially. For the other cases the three waves spread out leading to lower recovery. These statements hold true irrespective of whether the volatile oil is present initially or coinjected.

For the high boiling-point alkane共dodecane兲 the conden-sation wave collapses on the SCF, whereas the evaporation wave separates from the SCF. This leads to some positive effect on the oil recovery.

For the low boiling-point alkane共cyclobutane兲 the evapo-ration wave collapses on the SCF, whereas the evapoevapo-ration wave separates from the steam condensation shock. This has only a small effect on the oil recovery.

In summary, resonance, i.e., equality of wave speeds, leads to high amplitude waves, i.e., favorable recovery. This occurs for medium boiling-point alkanes.

We use these ideas to find the bifurcation loci in the pres-sure or carbon-number plane. We used an injected mass frac-tion 共mass volatile-oil or total-mass兲 of 20%. We derived polynomial expansions for the properties of volatile oil in terms of the carbon number, i.e., for the viscosity, the liquid heat capacity, the evaporation heat, the molar weight, the liquid density, and the boiling point. In this way we find not only the properties of the alkanes with integer carbon num-bers but also of any pseudocompound, characterized by any real value for the carbon number. We carried out enough simulations to isolate the wave sequences typical of high, medium, or low boiling-temperature volatile oils. The result-ing curves are shown in Fig. 10. The curves are accurate within a carbon-number change of 0.1. The range of favor-able medium temperature boiling point, where resonance oc-curs, becomes larger with increasing pressure.

In our simulation we injected 20 w / w % volatile oil in the steam rather than 5% as we did in the laboratory experiments for reasons of clear illustration. For 5% simulations the volatile-oil bank is thinner as expected, but the overall pic-ture does not change共see Fig.11兲.

FIG. 9. Steam displacement of a dead-oil or dodecane mixture. The water saturation, oil saturation, and steam saturation are shown as the dashed, solid, and dashed-dotted curves. The volatile-oil frac-tion in the oil phase and the reduced temperature are shown as the dotted and filled-square curves.

VIII. CONCLUSIONS

We developed a model that captures the main physical features of thermal three-phase flow, involving water, dead oil, and volatile oil. The numerical solution for different in-jected mixtures and initial oil composition reveals its struc-ture in terms of rarefaction and shock waves. These waves are validated by verifying that they satisfy all properties pre-dicted by mathematical analysis based on the mathematical theory of nonlinear conservation laws. In the solution found computationally numerical diffusion effects are controlled by using extremely small grid cells.

In the 1D setting, coinjection of medium boiling tempera-ture volatile oil in steam leads to 100% recovery of oil共Fig. 6兲. This improvement is due to the formation of an increas-ingly long volatile-oil bank displacing the oil in place. The initial presence of medium boiling-temperature volatile oil also improves oil recovery共see Fig. 7兲. This is due to the formation of a thin volatile-oil bank displacing the oil in place. Clearly the volatile-oil bank displaces all the dead oil because they are in the same phase. This solution agrees with

the observations found previously in the laboratory experi-ments. There is agreement between analytical, numerical, and experimental results.

As far as the recovery efficiency is concerned, the initial presence of medium boiling-temperature volatile oil has a positive effect. The initial presence of high boiling-temperature volatile oil has a much smaller effect 共see Fig. 9兲. Coinjection of high boiling-temperature volatile oil in steam has a negligible effect共see Fig.8兲. Coinjected or ini-tially present low boiling-temperature volatile oil has no ef-fect共see Figs.4and5兲.

This model reveals that the essential mechanism for good recovery is that all the volatile oil condenses at the same point where the steam condenses. In mathematical language, this occurs when the evaporation and condensation shock speeds coincide with the speed of the steam condensation front. In physical language, high boiling-temperature volatile oil finds it difficult to evaporate and therefore the evapora-tion wave is slower than the steam condensaevapora-tion shock. The low temperature boiling-point volatile oil finds it difficult to condense and therefore the condensation wave is faster than the steam condensation shock. For the medium oil the evapo-ration and the condensation wave collapse on the steam con-densation shock.

ACKNOWLEDGMENTS

This work follows the brilliant inspiration by Daan N. Dietz. Experimental evidence was also gathered in the MSc. theses of A. Emke, G. Metselaer, J. W. Scholten, D. W. van Batenburg, R. Quak, and C. T. S. Palmgren. This paper was initiated at the Institute for Mathematics and its Applications in Minnesota with NSF support and completed with support from NWO-WOTRO 共Stichting voor Wetenschappelijk Onderzoek van de Tropen兲. We thank C. J. van Duijn for his early participation in this effort. We are grateful to our librar-ian, Paul Suijker, for discovering details of the references. We thank the anonymous referees for valuable comments. We thank TU-Delft and IMPA for their hospitality and sup-port.

APPENDIX A: STEAM OR HEPTANE MIXTURE DISPLACING A DEAD OIL

WITH LESS VOLATILE OIL

To illustrate the effect of low共5 w/w %兲 concentration of volatile oil in the injected mixture we made the run shown at the bottom of Fig.11. Clearly the wave structure is identical to that found in Fig.6. The difference is only that the dis-solved volatile-oil peak is much thinner. Even for this low-injection concentration the recovery is high, as shown by the small amount of oil left behind. However, we observe that the initialization effect, i.e., the transient bump of oil left behind near the injection point is larger because the volatile-oil bank takes longer to build up.

The top part of Fig.11includes the fluxes for water, vola-tile oil, and dead oil, and also the pressure. We rescaled the values for clear illustration. However, all mass fluxes are on

the same scale. Notice that the pressure gradients are much larger in the liquid zone.

APPENDIX B: PHYSICAL QUANTITIES—SYMBOLS AND VALUES

In this appendix we summarize the values and units of the various quantities used in the computation and empirical ex-pressions for the various parameter functions. All enthalpies per unit mass are are taken relative to the reference tempera-ture of the components in their standard form. All heat ca-pacities are specified at constant pressure. All enthalpies in their standard form are zero at the reference temperature.

1. Temperature-dependent properties of steam, water, and heptane

We use Refs.关40–42兴 to obtain all the following proper-ties. Properties for other volatile components, such as cy-clobutane and dodecane can be derived from the same refer-ences. A conventional choice for the reference temperature is

T

¯ =298.15 K.

The rock enthalpy Hris

共1 −␸兲Cr共T − T¯兲, Cr= 3 274 000共J/m3K兲. 共B1兲

The liquid water enthalpy hW共T兲 关J/kg兴 is approximated

by

hW共T兲 = cW共T − T¯兲, cW= 4184共J/kg K兲. 共B2兲

The heptane enthalpy hoV 关J/kg兴 and the dead-oil enthalpy

hoD are approximated by

hoV共T兲 = coV共T − T¯兲, hoD共T兲 = coD共T − T¯兲. 共B3兲

The values for the heat capacities of heptane and dead oil are

coV= 2242, coD= 1914.1共J/kg K兲. 共B4兲

The liquid volatile-oil and dead-oil heat capacities are not the exact heat capacities, but slightly adapted so that the en-thalpy of oil per unit volume is independent of composition. Therefore the heat capacity of the oleic phase per unit vol-ume can also be defined independently of composition, lead-ing to an oleic phase heat capacity per unit volume of

Co= 1.531⫻ 106共J/m3K兲.

The steam enthalpy hgWis given by

hgW共T兲 = hgW s _{共T兲 + ⌳}

W共T¯兲, 共B5兲

and the sensible steam enthalpy is approximated as

hgW s _{共T兲 = c}

pgw共T − T¯兲, cpgw= 1964共J/kg K兲. 共B6兲

The volatile-oil vapor enthalpy hgVas a function of

tem-perature is given by

hgV共T兲 = hgV s _{共T兲 + ⌳}

V共T¯兲 共B7兲

and the sensible heptane enthalpy is approximated as

hgV s _{共T兲 = c}

pgv共T − T¯兲, cpgv= 1658共J/kg K兲. 共B8兲

For the latent heat ⌳W共T兲 关⌳V共T兲兴 共J/kg兲 or evaporation

heat of water共heptane兲 we use

⌳W= 3 105 600 − 2220T, ⌳V= 538 830 − 584T.

共B9兲 The liquid-water viscosity␮w关Pa s兴 is approximated by

␮w= exp共− 12.06 + 1509/T兲. 共B10兲

The viscosity of the dead oil␮odand heptane␮ovare written

as

␮od= e−13.79+3781/T, ␮ov= e−10.813+880.2/T, 共B11兲

and the viscosity of the oleic phase is approximated by the quarter power rule

␮o=

冉

␳ov ␳V ␮ov 1/4 +␳od ␳D ␮od 1/4

冊

4 . 共B12兲

We assume that the viscosity of the gas is independent of composition

␮g= 1.8264⫻ 10−5共T/300兲0.6. 共B13兲

The water saturation pressure is given by Eq. 共3兲. The pure phase densities of steam and volatile-oil vapor are given by Eqs.共2兲 and the corresponding concentrations␳gw,␳gvare

given in Table.I

2. Three-phase relative permeabilities

We used Stone’s expressions关31兴 for three-phase perme-ability: Eqs.共B14兲 and 共B15兲 describe the water relative per-meability krw, the gas-phase relative permeability krg, and the

oil relative permeability kro, respectively. For convenience

we have taken the residual oil parameter Somused by Fayers

关31兴 equal to zero. The relative permeabilities krw, krg are

functions solely of the water saturation Swand the gas

satu-ration Sg, respectively. krw= krw

⬘

Swe3+2/␭, krg= krg

⬘

共1 − Sge兲2共1 − Sge1+2/␭兲, 共B14兲 kro= So共1 − Swc兲 krcow共1 − Sw兲共1 − Swc− Sg兲 krowkrog, Swe= Sw− Swc 1 − Swc− Sor , Sge= 1 − Sg− Swc 1 − Swc− Sor , krow= krg

⬘

共1 − Swe兲2共1 − Swe 1+2/␭_{兲, k} rog= krw

⬘

Sge 3+2/␭_. 共B15兲 We took k_rw

⬘

= 1 / 2, k_rg

⬘

= 1, and krcow= 1. Here ␭=0.5 is the

sorting factor, Swcgiven in the table and Sor= 0 are the

con-nate water saturation and the residual oil saturation, respec-tively.

f_␣=共kr␣/␮␣兲/共krw/␮w+ kro/␮o+ krg/␮g兲, 共B16兲

where f_␣ is the fraction of the volume flux of phase␣ 关26兴. APPENDIX C: ANALYSIS OF THE CHARACTERISTIC

EQUATIONS

We observe that making the determinant of Eq.共17兲 equal to zero leads to a polynomial equation of third order in ␭ =␩/ u and thus we get 共if all solutions are real兲 a slow, a medium, and a fast wave solution.

We want to eliminate the first two elements of the fourth row of Eq. 共17兲. To do so we find x and y such that xF11 + yF31+ F41= 0 with Fijgiven in Eq. 共D4兲 and we obtain

y = Hg− Ho+ ␳gw ␳W 共Ho− Hw兲 1 −␳gv ␳V −␳gw ␳W , x = Hg− Hw− ␳gv ␳V 共Ho− Hw兲 1 −␳gv ␳V −␳gw ␳W .

Remembering that Gij are given in Eqs. 共D3兲, it is easy to

verify for␣= 1 , 2,

TABLE I. Summary of physical input parameters and variables.

Physical quantity Symbol Value Unit

Water, gas, oil fractional flows fw, fg, fo Eq.共B16兲. m3/ m3

Steam, volatile-oil enthalpy/unit mass hgW, hgV Eqs.共B5兲 and 共B7兲 J/kg

Sensible enthalpy/unit mass h_gWs , h_gVs Eqs.共B6兲 and 共B8兲 J/kg

Volatile-oil, oil enthalpy/unit mass hoV, hoD Eqs.共B3兲 J/kg

Gas enthalpy, oil enthaply H_g, H_o H_gw+ H_gv, H_ov+ H_od J / m3

Steam, volatile-oil enthalpy HgW, HgV ␳gW共T兲hgW共T兲,␳gV共T兲hgV共T兲 J / m3

Sensible steam, volatile-oil enthalpy H_gWs , H_gVs ␳gW共T兲hgWs 共T兲,␳gV共T兲hgVs 共T兲 J / m3

Partial steam, volatile-oil enthalpy Hgw, Hgv ␳gw共T兲hgW共T兲,␳gv共T兲hgV共T兲 J / m3

Pure volatile-oil, dead-oil enthalpy H_oV, H_oD ␳V共T兲hoV共T兲,␳oD共T兲hoD共T兲 J / m3

Volatile-oil, dead-oil enthalpy Hov, Hod ␳ov共T兲hoV共T兲,␳od共T兲hoD共T兲 J / m3

Rock enthalpy H_{r Cr共T−T¯兲, Eq. 共B1兲} J / m3

Water enthalpy HW ␳W共T兲hW共T兲 J / m3

Porous rock permeability k 1.0⫻10−12 _m2

Water, gas, oil relative permeabilities krw, krg, kro Eqs.共B14兲 and 共B15兲 m3/ m3

Molar weight, H2O, C7H16, dead-oil MW, MV, MD 0.018, 0.10021, 0.4 kg/mole

Total pressure P 1.0135⫻105 Pa

Atmospheric pressure Po 1.0135⫻105 Pa

Partial pressures P_w, P_v Eqs.共3兲 and 共4兲 Pa

Water, vapor, oil saturations Sw, Sg, So independent variables m3/ m3

Residual oil, connate water saturation S_or, S_wc 0, 0.15 m3_{/ m}3

Injection saturations S_winj, S_oinj input m3/ m3

Temperature T independent variable K

Three-phase temperature T共xov= 1兲 Eq.共5兲 K

Reservoir, injection temperature T_o, Tinj _{293, 370–373 K}

Boiling point of water, volatile oil T_bw, T_bv 373.15, 371.57 for heptane K

Total Darcy velocity u volume flux of all phases m3_/共m2_s兲

Total injection velocity uinj injected volume flux m3/共m2s兲

Water, volatile-oil evaporation heat ⌳W,⌳V see Eqs.共B9兲 J/kg

Water, steam, oil viscosity ␮w,␮g,␮o Eqs.共B10兲–共B13兲 Pa s

Water, steam, volatile-oil vapor density ␳W,␳gW,␳gV 998.2, Eqs.共2兲 kg/ m3

Pure heptane, dead-oil densities ␳V,␳D 683, 800 kg/ m3

Steam, volatile-oil vapor concentrations ␳gw,␳gv ␳gWPw/ P ,␳gVPv/ P, Eq.共2兲 kg/ m3

Liquid volatile-oil concentrations ␳ov,␳od obtained from Eqs.共6兲 and 共1兲 kg/ m3

Molar fraction volatile-oil in liquid-oil xov Eqs.共5兲 and 共4兲 ⫺

xF12+ yF32+ F42= 0, xG1␣+ yG3␣+ G4␣= 0. NowF43= xF13+ yF33+ F43in the system 共17兲 becomes

F43=⌬HTfg+ Ho

⬘

+共Hw

⬘

− Ho

⬘

兲fw+共Hg

⬘

− Ho

⬘

兲fg, where ⌬HT= x

冉

␳gw ␳W

冊

⬘

+ y

冉

␳gv ␳V

冊

⬘

.

The primes indicate differentiation relative to temperature. Analogously,G43= xG13+ yG33+ G43 becomes Hr

⬘

+␸兵⌬HTSg+关Ho

⬘

+共Hw

⬘

− Ho

⬘

兲Sw+共Hg

⬘

− Ho

⬘

兲Sg兴其. FinallyF4= xF1+ yF3+ F4 becomes F4= Hg− ␳gw ␳W Hw− ␳gv ␳V Ho 1 −␳gv ␳V −␳gw ␳W .

Now we are ready for the Gaussian elimination. In Eq.共17兲 we add to the fourth row the first row multiplied by x and the third row multiplied by y leading to

冢

F11−␭G11 F12−␭G12 F13−␭G13 F1

F21−␭G21 F22−␭G22 F23−␭G23 F2

F31−␭G31 F32−␭G32 F33−␭G33 F3

0 0 F43−␭G43 F4

冣

. 共C1兲

The rest of the Gaussian elimination is just as tedious and straightforward. The result is the matrix on the left-hand side of Eq.共18兲, where A13共␭兲 ⬅ A₁₃0 −␭␸A₁₃1 = fw

⬘

+

冢

冉

␳gw ␳W

冊

⬘

−␳gw ␳W

冉

␳gv ␳V

冊

⬘

+

冉

␳gw ␳W

冊

⬘

␳gv ␳V +␳gw ␳W − 1

冣

共fg−␭␸Sg兲, A33共␭兲 =

冢

冉

␳gv ␳V

冊

⬘

+

冉

␳gw ␳W

冊

⬘

␳gv ␳V +␳gw ␳W − 1

冣

共fg−␭␸Sg兲 + fg

⬘

, A3= 1 ␳gv ␳V +␳gw ␳W − 1 + fg, A₁= fw+ ␳gw ␳W 1 −␳gv ␳V −␳gw ␳W . 共C2兲 Usingvod= 1 −vov we write A₂₃0 =vodfg

冉

␳gv ␳V

冊

⬘

冉

1 −␳gw ␳W

冊

+

冉

␳gw ␳W

冊

⬘

␳gv ␳V 1 −␳gv ␳V −␳gw ␳W −␷_ov

⬘

fo, 共C3兲 A₂₃1 =vod␸Sg

冉

␳gv ␳V

冊

⬘

冉

1 −␳gw ␳W

冊

+

冉

␳gw ␳W

冊

⬘

␳gv ␳V 1 −␳gv ␳V −␳gw ␳W +␸␷_ov

⬘

So, 共C4兲 A2= F2+vovA1− A3

冉

␳gv ␳V −␷o,v

冊

. 共C5兲

These results are used in Sec. IV B 1

APPENDIX D: FLUX, ACCUMULATION FUNCTIONS, AND DERIVATIVES

We use Eqs.共8兲 and 共7兲. We chose to divide the first and second of Eqs.共8兲 by␳Wand␳V, respectively. Thus we

ob-tain Eq.共9兲. The accumulation functions Gitake the

follow-ing form for water, volatile oil, total oil, and energy, respec-tively, Gw=␸

冉

Sw+ ␳gw ␳W 共T兲Sg

冊

, Gv=␸␷ov共T兲共1 − Sw− Sg兲 + ␳gv ␳V 共T兲Sg, Go=␸

冋

1 − Sw+

冉

␳gv ␳V 共T兲 − 1

冊

Sg

册

, GT= Hr+␸关Ho+共Hw− Ho兲Sw+共Hg− Ho兲Sg兴, 共D1兲

where all enthalpies are functions of the temperature. Simi-larly, the flow functions Fitake the form

Fw= fw+ ␳gw ␳W 共T兲fg, Fv=␷ov共T兲共1 − fw− fg兲 + ␳gv ␳V 共T兲fg, Fo=共1 − fw− fg兲 + ␳gv ␳V fg, FT= Ho+共Hw− Ho兲fw+共Hg− Ho兲fg. 共D2兲