Reservoir Simulation of Cyclic Steam

Stimulation in the Cold Lake Oil Sands

C.I. Beattie, SPE, Esso Resources Canada Ltd.; T.C. Boberg, SPE, Exxon Production Research Co.;
and 0.5. McNab, Esso Resources Canada Ltd.

Summary. Steam injectivity during cyclic steam stimulation (CSS) at Cold Lake can be achieved only by injecting at pressures high
enough to fail the formation mechanically. The reservoir also exhibits water/oil relative permeability hysteresis. This paper describes
enhancements made to a thermal reservoir simulator to incorporate these Cold Lake physics. The geomechanical model allows both
localized fracturing and global reservoir deformation (dilation and history-dependent recompaction). This representation allows the simulator
to match injection and production pressures that are otherwise difficult to reproduce. The history-dependent relative permeability hyste-
resis model calculates gridblock relative permeabilities that always lie on or between input imbibition and drainage bounding curves.
The model makes it possible to use laboratory-derived relative permeabilities in a simulation and still match field WOR's.

Introduction
The Cold Lake oil-sand deposit in northern Alberta contains over Dilation. Two observations lead to the conclusion that significant
40x 10 9 m 3 of bitumen in place. The extremely high oil viscosity dilation of the reservoir occurs during steam injection at Cold Lake.
and low native water mobility in this reservoir result in negligible First, surface benchmark arrays at Cold Lake have been used to
initial injectivity or productivity. Steam injection at commercial rates measure surface uplifts as great as 45 cm during injection, far larger
requires injection pressures high enough to cause both localized than can be attributed to thermal expansion or tensile fracturing
fracturing and widespread PV increases in the formation. The re- of the formation. Second, steam injectivity at Cold Lake is greater
sulting complex geomechanical behavior determines the initial steam than might be expected from native reservoir properties. A com-
flow paths and provides significant drive energy during CSS oper- mon problem when simulating steam injection into tar sands is
ations. 1 The reservoir also exhibits relative permeability hystere- matching relatively high observed injectivities when reasonable frac-
sis. These phenomena must be represented adequately in a reservoir
ture lengths and rock compressibilities are used. To match observed
simulator to properly model Cold Lake reservoir performance dur-
injectivities, most published simulations of Cold Lake CSS7-9 have
ing CSS.
This paper describes the deformation, fracturing, and relative per- used a rock compressibility one or two orders of magnitude larger
meability hysteresis behavior in the Cold Lake Clearwater reser- than the actual value of about 1. 0 GPa -I at reservoir condi-
voir and the methods used by Esso Resources Canada Ltd. (ERCL) tions.lO,ll Coats et al. ,12 who used this high-compressibility
to include these phenomena in our simulators. The physics and method for a California heavy-oil reservoir, called it the "spongy-
modeling aspects of each of these phenomena are described, and rock" approach. The only published Cold Lake simulation that did
methods to determine values for the additional input data required not use an enhanced-compressibility approach required extremely
by the simulator are discussed. Examples illustrate how the model long fractures to match the observed injectivity. 13
enhancements improve the simulator matches of field observations. The major problem with the spongy-rock approach is that it pre-
The enhancements described have been incorporated into a com- dicts a steady increase in injection pressure with time, whereas the
mercially available thermal simulator 2 and a thermal version of field pressures increase rapidly initially and then level off for most
Exxon Production Research Co. 's MARS simulator. 3 Although ofthe cycle. This is illustrated in Fig. 1, which shows typical field
different in formulation, both are fully implicit models and give and simulation results for first-cycle injection. The spongy-rock
virtually identical results when used to simulate the same problem. simulation used a rock compressibility of 35 GPa - 1 .
Surface heave and the high observed injectivity can be explained
Reservoir and Operations Description by two mechanisms that increase porosity. First, oil sands demon-
ERCL's Cold Lake operations target the Clearwater formation, a strate nonlinear compressibility behavior. Fig. 2 shows how the
near-shore deltaic sand of Cretaceous Age at a depth of about 450 compressibility increases dramatically as the effective stress nears
m. The sands are thick, often >40 m, with a high net/gross thick- zero. 14 In addition, shear failure can be induced in the formation
ness ratio. Porosity ranges from 30 to 35 %, with oil saturations at sufficiently low effective stresses in the presence of anisotropic
that average 70% PV. At the initial reservoir temperature of 13°C, stresses. 14,15 Shear failure caused by increasing pore pressure re-
the oil viscosity is about 100 Pa' s; oil viscosity decreases to 0.002
sults in dilation of the pore system.
Pa's at 250°C.
ERCL uses CSS at Cold Lake. Steam is injected at 225 m 3 /d As steam is injected into a relatively incompressible reservoir,
for 30 to 60 days, followed by a production period of 120 to 400 pore pressure increases and the effective stress decreases. At pore
days, depending on cycle number. Average producing day oil rates pressures corresponding to low effective stresses, the formation
decline from> 20 m 3 /d in the first cycle to about 7 m 3 /d in the compressibility increases by about two orders of magnitude, and
eighth cycle; WOR's increase from 1 to 2 to more than 4 during shear failure may also occur. Both these phenomena result in dila-
the same period. Total daily production from almost 1,800 pilot tion, with rapid increases in porosity and permeability. This great-
and commercial wells averaged 14000 m 3 /d in 1989. More de- ly increases the steam injectivity, and injection pressures increase
tailed descriptions of pilot 4 and commercial 5,6 operations are slowly thereafter. The dilation is reflected by uplift of the ground
available. • surface above the injection location.
We model this phenomenon by specifying a dilation pressure,
Deformation Modeling PD, below which behavior is elastic and a low original compressi-
High-pressure injection into a tar-sand reservoir causes fracturing bility value is used. Above the dilation pressure, a higher dilation
(which is discussed later) and widespread PV changes (deforma- compressibility is used, enabling the porosity to increase rapidly
tion) in the formation. Deformation involves both dilation and sub- with pressure, as shown by the line labeled "dilation" in Fig. 3.
sequent recompaction. Details of thc modeling of these phenomena A maximum porosity is also specified, above which further dila-
are discussed below. tion is not permitted and the compressibility reverts to a low value.
Fig. 1 shows the improved match of injection pressure that results
Copyright 1991 Society of Petroleum Engineers from use of this dilation model.

200 SPE Reservoir Engineering, May 1991

 1.4
12r-----------------------------------~_,

1.2

11
... .. ...
----..
1.0
'0
0..
3 0.8
VI
VI
w
'"
I-
VI 0.6
SPONGY ROCK w
DILATION >
;::: 0.4
FIELD DATA u
w
u..
u..
w
0.2

0.0
10 20 30 40 0.1 10 100 1000
INJECTION TIME (DAYS)
COMPRESSIBILITY (GPa ")

Fig. 1-Slmulator match of first-cycle Injection pressure. Fig. 2-0ll-sand compresslblllty.14

The functional form used for all the porosity/pressure relation- Recompaction. Although the recompaction behavior of an oil sand
ships in our deformation model is that has undergone dilation is not well understood, the general be-
havior can be inferred. Surface heave at Cold Lake is cyclical, in-
cf>=cf>r exp[c(p-Pr)]' .............................. (1)
creasing during injection and decreasing during production.
where cf>r occurs at reference pressure Pn cf> is a porosity at some However, some of the surface uplift is permanent. The oldest area
pressure P, and c is compressibility. The linearized approximation of our Leming pilot, for example, has a permanent residual uplift
to Eq. 1 often used in reservoir simulators is not valid for the large of 15 cm that has remained unchanged for several years, although
compressibilities typical of dilation. injection and production have continued during that time. Hence,
Use of the dilation model requires estimates of the dilation pres- while partial reversal of the dilation occurred, recompaction to the
sure, PD, the dilation compressibility, and the maximum allowa- original porosity is not occurring.
ble porosity. The dilation pressure is close to or slightly lower than Field observations also suggest that recompaction does not be-
the fracture pressure. 14,15 The dilation compressibility can be es- gin immediately when the pressure begins to decline following di-
timated from Fig. 2 to be 100 to 1000 GPa -I. Experience shows lation. When a Cold Lake CSS well is initially placed on production
that results are relatively insensitive to variations within this range. following steaming, the bottomhole pressure (BHP) is sufficiently
We generally use a large value for the maximum allowable porosi- high that the well flows without pumping. After flowback ends at
ty and find that other regions of the reservoir become the a BHP of about 4 MPa, the pump is seated to continue production.
predominant failure zone before the most highly pressured grid- In a first-cycle well, flowback normally lasts for less than 2 weeks.
blocks reach their maximum porosity. Using this approach, we have During this time, the pressure declines from more than 10 to 4 MPa.
observed porosity increasing to a maximum of 110 to 120 % of the A thermal simulator cannot compute a pressure change of this mag-
original value, which compares favorably to a published predic- nitude with a relatively small fluid withdrawal unless an initial period
tion of 123 %. 15 of elastic behavior with little recompaction to support the pressure
A similar treatment of porosity increases with pressure during is assumed.
injection was proposed by Ito 16 as part of his "sand deformation Recompaction, therefore, has two phases: an initial elastic peri-
concept. " Ito limited the maximum porosity to 104% of the origi- od when no recovery of dilation occurs and porosity is changing
nal value, however, and did not assume the existence of a fracture only as a result of the original compressibility, and a recompaction
in his simulation. In our approach, both fracturing and deforma- period with an enhanced compressibility that allows recovery of
tion are assumed to occur because this best fits available field ob- some of the dilation that occurred during injection. We model this
servations and published theory. 14 Ito also assumed that dilation behavior as shown schematically in Fig. 3. After the pressure be-
was reversible. As the next section shows, we found that this is gins to decline but before it reaches the recompaction pressure, PR,
not a good assumption. the slope of the porosity/pressure function is determined by the low

jr--------- • INITIAL
°MAX " SAND
/ DEPOSITION
: CONSOLIDATION ,,'
ELASTIC :

1
CURVE,)""

7" " "

>-
I-
Ui END OF

~
0
'"0
,,-- ---------------------- §~ A
0..
END OF
PRODUCTION .

J . . ~=====~~~~- ------ ~-------'

A.
'1'0
8__ ----------

J INITIAL
RESERVOIR
CONDITIONS
...______...!E;:LA""S:o,T:..:.IC=--__- . J
Q

1 °0 POINT OF
MAXIMUM
COMPACTION

~
VIRGIN
RESERVOIR
CONDITIONS

INCREASING PORE PRESSURE - -

INCREASING EFFECTIVE STRESS
PRESSURE

Fig. 4-Representatlon of Cold Lake sand behavior basad on

Fig. 3-Reservolr deformation model. Ref. 17.

SPE Reservoir Engineering, May 1991 201

 7500.,---------------------------------,
12,-------------------------------------,

10
PR = 3.0 MPa ~

'"E
PR = 4.1 MPa
PR = 6.5 MPa ...
~7000

PR = 9.0 MPa o
'"3
fiELD DATA ;; 6500
...
CJ
~
C
o
g: 6000

5500~----~------~----~------~----~

0.0 0.2 0.4 0.6 0.8 1.0

30 60 90 120 150 RESIDUAL DILATION fRACTION. fR
PRODUCTION TIME (DAYS)

Fig. 6-Slmulator predictions of first-cycle produced fluid

Fig. 5-Slmulator match of flrst-cycle production pressure. volumes.

original compressibility, as illustrated by the upper line labeled The pressure at which recompaction begins, PR, andfR must be
"elastic. " The simulator tracks the maximum pressure and porosity obtained from history matching because no published estimates of
reached by each gridblock during dilation, and these values are used their values exist. Fortunately, each parameter has a unique effect
as the reference values in Eq. 1 during the initial elastic period. on the simulated results, and it is therefore relatively straightfor-
Below PR, the reservoir begins to recompact. We assume that ward to determine their values.
blocks that undergo large dilations will recompact more than blocks It has been estimated with our model that recompaction provides
with small dilations, and we defme the residual dilation fraction, about 60% of the drive energy in early cycles of Cold Lake CSS.!
fR, to be the fraction of the total dilation in a gridblock that is per- The value used for PR determines the pressure range over which
manent and unrecoverable. When pressure declines to PR, the this drive is obtained and hence controls the simulated production
minimum allowable porosity for each block is calculated from its pressures. Fig. 5 shows the simulated BHP during the first pro-
historical maximum dilation porosity andfR. The minimum porosi- duction cycle for different values of PRo Also shown are represen-
ty applies atp=O. With the coordinates of both ends of the recom- tative field data. A relatively low PR must be used to reproduce
paction function thus specified, Eq. 1 can be used to calculate the the observed pressure decline, which supports the need to model
recompaction compressibility for each block. Because a single fR
an elastic period before recompaction begins following dilation.
value is used for all blocks, the method allows more recompaction
Failure to match the production pressure decline will result in the
in blocks that dilated more, while guaranteeing that the porosity
simulator incorrectly predicting interwell communication events dur-
will never fall below the original value. The dashed lines in Fig.
ing multiwell simulations because the pressure difference between
3 show recompaction behavior for two gridblocks where less dila-
injecting and producing wells significantly affects this communi-
tion occurred than in the case represented by the solid lines. The
cation. Interwell communication can have a significant impact on
different recompaction compressibilities are readily apparent. (All
recompaction cases shown in Fig. 3 usedfR=0.50.) production performance 13 .!8 and must therefore be modeled ac-
The quadrilateral shape of the deformation model illustrated in curately. In addition, the pressure determines the amount of steam
Fig. 3 is consistent with the experimental observations of Roscoe flashing and solution-gas drivc< that occurs in the reservoir, and poor
et al. 17 Fig. 4 is schematically similar to the behavior their model pressure matches will therefore correspond to poor representations
would describe for a constant total stress condition in the Cold Lake of these important drive mechanisms.
Clearwater formation. The high initial porosity at deposition de- The residual dilation fraction,jR' largely controls the amount of
creased continually with burial until the point of maximum com- fluid the reservoir can produce. Fig. 6 shows how total fluid pro-
paction was reached during the last Ice Age, when an ice sheet 1800 duction varies with fR during a first-cycle simulation. Although
m thick covered the area. The current virgin reservoir conditions other parameters affect the total fluid production, fR has the larg-
were reached after melting of the ice reduced the effective stress est impact without substantially changing the produced WOR. The
acting on the sands. The path that is then followed during injection value used for fR is therefore chosen to provide a match of the ob-
and subsequent production is shown in Fig. 4. The similarity to served total fluid volume. We have found that anfR value in the
our model (Fig. 3) is apparent. range of 0.4 to 0.5 usually provides the best match for Cold Lake.

0.12~----------------------------------,

END Of
CYCLE 2 , ________________________ , 0.10 IMBIBITION ,,1
BOUNDING CURVE .. "
\ ,it
i" ........... ' .. /,'

mil
END Of ... ~-.- ...-.--.-.-.-. 0.08
,,
~ CYCLE ~\~// •
• 0.06
~l>.
,,
,,
~ ggf. -- m-____ -_m: m----;-= """ 0.04
/ DRAINAGE
0.02 START , .. BOUNDING CURVE

INITIAL RESERVOIR
CONDITIONS
CYCLE 2
CYCLE 3
0.1 0.2
" 0.3 0.4 0.5 0.6
"
0.00 ~--...-"--__,_-='--=rc=..::..--__--_,__---.---_,_--__I
0.7 0.8 0.9
PRESSURE WATER SATURATION

Fig. 8-Grldblock water relative permeability history over

Fig. 7-Grldblock porosity history over three cycles. three cycles.

202 SPE Reservoir Engineering. May 1991

Behavior in Subsequent Cycles. When the pressure begins to in-
crease following recompaction (e.g., during the next injection 0.7,-------------------,
,
cycle), the reservoir initially behaves elastically and the original 0.6 \ \ NO HYSTERESIS
compressibility is used. This continues until the dilation function ~ 0.5 ,,
VI
is intersected, at which time dilation again occurs. The subsequent
dilation and recompaction behavior is similar to that described previ-
0.4
------
............ -------WITH -HYSTERESIS
',
0.3 '"
ously. If the pressure begins to decline again before the dilation
line is intersected, then the elastic line is reversible and the grid- 0.451-====::-;:;;:;-;~~~;======:==~
t NO HYSTERESIS
block will move back along it until the recompaction line is inter- ~0.42
sected and recompaction continues. Fig. 7 shows the behavior of 51o 0.39 '" WITH HYSTERESIS .
a gridblock over three cycles. In Cycle 2, the block dilates more 0.. 0.36 " ,
than it did in the previous cycle. In Cycle 3, the block does ·not 0.33
reach the dilation envelope and returns along the elastic line until 0.30+----.----~----r---___.--__,r--__1
it intersects the recompaction line and recompaction resumes.
In summary, the deformation model consists of nonreversible di-
o 5 10 15 20 25 30
DISTANCE FROM FRACTURE (m)
lation and recompaction lines that are joined by reversible lines
describing elastic behavior. All these paths are history-dependent Fig. 9-Water saturation and porosity normal to a vertical
and behavior is therefore unique for each gridblock. fracture.

Permeability Changes. An increase in the absolute permeability

of the reservoir will accompany the porosity increase that results When the block pressure exceeds this value, dilation occurs. Large
from dilation. We model the permeability changes with values of KMUL are specified for fracture plane blocks so that the
permeability calculated by Eq. 2 reaches values much greater than
kxyz =kOxyz exp[KMULxyz (4) -4>0)/(1-4>0)1. ............. (2) in nonfracture plane blocks when the pressure is above fracture
where k and 4>.= permeability and porosity, respectively, ko and pressure.
4>0 = original permeability and porosity, respectively, and K MUL = Use of this fracture model requires the user to specify the orien-
user-defined multiplier. Different values of KMUL may be speci- tation, initiation pressure, and K MUL value for the fracture. The
fied in each gridblock. orientation (vertical or horizontal) is usually determined by stress
The relationship of permeability to porosity is not well defined, testing 27 or from field temperature profIles. 28 Fracture pressure
but Kozeny 19 predicts that the permeability of an unconsolidated can also be obtained by in-situ stress tests or, for horizontal frac-
sand is proportional to the cube of the porosity. This yields a per- tures, by integration of the density log. KMUL values are obtained
meability increase of 70% for a porosity increase of 20%. Dus- by history matching. A low KMUL value will keep injected fluid
seault and Rothenburg l5 predict a permeability increase of 50% close to the well, resulting in unreasonably high injection pressures
for a porosity increase of 25 % during shear failure. We found that and high produced WOR's. Larger KMUL values will cause steam
the simulation result is not particularly sensitive to changes on the to move farther away from the well, lowering both the injection
order of the difference between these two predictions. We use a pressures and the computed WOR's during production. We found
KMUL value that results in about a 70% increase in permeability that single-well simulation results are relatively insensitive to
at the highest porosity that our simulations reach. Permeabilities KMUL values as long as sufficiently large values are used. Better
decrease with porosity during recompaction, but a residual perme- estimates of KMUL values are obtained from multiwell models, in
ability increase exists as a result of the residual porosity increase which proper KMUL values are critical to matching observed in-
described earlier. terwell communication events. 18 We generally use K MUL values
that result in fracture permeabilities (at the injection pressure) that
Fracture Modeling are 100 to 500 times those of the reservoir permeability.
Steam injection into the Cold Lake reservoir is known to cause ten- Ito l6 used "super upstream convection weighting" in his simu-
sile fracturing of the formation. 5,8 Fracturing must be included in lations to move heat away from the well during injection and to
reservoir simulators to obtain a reasonable representation of reser- achieve a match of producing bottornhole temperatures (BHT's).
voir behavior. This is particularly important to match interwell com- However, he did not include a fracture in his model. We find that
munication observations from the field. 13,18 Published approaches the fracture is able to move heat away from the well efficiently,
to fracture modeling have ranged from simple manual methods to and our model matches BHT's without the use of any super weight-
attempts at fully rigorous representations. The most ambitious ap- ing factors.
proaches involve direct computation of fracture dimensions in ad- Calculating fracture permeability from Eqs. 1 and 2, rather than
dition to multiphase fluid flow and heat-transfer calculations. 20,21 as a step25 or bilinear 26 function of pressure, has two advantages.
Less advanced approaches incorporate fracture behavior with First, the fracture gridblocks undergo deformation phenomena simi-
sourcelsink methods 22-24 or apply an enhanced pressure-dependent lar to all other blocks and will have the permanent residual perme-
fracture permeability in a specified plane of gridblocks. 25,26 The ability increase that occurs because not all dilation is recoverable.
simplest approach is to impose an increased transmissibility dur- Second, Eq. 2 results in a smoothly increasing permeability with
ing injection on a predetermined group of blocks in the fracture no discontinuity at the fracture pressure, which significantly im-
plane, returning the transmissibility to its original value during pro- proves simulator stability. We obtained reductions of up to 40%
duction. 7,8,13 in CPU time by using Eq. 2 rather than a bilinear permeabilityIpres-
We have run multiwell simulations with as many as 4,500 grid- sure function to model the fracture plane permeability.
blocks and therefore do not wish to slow the performance of our
simulator by use of complex fracture modeling approaches. We Relative Permeability Hysteres.s
simulate fracturing by use of porosity-dependent transmissibility Many simulation studies of CSS in heavy-oil reservoirs result in
multipliers in a specified plane of gridblocks. This approach al- the observation that matching the field WOR requires reduction of
lows the fracture length to vary from cycle to cycle depending on the water relative permeability to well below experimentally meas-
leakoff conditions. The method also allows the fracture to stop grow- ured values, especially at low water saturations. 12,29,30 It then be-
ing if breakthrough to a higher-mobility region near a neighboring comes difficult, however, to obtain simulated injectivities that are
well occurs. as high as observed values. 12,29 It has been concluded from these
The fracture plane permeability is increased by use ofEqs. 1 and observations that heavy-oil reservoirs exhibit relative permeability
2 when the pressure exceeds the fracture pressure. PD is set equal hysteresis. 12 ,29,30 Hysteresis has been experimentally observed
to the desired fracture pressure for gridblocks in the fracture plane. (Ref. 30 provides a review of experimental results). Hysteresis in
SPE Reservoir Engineering, May 1991 203
both oil and water phases has been reported for a heavy-oil reser- relative permeabilities actually used therefore tend to be closer to
voir. 31 We observed relative permeability hysteresis experimen- the drainage bounding curve than the imbibition bounding curve
tally for Cold Lake cores and therefore include it in our simulations. throughout much of the simulation.
Several methods of modeling relative permeability hysteresis have
been reported. 12,29,30,32 Some methods switch manually from the Conclusions
imbibition to the drainage curve at the end of injection, 29,30 while
1. A reservoir deformation model has been developed that ap-
others make relative permeability a function of pressure, assum-
propriately represents the main features of oil-sand dilation and
ing that imbibition occurs only during injection and drainage only
recompaction occurring during CSS at Cold Lake. The model al-
during production. 12 However, neither approach is physically ac-
lows the simulator to match injection and production pressures and
curate. In addition, the curve-switching method can cause numeri-
flowback times that are otherwise difficult to reproduce.
cal instability.
2 .. The relative-permeability-hysteresis model presented provides
Because imbibition and drainage correspond to increasing and
decreasing water saturation, respectively (for a water-wet sand), a simple and effective means of modeling history-dependent be-
relative permeability hysteresis should be modeled as a function havior.
of saturation. Our approach is to define bounding imbibition and 3. The reservoir simulator is an important tool in developing an
drainage curves with common endpoints. When a saturation reversal understanding of the complex physical phenomena occurring dur-
occurs, the simulator calculates a scanning curve that defines how ing CSS in tar-sand formations.
the relative permeability will move toward the appropriate bound-
ing curve. The scanning curves must be determined by a method Nomenclature
that guarantees that calculated values always remain on or between A,B = constants in Eq. A-4
the bounding curves. Our approach satisfies this need; the Appen- c = compressibility, MPa- 1
dix provides mathematical details. We use hysteresis only for the fR = residual dilation fraction
water/oil two-phase relative permeabilities. Three-phase permea- k = absolute permeability, m 2
bilities are calculated using the gas/oil two-phase curves, the hys- krw = relative permeability to water
teretic water/oil curves, and Stone's first model. 33 K = ratio for scanning curve calculation (Eq. A-3)
Fig. 8 shows bounding curves derived from experimental meas-
KMUL = permeability multiplier (constant in Eq. 2)
urements on Cold Lake cores along with the water relative perme-
ability paths followed by a gridblock over several cycles in a CSS n = scanning-curve exponent
simulation. The gridblock water saturation initially increases from P = pressure, MPa
0.25 to 0.49, and the water relative permeability follows the imbi- PD = pressure at which dilation begins, MPa
bition bounding curve. At a saturation of 0.49, a saturation rever- PR = pressure at which recompaction begins, MPa
sal occurs and the relative permeability moves along a scanning Sw = water saturation, fraction
curve toward the drainage bounding curve. When the saturation c/> = porosity, fraction
reaches 0.34, a second reversal occurs and a new scanning curve
that approaches the imbibition curve is followed. Further scanning Subscripts
curves are followed after saturation reversals occur in later cycles. d = drainage
This example demonstrates the history-dependent nature of the rela- i = imbibition
tive permeabilities. Gridblocks with the same saturation can have ir = value at irreducible water saturation
significantly different oil or water permeabilities as a result of differ- max = maximum value
ent saturation histories.
P = value when saturation reversal occurs
The improved simulator match of produced WOR's made possi-
r = reference value
ble by the use of relative permeability hysteresis is well-
documented. 12 ,29-31 The use of hysteresis also affects the calcu- ro = value at residual oil saturation
lated fluid distributions and failed zone geometry. Fig. 9 shows x,y,z = flow direction
the effect of relative permeability hysteresis on calculated proftles o = value at original reservoir conditions
of water saturation and porosity normal to a vertical fracture fol-
lowing first-cycle steam injection. The fluid front and dilated zone Superscript
(as reflected by increased porosity) move farther into the reservoir *= normalized value
when hysteresis is included. This results in increased injectivity and
less extensive fractures during steam injection than when relative Acknowledgments
permeability hysteresis is not used. Because the advance of the fluid
The comments and suggestions of P.R. Kry, T.W. Miller, M.Y.
and porosity fronts determines the timing and magnitude of inter-
Kwad, S.Y. Bharatha, R.C.K. Wong, and E.S. Denbina during
well communication events, including relative permeability hyste-
resis is particularly important in multiwell simulations. the course of this work are appreciated. The permission of ERCL
The imbibition and drainage bounding curves and the scanning- to publish this paper is gratefully acknowledged.
curve exponents (see the Appendix) must be defined to use this
relative-permeability-hysteresis model. The imbibition curve is best References
obtained experimentally. The experimental curve is adjusted, if nec- 1. Denbina, E.S., Boberg, T.C., and Rotter, M.B.: "Evaluation of Key
essary, to honor any available data for the effective permeability Reservoir Drive Mechanisms in the Early Cycles of Steam Stimulation
to water at original reservoir conditions, which can be derived from at Cold Lake," SPERE (May 1991) 207-211.
analysis of observation-well pressure responses to nearby steam in- 2. Coats, K.H.: "A Highly Implicit Steamtlood Model," SPEl (Oct. 1978)
jection. 14 The bounding drainage curve also can be obtained ex- 369-83.
perimentally. History matching of the WOR's observed in the field, 3. Kendall, R.P. et aI.: "Development of a Multiple Application Reser-
however, will normally require adjustment of the experimental voir Simulator for Use on a Vector Computer," paper SPE 11483
presented at the 1983 SPE Middle East Oil Technical Conference and
drainage bounding curve.
Exhibition, Manama, Bahrain, March 14-17.
The scanning-curve exponents determine how rapidly the rela- 4. Oallant, R.I. and Dawson, A.O.: "Evolution of Technology for Com-
tive permeability approaches the bounding curve after a flow rever- mercial Bitumen Recovery at Cold Lake," paper 227 presented at the
sal. 'These values are extremely difficult to determine experimentally 1988 UNITARlUNDP Inti. Conference on Heavy Crude and Tar Sands,
and are normally used as history-match parameters in the simula- Edmonton, Aug. 7-12.
tion. We found it best to use scanning-curve exponents that move 5. Mainland, 0.0. and Lo, H.Y.: "Technological Basis for Commercial
the relative permeabilities rapidly toward the drainage bounding In-Situ Recovery of Cold Lake Bitumen," paper RTD3(l) presented
curve but more slowly toward the imbibition bounding curve. The at the 11th World Pet. Cong., London, Aug. 28-Sept. 2.

204 SPE Reservoir Engineering, May 1991

 6. de Souza, J.F.C.: "Commercial Success for Esso at Cold Lake," paper
presented at the 1986 AOSTRA Seminar on Advances in Petroleum
Recovery and Upgrading Technology, Calgary, June 12-13.
7. Dietrich, J .K.: "Cyclic Steaming of Tar Sands Through Hydraulically
Induced Fractures," SPERE (May 1986) 217-29. ::::E
8. Duerksen, J.H., Cruickshank, G.W., and Wasserman, M.L.: "Per- ...'"
0-
formance and Simulation of a Cold Lake Tar Sand Steam-Injection Pi-
lot," JPT(Oct. 1984) 1781-90. ...
...J

9. Settari, A. and Raisbeck, J.M.: "Analysis and Numerical Modeling ...'"'"

l-
(k~, )p.•.•..•.••. :
of Hydraulic Fracturing During Cyclic Steam Stimulation in Oil Sands," e(

JPT (Nov. 1981) 2201-12.

~ I INITIAL

10. Scott, J.D. and Kosar, K.M.: "Geotechnical Properties of Athabasca ...
0
N (k ~w ) ••••••••• •••••••••••••••••.
i POINT, P
i DRAINAGE
Oil Sands," paper 2-3 presented at the 1984 WestemResearchInst.IU.S.
:J
e(
P i BOUNDING CURVE
::::E :
DOE Tar Sand Symposium, Vail, Colorado, June 26-29. :
'"0z : ···········(k· d)
11. Agar, J.G.: "Behavior of Oil Sands at Elevated Temperatures," PhD i rw p
dissertation, U. of Alberta, Edmonton (1984). 0
12. Coats, K.H., Ramesh, A.B., and Winestock, A.G.: "Numerical Model- 0 (S~ )p
ing of Thermal Reservoir Behavior," The Oil Sands of Canada- NORMALIZED WATER SATURATION
Venezuela, 1977, Special Volume, CIM, Montreal (1977) 17, 399-410.
13. Pethrick, W.O., Sennhauser, E.S., and Harding, T.G.: "Numerical
Modeling Optimization of Cyclic Steam Stimulation in Cold Lake Oil
Sands," paper CIM 86-37-21 presented at the 1986 Petroleum Soc. Fig. A-1.,.-Normallzed water relative permeability curves.
of CIM Annual Technical Meeting, Calgary, June 8-11.
14. Settari, A., Kry, P.R., and Yee, C.-T.: "Coupling of Fluid Flow and
Soil Behavior To Model Injection Into Unconsolidated Oil Sands, .. paper
Appendix-Derivation of Equations for
CIM 88-39-72 presented at the 1988 Petroleum Soc. of CIM Annual the Relatlve.Permeablllty·Hysteresls
Technical Meeting, Calgary, June 12-16. Scanning Curve
15. Dusseault, M.B. and Rothenburg, L.: "Shear Dilatency and Permea- Fig. A-I shows a pair of normalized water relative permeability
bility Enhancement in Oil Sands," paper 32 presented at the 1988 UNI-
bounding curves, one for increasing water saturation (imbibition)
TARIUNDP Inti. Conference on Heavy Crude and Tar Sands,
Edmonton, Aug. 7-12. and one for decreasing water saturation (drainage). These curves
16. Ito, Y.: "The Introduction of the Microchanneling Phenomenon to Cy- are assumed to depend only on water saturation. Point P refers to
clic Steam Stimulation and Its Application to the Numerical Simulator conditions in a block when a saturation reversal occurs. The nor-
(Sand Deformation Concept)," SPE/ (Aug. 1984) 417-30. malized water relative permeabilities and water saturations in Fig.
17. Roscoe, K.H., Schofield, A.N., and Wroth, C.P.: "On Yielding of A-I are defmed as
Soils," Geotechnique (1958) 8, 22-53.
18. Vittoratos, E., Scott, G.R., and Beattie, C.I.: "Cold Lake Cyclic Steam k~ =krw/(krw)ro ................................ (A-I)
Stimulation: A Multiwell Process," SPERE (Feb. 1990) 19-24.
19. Scheidegger, A.E.: "The Physics of Flow Through Porous Media," and S;=[Sw-(Sw)ir]/[(Sw)ro -(Sw)ir]' ............... (A-2)
third edition, U. of Toronto Press, Toronto (1974) 103-04.
As the water saturation increases, the block's water relative per-
20. Farouq Ali, S.M. and Blunschi, J.: "Cyclic Steam Stimulation with
Formation Parting," paper CIM 83-34-45 presented at the 1983 Pe- meability will move toward the imbibition bounding curve (Fig.
troleum Soc. of CIM Annual Technical Meeting, Banff, May 10-13. A-I). We define a ratio, K, that compares the distance between
21. Settari, A. and Cleary, M.P.: "Three-Dimensional Simulation ofHy- the block's permeability and this bounding curve to the total dis-
draulic Fracturing," JPT (July 1984) 1177-90. tance between the bounding curves, all at the same water saturation:
22. Nghiem, L.X.: "Modeling Infinite-Conductivity Vertical Fractures
Using Source and Sink Terms," SPE/ (Aug. 1983) 633-44. K=[(k:'w; -k~)]/[(k~i -k!d)]' .................... (A-3)
23. Geshelin, B.M., Grabowski, J.W., and Pease, E.C.: "Numerical Study
Now, k~ can be calculated from Eq. A-3 if the bounding curves
of Transport of Injected and Reservoir Water in Fractured Reservoirs
During Steam Stimulation," paper SPE 10322 presented at the 1981 SPE and K are known. Assume a form for K [over the range 1 ~S;~
Annual Technical Conference and Exhibition, San Antonio, Oct. 5-7. (S;)p] of
24. Lin, C.Y.: "A New Approach for Simulation of Cyclic Steam Stimu-
lation Above Fracture Pressure by Modifying a Thermal Simulator," K =A +B(1-s;)n, .............................. (A-4)
paper SPE 18077 presented at the 1988 SPE Annual Technical Con- where A, B, and n are constants.
ference and Exhibition, Houston, Oct. 2-5.
There are three boundary conditions. First,
25. Soni, Y. and Harmon, R.A.: "Simulation of the Saner Ranch Fracture
Assisted Steamflood Pilot," J. Cdn. Pet. Tech. (Jan.-Feb. 1986) 25, K=O at S;=I, ................................. (A-5)
No. 1,57-70.
26. Coats, K.H.: "Thermal Simulation of Tar Recovery in Hydraulically which forces the scanning curve to approach the imbibition bound-
Fractured Formations, " paper presented at the 1983 Boeing Computer ing curve. Thus, A=O and n>O.
Services Spring Colloquium for the Geosciences, Seattle. The second boundary condition forces the scanning curve to pass
27. Gronseth, J.M. and Kry, P.R.: "Instantaneous Shut-In Pressure and
through the initial Point P:
Its Relationship to the Minimum In-Situ Stress," Hydraulic Fractur-
ing Stress Measurements, M.D. ZobackandB.C. Hainson (eds.), NatI. K=Kp at S;=(S;)p ............................. (A-6)
Academy Press, Washington, DC (1983) 55-60.
28. Vittoratos, E.S.: "Interpretation of Temperature ProfJIes From the Therefore, K=Kp{(I-S;)/[l-(S;)p]}n . ............. (A-7)
Steam-Stimulated Cold Lake Reservoir," paper SPE 15050 presented
at the 1986 SPE California Regional Meeting, Oakland, April 2-4. The third boundary condition requires that the derivative of K
29. Dietrich, J .K.: "Relative Permeability During Cyclic Steam Stimulation with respect to S; be finite in the range I ~S;~(S;)p. At S;=I,
of Heavy-Oil Reservoirs," JPT (Oct. 1981) 1987-89. this derivative is infinite when n < I. Therefore, n ~ 1.
30. Bang, H.W.: "Simulation Study Shows Hysteresis Effect on Oil Recov- So, by combining Eqs. A-3 and A-7 and defining Kp with Eq.
ery During a Cyclic Steam Process," Oil & Gas J. (Feb. 27, 1984) A-3, we obtain the expression for k~:
83-86.
31. Bennion, D.W., Moore, R.G., and Thomas, F.B.: "Effect of Rela-
tive Permeability on the Numerical Simulation of the Steam Stimula- k* =k* .- (k*) (k*) ][ 1 - S*
rwi p- rw P w
]n (k* .-k* )
tion Process, " paper CIM 83-34-46 presented at the 1983 Petroleum rw rwz [ * * * rwz rwd,
(krwi)P-(krwd)P 1-(Sw)p
Soc. of CIM Annual Technical Meeting, Banff, May 10-13.
32. Killough, J.E.: "Reservoir Simulation With History-Dependent Satu- n~ 1. ......................................... (A-8)
ration Functions," SPE/ (Feb. 1976) 37-48; Trans., AIME, 261.
33. Stone, H.L.: "Probability Model for Estimating Three-Phase Relative Because the other derivations are similar, only the fmal equa-
Permeability," JPT(Feb. 1970) 214-18; Trans., AIME, 249. tions are presented.

SPE Reservoir Engineering, May 1991 205

 Authors Normalized water relative permeability, normalized water satu-
ration decreasing.

k* =k * + [ (k~)p
-(kr'::.vt)p][ S! ]n(k*. -k * )
rw rwd * *
(krwi)P-(krwd)P *
(Sw)p TWI rwd.

n2: 1. ......................................... (A-9)

Normalized oil relative permeability, normalized water satura-
tion increasing.

* =k*. [(k:'c,)p-(k:'c,i)P][ I-S! ]n(k * -k*.)

k TO ro.+ * *
(krod)P-(kroi)P
*
1-(Sw)p
rod ro.'
Beattie Boberg McNab
n2: 1. ........................................ (A-I0)
Craig I. "aUla Is a reservoir engineering specialist at Esso
Resources Canada Ltd. In Calgary, where Is responsible for where the normalized oil relative permeability is
simulation studies of secondary and tertiary recovery proc-
esses. He previously worked on reservoir technology devel- lC';.o =kro/(kro)ir' ................................ (A-H)
opment and thermal simulation applications for the Cold Lake Normalized oil relative permeability, normalized water satura-
oll-sands project. He holds BASe and MASc degrees In chem·
lcal engineering from the U. of Toronto. Thom. . C. Boberg tion decreasing.
Is a research adviser and group leader for Cold Lake reser-
voir simulation with Exxon Production Research Co. In k* =k * _[( k*)
rod P- (k*)][
ro P s*]n (k * -k*.)
_w_
Houston. He Joined Exxon In 1959 and has spent 16 years TO rod * *
(krod)p-(kroi)P
*
(Sw)p
rod ro.,
working In thermal recovery, Including a year with Imperial
011 Ltd. working on Cold Lake. He previously was supervisor n2: 1. ........................................ (A-12)
of the Thermal Recovery Research Section and directed
studies of the application of steam stimulation to several large SI Metric Conversion Factors
Venezuelan reservoirs. He holds a as degree In physics from
William and Mary, BS and MS degrees In chemical engineer- bbl x 1.589 873 E-Ol m3
Ing from MIT, and a PhD In chemical engineering from the ep x 1.0* E-03 Pa's
U. of Michigan. Gordon S. McNab researches geostatlstlcs ft x 3.048* E-Ol m
applications at Esso Resources' Research Centre In Celgary. ft3 x 2.831 685 E-02 m3
He previously conducted research Into recovery processes OF (OF-32)/1.8 °C
for tar-sands reservoirs. He holds BASe and MASc degrees in. x 2.54* E+OO em
In chemical engineering from the U. of British Columbia, and psi x 6.894 757 kPa
E+OO
an ScD In engineering science from Oxford U. In the U.K. kPa- 1
psi -1 x 1.450 377 E-Ol
'Converslon factor Is exact. SPERE
Original SPE manuscript received for review March 23, 1989. Paper accepted for publica-
tion Oct. 31, 1990. Revised manuscript received July 13, 1990. Paper (SPE 18752) first
presented attha 1989 SPE California Ragional Meeting held In Bakersfield, April 5-7.

206 SPE Reservoir Engineering. May 1991