This preview shows page 1. Sign up to view the full content.
Unformatted text preview: INTEGRATION OF ROCK PHYSICS AND RESERVOIR
SIMULATION FOR THE INTERPRETATION OF TIMELAPSE
SEISMIC DATA AT WEYBURN FIELD, SASKATCHEWAN By
Leo Thomas Brown Copyright by Leo Thomas Brown 2002
All Rights Reserved ii A thesis submitted to the Faculty and the Board of Trustees of the Colorado School of
Mines in partial fulfillment of the requirements for the degree of Masters of Science
(Geophysics). Golden, Colorado
Date
Signed:
Leo Thomas Brown Approved:
Dr. Thomas L. Davis
Thesis Advisor Golden, Colorado
Date
Dr. Terence K. Young
Professor and Head,
Department of Geophysics iii ABSTRACT This thesis research integrates reservoir simulation with timelapse (4D) seismic
monitoring of reservoir processes, through rock and fluid physics modeling. During the CO2
injection program at Weyburn Field, changes in reservoir fluid pressure, fluid composition and
saturation are expected. Fluid models are developed for the acoustic properties of brine, oil, and
oil – CO2 mixtures. An anisotropic rock physics model is used to calculate the sensitivity of the
seismic properties of the reservoir to fluid and stress changes. Reservoir simulation of the
enhanced oil recovery operations provides estimates of the changes in pore pressure, saturation,
and fluid composition. The reservoir simulation output is combined with the rock and fluid
physics models to estimate the change in seismic properties of the reservoir. These predicted
changes are compared to the timelapse difference anomalies in the Pwave seismic data.
The fluid physics models are based on existing empirical relations, laboratory
measurements and equation of state modeling. The pressure and porositydependent anisotropic
model for the reservoir zones is built from ultrasonic measurements on core samples, analysis of
geophysical logs, and effective medium modeling for saturated, fractured rocks. This model can
be used to calculate the density and elastic stiffness matrix of a transversely isotropic rock with a
horizontal (HTI) symmetry axis. Synthetic seismic modeling shows that changes in the upper and
lower (Marly and Vuggy) reservoir zones are not independently resolved in the seismic data. The
magnitude of the expected changes in Pwave reflection amplitude due to CO2 injection is 15% to
20%, and should be detected in the timelapse seismic data.
Through interpretation of Pwave seismic data volumes, areas effectively contacted by
CO2 are identified. The observed timelapse anomalies correlate strongly with the modeled CO2
movement and Pimpedance decrease. The differences in the seismic data include more spread
out anomalies, differences in location of anomalies, and evidence for CO2 fingering along
fractures. This thesis research provides the forward model for calculating changes in seismic
properties from reservoir processes. It can be used in future research in integrated reservoir
inversion to refine the reservoir model and the reservoir simulation process. iv TABLE OF CONTENTS ABSTRACT .......................................................................................................................iv
TABLE OF CONTENTS .................................................................................................... v
LIST OF FIGURES............................................................................................................ix
LIST OF TABLES .........................................................................................................xviii
ACKNOWLEDGMENTS................................................................................................. xx
CHAPTER 1 INTRODUCTION ........................................................................................ 1
1.1 Introduction............................................................................................................................1
1.2 Weyburn Field .......................................................................................................................2
1.3 Weyburn Field Geology and Reservoir Properties ................................................................3
1.4 Enhanced Oil Recovery Operations.......................................................................................5
1.5 RCP Phases VIII and IX Study..............................................................................................6 CHAPTER 2 ACOUSTIC PROPERTIES OF RESERVOIR FLUIDS.............................. 9
2.1 Introduction............................................................................................................................9
2.2 Weyburn Reservoir Fluids ...................................................................................................10
2.2.1 Brine........................................................................................................................................... 10
2.2.2 Oil............................................................................................................................................... 12
2.2.3 CO2 ............................................................................................................................................. 14
2.2.4 Oil + CO2 Mixtures .................................................................................................................... 16 2.3 Summary..............................................................................................................................19 CHAPTER 3 ULTRASONIC MEASUREMENTS ON CORE SAMPLES AND DRY
ROCK ACOUSTIC PROPERTIES .................................................................................. 21
3.1 Introduction..........................................................................................................................21
3.2 General Fluid Substitution Theory.......................................................................................22
3.3 Ultrasonic Testing Procedure...............................................................................................23
3.3.1 Sample Preparation..................................................................................................................... 23
3.3.2 Ultrasonic Testing Equipment.................................................................................................... 24
3.3.3 Ultrasonic Measurements ........................................................................................................... 26 v 3.4 Dry Rock Ultrasonic Measurements....................................................................................27
3.4.1 Marly Unit .................................................................................................................................. 27
3.4.2 Vuggy Unit................................................................................................................................. 31 3.5 Saturated Rock Ultrasonic Measurements ...........................................................................34
3.5.1 Constant Differential Pressure.................................................................................................... 35
3.5.2 Constant Confining Pressure ...................................................................................................... 37
3.5.3 Seismic Attributes ...................................................................................................................... 39 3.6 Summary of PanCanadian Rock Testing Results ................................................................42
3.7 Summary..............................................................................................................................43 CHAPTER 4 ISOTROPIC ROCK PHYSICS MODELING ............................................ 45
4.1 Introduction..........................................................................................................................45
4.2 Gassmann’s Theory .............................................................................................................45
4.2.1 Comparison of Gassmann’s Equation and Lab Results.............................................................. 46
4.2.2 Modeling with Gassmann’s Equation and Lab Data .................................................................. 48 4.3 Fluid Substitution on Geophysical Logs..............................................................................52
4.4 Integration of Ultrasonic Core and Geophysical Log Data..................................................63
4.4.1 Pressure Dependence of Elastic Moduli..................................................................................... 64
4.4.2 Porosity Dependence of Elastic Moduli ..................................................................................... 65 4.5 Summary..............................................................................................................................69 CHAPTER 5 ANISOTROPIC ROCK PHYSICS MODELING ...................................... 71
5.1 Introduction..........................................................................................................................71
5.2 Fractures ..............................................................................................................................71
5.2.1 Previous Fracture Characterization Studies................................................................................ 72
5.2.2 Hudson’s Model for Cracked Media .......................................................................................... 73
5.2.3 Single Fracture Set – Hudson’s Model....................................................................................... 74
5.2.4 Multiple Fracture Sets ................................................................................................................ 78 5.3 Stiffness Tensor of Saturated Anisotropic Rock – Brown and Korringa’s Equation...........79
5.4 PressureDependent, Anisotropic, Saturated Rock Physics Model .....................................80
5.5 Summary..............................................................................................................................88 CHAPTER 6 SYNTHETIC SEISMIC MODELING ....................................................... 90
6.1 Introduction..........................................................................................................................90 vi 6.2 Seismic Impedance Modeling..............................................................................................90
6.3 Synthetic Seismograms........................................................................................................95
6.3.1 Correlation of Reflectors and Reflections .................................................................................. 95
6.3.2 Sensitivity of Synthetics to Fluid Saturation .............................................................................. 98 6.4 Summary............................................................................................................................101 CHAPTER 7 ACOUSTIC PROPERTIES OF WEYBURN OIL – CO2 SYSTEM ....... 102
7.1 Introduction........................................................................................................................102
7.2 Compositional Reservoir Simulation.................................................................................103
7.3 Fluid Modeling Options.....................................................................................................104
7.3.1 FLAG 4.0 ................................................................................................................................. 104
7.3.2 STRAPP ................................................................................................................................... 104
7.3.3 Reservoir Simulator.................................................................................................................. 108
7.3.4 Cubic EOS built on STRAPP Database ................................................................................... 108 7.4 Summary............................................................................................................................114 CHAPTER 8 RESERVOIR SIMULATION AND SEISMIC ATTRIBUTE MODELING
......................................................................................................................................... 115
8.1 Introduction........................................................................................................................115
8.2 Reservoir Simulation .........................................................................................................116
8.2.1 Summary of Reservoir Model .................................................................................................. 116
8.2.2 Summary of Output Data ......................................................................................................... 122 8.3 Mapping Reservoir Simulation Data to Seismic Attributes...............................................132
8.4 Modeling Results ...............................................................................................................134
8.5 Discussion of Modeling .....................................................................................................146
8.6 Summary............................................................................................................................149 CHAPTER 9 INTEGRATED SEISMIC DATA INTERPRETATION ......................... 150
9.1 Introduction........................................................................................................................150
9.2 Seismic Reflection Surveys ...............................................................................................150
9.3 PWave Amplitude Interpretation......................................................................................153
9.4 Integration of Seismic and Reservoir Simulation Data......................................................160
9.5 Recommendations for Further Research............................................................................166
9.6 Summary............................................................................................................................167 vii REFERENCES CITED ................................................................................................... 170
APPENDIX: CUBIC EQUATION OF STATE MODEL FOR DENSITY AND BULK
MODULUS OF GAS, CONDENSATE, LIGHT AND MEDIUM OIL ........................ 173
A.1 Introduction......................................................................................................................173
A.2 Procedures........................................................................................................................175
A.3 Source of Data..................................................................................................................178
A.4 Results..............................................................................................................................188
A.4.1 Cubic EOS Coefficients ......................................................................................................... 188
A.4.2 Sample Code and Calibration Data ........................................................................................ 188
A.4.3 Example Results for Cubic EOS ............................................................................................ 190 A.5 Discussion........................................................................................................................197
A.5.1 Consideration of Phase Behavior ........................................................................................... 202
A.5.2 Relationship between Descriptive Properties of Oil and Molecular Weight .......................... 206 A.6 Summary ..........................................................................................................................207 viii LIST OF FIGURES Figure 1.1 Location of Weyburn Field in Saskatchewan, Canada. ................................................................ 2
Figure 1.2 Stratigraphic column for Weyburn Field. Left side is after Dietrich and Magnusson (1998).
Right side is after Wegelin (1984). ....................................................................................................... 4
Figure 1.3 Schematic cross section of reservoir showing EOR strategy........................................................ 6
Figure 1.4 Location of RCP 4D 9C seismic survey area. Modified from PanCanadian Petroleum. .......... 7
Figure 1.5 EOR infrastructure at RCP 4D 9C seismic survey area............................................................. 7
Figure 2.1 Bulk modulus of brine as a function of salinity and pressure, at reservoir temperature. Modeled
with FLAG 4.0. ................................................................................................................................... 11
Figure 2.2 Density of brine as a function of salinity and pressure, at reservoir temperature. Modeled with
FLAG 4.0. ........................................................................................................................................... 12
Figure 2.3 Bulk modulus of Weyburn crude oil as a function of composition and pressure, at reservoir
temperature. Modeled with FLAG 4.0. .............................................................................................. 13
Figure 2.4 Density of Weyburn crude oil as a function of composition and pressure, at reservoir
temperature. Modeled with FLAG 4.0. .............................................................................................. 14
Figure 2.5 Bulk modulus of CO2 as a function of pressure, at reservoir temperature. Modeled with
STRAPP. ............................................................................................................................................. 15
Figure 2.6 Density of CO2 as a function of pressure, at reservoir temperature. Modeled with STRAPP... 15
Figure 2.7 Comparison of ultrasonic measurement results of Pwave velocity with FLAG model............. 17
Figure 2.8 Bulk Modulus of oil and CO2 mixtures as a function of composition and pressure, at reservoir
temperature. Modeled with STRAPP and FLAG. .............................................................................. 17
Figure 2.9 Density of oil and CO2 mixtures as a function of composition and pressure, at reservoir
temperature. Modeled with STRAPP and FLAG. .............................................................................. 18
Figure 2.10 Bulk modulus of reservoir fluids as functions of composition and pressure, at reservoir
temperature. Modeled with STRAPP and FLAG. .............................................................................. 20
Figure 2.11 Density of reservoir fluids as functions of composition and pressure, at reservoir temperature.
Modeled with STRAPP and FLAG..................................................................................................... 20
Figure 3.1 Example of change in Pwave velocity of reservoir rock under waterflood due to changes in
fluid saturation and pressure................................................................................................................ 23
Figure 3.2 Schematic of ultrasonic testing setup for saturated cores. .......................................................... 24 ix Figure 3.3 Recording setup for ultrasonic testing ........................................................................................ 25
Figure 3.4 Variation of ultrasonic P and Swave velocities with pressure for dry Marly sample PC 1213
4626. Hysteresis can be seen in the variation between the increasing versus decreasingpressure
measurements. ..................................................................................................................................... 28
Figure 3.5 Variation of ultrasonic P and Swave velocities with pressure for dry Marly sample PC 1411
4659.5.................................................................................................................................................. 29
Figure 3.6 Variation of dry rock bulk modulus with pressure for Marly samples PC 1213 4626 and PC 1411 4659.5............................................................................................................................................. 29
Figure 3.7 Variation of Lame’s parameter and shear modulus with pressure for dry Marly sample PC 1213
4626..................................................................................................................................................... 30
Figure 3.8 Variation of Lame’s parameter and shear modulus with pressure for dry Marly sample PC 1411
4659.5.................................................................................................................................................. 30
Figure 3.9 Variation of ultrasonic P and Swave velocities with pressure for dry Vuggy sample PC 1411
4698..................................................................................................................................................... 32
Figure 3.10 Variation of ultrasonic P and Swave velocities with pressure for dry Vuggy sample PC 1213
4643..................................................................................................................................................... 32
Figure 3.11 Variation of dry rock bulk modulus with pressure for Vuggy samples PC 1213 4643 and PC
1411 4698. ......................................................................................................................................... 33
Figure 3.12 Variation of Lame’s parameter and shear modulus with pressure for dry Vuggy sample PC 1411 4698................................................................................................................................................ 33
Figure 3.13 Variation of Lame’s parameter and shear modulus with pressure for dry Vuggy sample PC 1213 4643................................................................................................................................................ 34
Figure 3.14 Variation of ultrasonic Pwave velocity with pore pressure at constant differential pressure for
saturated Vuggy sample PC 1411 4698. ............................................................................................ 36
Figure 3.15 Variation of ultrasonic Swave velocity with pore pressure at constant differential pressure for
saturated Vuggy sample PC 1411 4698. ............................................................................................ 36
Figure 3.16 Variation of ultrasonic Pwave velocity with pore pressure at constant confining pressure for
saturated Vuggy sample PC 1411 4698. ............................................................................................ 38
Figure 3.17 Variation of ultrasonic Swave velocity with pore pressure at constant confining pressure for
saturated Vuggy sample PC 1411 4698. ............................................................................................ 38
Figure 3.18 Variation of ultrasonic VP/VS with pore pressure at constant confining pressure for saturated
Vuggy sample PC 1411 4698............................................................................................................. 40 x Figure 3.19 Variation of shear splitting parameter with pore pressure at constant confining pressure for
saturated Vuggy sample PC 1411 4698. ............................................................................................ 41
Figure 3.20 Variation of shear splitting parameter with differential pressure for dry Marly samples PC 1213 4626 and 1411 4659.5................................................................................................................... 41
Figure 4.1 Comparison of Pwave velocities calculated with Gassmann’s equation and ultrasonic core
measurements. Confining pressure was fixed at 34.5 MPa. With increasing pore pressure,
differential pressure decreases............................................................................................................. 47
Figure 4.2 Variation of Marly dolostone Pwave velocity with fluid composition and pressure, calculated
with Gassmann’s equation................................................................................................................... 50
Figure 4.3 Variation of Marly dolostone VP/VS with fluid composition and pressure, calculated with
Gassmann’s equation........................................................................................................................... 50
Figure 4.4 Variation of Vuggy limestone Pwave velocity with fluid composition and pressure, calculated
with Gassmann’s equation................................................................................................................... 51
Figure 4.5 Variation of Vuggy limestone VP/VS with fluid composition and pressure, calculated with
Gassmann’s equation........................................................................................................................... 51
Figure 4.6 Location of wells 02042364 and 41041663 in relation to RCP study area. The locations of core
samples are also shown. ...................................................................................................................... 53
Figure 4.7 Geophysical logs from well 02042364. Marly zone is from approximately 1398 m to 1405 m,
Vuggy zone is from approximately 1405 to 1423 m. .......................................................................... 54
Figure 4.8 Dry rock bulk modulus versus depth for well 02042364. Marly zone is approximately 1398 m
to 1405 m; Vuggy zone is approximately 1405 m to 1423 m. ............................................................ 56
Figure 4.9 Dry rock bulk modulus versus porosity for Vuggy Zone, well 02042364. ................................ 57
Figure 4.10 VP versus depth, as a function of fluid saturation for well 02042364. Marly zone is 1398 m to
1405 m; Vuggy zone is 1405 to 1423 m.............................................................................................. 59
Figure 4.11 Percent change in VP from brine saturation, for well 02042364. Marly zone is 1398 m to 1405
m; Vuggy zone is 1405 to 1423 m....................................................................................................... 59
Figure 4.12 VP/VS versus depth, as a function of fluid saturation, for well 02042364. Marly zone is 1398
m to 1405 m; Vuggy zone is 1405 to 1423 m. .................................................................................... 60
Figure 4.13 Percent change in VP/VS from brine saturation, for well 02042364. Marly zone is 1398 m to
1405 m; Vuggy zone is 1405 to 1423 m.............................................................................................. 60
Figure 4.14 Density versus depth, as a function of fluid saturation, for well 02042364. Marly zone is 1398
m to 1405 m; Vuggy zone is 1405 to 1423 m. .................................................................................... 61 xi Figure 4.15 Percent change in density from brine saturation, for well 02042364. Marly zone is 1398 m to
1405 m; Vuggy zone is 1405 to 1423 m.............................................................................................. 61
Figure 4.16 Blocked Pwave impedance versus depth, as a function of fluid saturation, for well 02042364.
............................................................................................................................................................. 62
Figure 4.17 Blocked Swave impedance versus depth, as a function of fluid saturation, for well 02042364.
............................................................................................................................................................. 62
Figure 4.18 Variation of 1/KDRY with porosity for the Marly zone. ............................................................ 67
Figure 4.19 Variation of 1/µDRY with porosity for the Marly zone. ............................................................. 67
Figure 4.20 Variation of 1/KDRY with porosity for the Vuggy zone. ........................................................... 68
Figure 4.21 Variation of 1/µDRY with porosity for the Vuggy zone............................................................. 68
Figure 5.1 Variation of c33 (Pwave modulus) with crack density at constant differential pressure,
calculated using Hudson’s crack model. ............................................................................................. 76
Figure 5.2 Variation of c44 (S1wave modulus) and c66 (S2wave modulus) with crack density at constant
differential pressure, calculated using Hudson’s crack model. ........................................................... 76
Figure 5.3 Variation of c33 (Pwave modulus) with differential pressure at constant crack density,
calculated using Hudson’s crack model. ............................................................................................. 77
Figure 5.4 Variation of c44 (S1wave modulus) and c66 (S2wave modulus) with differential pressure at
constant crack density, calculated using Hudson’s crack model. ........................................................ 78
Figure 5.5 Modeled variation of c33 (Pwave modulus) with fluid bulk modulus at constant crack density
and differential pressure. ..................................................................................................................... 80
Figure 5.6 Marly bulk density as a function of fluid saturation and pore pressure. Fluid density values are
shown in Figure 2.11........................................................................................................................... 81
Figure 5.7 Vuggy bulk density as a function of fluid saturation and pore pressure. Fluid density values are
shown in Figure 2.11........................................................................................................................... 82
Figure 5.8 Marly Pwave velocity as a function of fluid saturation and pore pressure at constant crack
density and confining pressure. ........................................................................................................... 82
Figure 5.9 Marly S1wave velocity as a function of fluid saturation and pore pressure at constant crack
density and confining pressure. ........................................................................................................... 83
Figure 5.10 Marly S2wave velocity as a function of fluid saturation and pore pressure at constant crack
density and confining pressure. ........................................................................................................... 83
Figure 5.11 Vuggy Pwave velocity as a function of fluid saturation and pore pressure at constant crack
density and confining pressure. ........................................................................................................... 84 xii Figure 5.12 Vuggy S1wave velocity as a function of fluid saturation and pore pressure at constant crack
density and confining pressure. ........................................................................................................... 84
Figure 5.13 Vuggy S2wave velocity as a function of fluid saturation and pore pressure at constant crack
density and confining pressure. ........................................................................................................... 85
Figure 5.14 Marly VP/VS1 as a function of fluid saturation and pore pressure at constant crack density and
confining pressure. .............................................................................................................................. 86
Figure 5.15 Vuggy VP/VS1 as a function of fluid saturation and pore pressure at constant crack density and
confining pressure. .............................................................................................................................. 86
Figure 6.1 Variation of Marly Pwave impedance with fluid saturation and pressure at constant crack
density. ................................................................................................................................................ 91
Figure 6.2 Variation of Marly S1wave impedance with fluid saturation and pressure at constant crack
density. ................................................................................................................................................ 91
Figure 6.3 Variation of Vuggy Pwave impedance with fluid saturation and pressure at constant crack
density. ................................................................................................................................................ 92
Figure 6.4 Variation of Vuggy S1wave impedance with fluid saturation and pressure at constant crack
density. ................................................................................................................................................ 92
Figure 6.5 Variation of Pwave impedance with depth, Weyburn reservoir model. Crack density of Marly
= 0.03, Vuggy = 0.1, PC = 22 MPa, PP = 15 MPa................................................................................ 94
Figure 6.6 Variation of S1wave impedance with depth, Weyburn reservoir model. Crack density of Marly
= 0.03, Vuggy = 0.1, PC = 22 MPa, PP = 15 MPa................................................................................ 94
Figure 6.7 Pwave wavelet for RCP seismic survey (from Ida Herawati). .................................................. 96
Figure 6.8 Pwave synthetic seismogram from wells 02042364 and 410416613, wavelet in Figure 6.7. ... 97
Figure 6.9 Change in Pwave synthetic trace due to fluid change in both Marly and Vuggy zones............ 99
Figure 6.10 Change in Pwave synthetic trace due to fluid change in Marly zone. ..................................... 99
Figure 6.11 Change in Pwave synthetic trace due to fluid change in Vuggy zone................................... 100
Figure 7.1 Variation of effective bulk modulus with composition and pressure for Weyburn oil – CO2
system. Temperature is 63 °C........................................................................................................... 106
Figure 7.2 Variation of average density with composition and pressure for Weyburn oil – CO2 system.
Temperature is 63 °C......................................................................................................................... 106
Figure 7.3 Variation of singlephase bulk modulus with composition and pressure for Weyburn oil – CO2
system. Temperature is 63° C. (Based on STRAPP modeling)........................................................ 109
Figure 7.4 Variation of bulk modulus with average molecular weight for Weyburn oil – CO2 system.
Temperature is 63 °C, pressure varies from 3 to 27 MPa.................................................................. 110 xiii Figure 7.5 Variation of bulk modulus with pressure for Weyburn oil – CO2 system. Temperature is 63 °C,
average molecular weight varies from 44 to 205. ............................................................................. 110
Figure 7.6 Range of average molecular weights and pressures that can exist for single fluid phases in the
Weyburn oilCO2 system. Calculated with STRAPP. ...................................................................... 111
Figure 7.7 Comparison of predicted bulk modulus from cubic EOS with bulk modulus calculated in
STRAPP for Weyburn oil – CO2 system........................................................................................... 113
Figure 8.1 Marly zone thickness in RCP study area, from reservoir simulation model............................. 118
Figure 8.2 Vuggy zone thickness in RCP study area, from reservoir simulation model............................ 118
Figure 8.3 Marly zone porosity, from reservoir simulation model. ........................................................... 119
Figure 8.4 Vuggy zone porosity, from reservoir simulation model. .......................................................... 119
Figure 8.5 Marly zone xdirection permeability, from reservoir simulation model................................... 120
Figure 8.6 Vuggy zone xdirection permeability, from reservoir simulation model.................................. 120
Figure 8.7 Marly zone ydirection permeability, from reservoir simulation model................................... 121
Figure 8.8 Vuggy zone ydirection permeability, from reservoir simulation model.................................. 121
Figure 8.9 Reservoir simulation area with injection and production wells shown. ................................... 123
Figure 8.10 Marly zone pore pressure at time of repeat survey, from reservoir simulation. ..................... 124
Figure 8.11 Vuggy zone pore pressure at time of repeat survey, from reservoir simulation. .................... 124
Figure 8.12 Marly zone change in pore pressure, from reservoir simulation model.................................. 125
Figure 8.13 Vuggy zone change in pore pressure, from reservoir simulation model. ............................... 125
Figure 8.14 Marly zone change in fractional oil saturation, from reservoir simulation............................. 126
Figure 8.15 Vuggy zone change in fractional oil saturation, from reservoir simulation............................ 126
Figure 8.16 Marly zone change in fractional water saturation, from reservoir simulation. ....................... 127
Figure 8.17 Vuggy zone change in fractional water saturation, from reservoir simulation. ...................... 127
Figure 8.18 Marly zone repeat survey fractional “vapor” saturation, from reservoir simulation. ............. 129
Figure 8.19 Vuggy zone repeat survey fractional “vapor” saturation, from reservoir simulation. ............ 129
Figure 8.20 Marly zone change in CO2 mole fraction in hydrocarbons, from reservoir simulation. ......... 130
Figure 8.21 Vuggy zone change in CO2 mole fraction in hydrocarbons, from reservoir simulation. ......... 130
Figure 8.22 Marly zone change in fluid density, from reservoir simulation.............................................. 131
Figure 8.23 Vuggy zone change in fluid density, from reservoir simulation............................................. 131
Figure 8.24 Process for mapping reservoir simulation data to predicted seismic difference anomaly. ..... 133
Figure 8.25 Marly zone change in fluid bulk modulus, from reservoir simulation.................................... 135
Figure 8.26 Vuggy zone change in fluid bulk modulus, from reservoir simulation. ................................. 135
Figure 8.27 Marly zone change in Pwave velocity, from reservoir simulation. ....................................... 136 xiv Figure 8.28 Vuggy zone change in Pwave velocity, from reservoir simulation. ...................................... 136
Figure 8.29 Marly zone change in Pwave impedance, from reservoir simulation.................................... 137
Figure 8.30 Vuggy zone change in Pwave impedance, from reservoir simulation................................... 137
Figure 8.31 Marly zone change in S1wave velocity, from reservoir simulation....................................... 138
Figure 8.32 Vuggy zone change in S1wave velocity, from reservoir simulation...................................... 138
Figure 8.33 Marly zone change in S1wave impedance, from reservoir simulation. ................................. 139
Figure 8.34 Vuggy zone change in S1wave impedance, from reservoir simulation. ................................ 139
Figure 8.35 Marly zone change in VP/VS1, from reservoir simulation....................................................... 140
Figure 8.36 Vuggy zone change in VP/VS1, from reservoir simulation...................................................... 140
Figure 8.37 Weyburn reservoir change in pore pressure from baseline to repeat survey, based on reservoir
simulation.......................................................................................................................................... 142
Figure 8.38 Weyburn reservoir change in CO2 content, as approximate fraction of pore volume, from
baseline to repeat survey, based on reservoir simulation. ................................................................. 143
Figure 8.39 Weyburn reservoir percent change in Pimpedance from baseline to repeat survey, based on
reservoir simulation and rock physics modeling. .............................................................................. 144
Figure 8.40 Weyburn reservoir percent change in S1impedance from baseline to repeat survey, based on
reservoir simulation and rock physics modeling. .............................................................................. 145
Figure 8.41 Weyburn reservoir percent change in VP/VS1, from baseline to repeat survey, based on
reservoir simulation and rock physics modeling. .............................................................................. 146
Figure 9.1 Baseline Pwave RMS amplitude for +/ 2 ms window around Marly trough reflection. ........ 154
Figure 9.2 Difference (Repeat – Baseline) in Pwave RMS amplitudes, window of +/2 msec around Marly
trough reflection. ............................................................................................................................... 155
Figure 9.3 Percent difference (Repeat – Baseline) in Pwave RMS amplitudes, window of +/2 msec
around Marly trough reflection. ........................................................................................................ 156
Figure 9.4 Difference (Repeat – Baseline) in Pwave RMS amplitudes, window of 3 msec above and 9
msec below Marly trough reflection.................................................................................................. 157
Figure 9.5 Percent difference (Repeat – Baseline) in Pwave RMS amplitudes, window of 3 msec above
and 9 msec below Marly trough reflection........................................................................................ 158
Figure 9.6 Process for interpretation of timelapse seismic data, using an integrated reservoir simulation
and rock and fluid physics model...................................................................................................... 161
Figure 9.7 Weyburn reservoir percent change in Pimpedance, predicted from reservoir simulation and
rock physics modeling. Horizontal CO2 injection wells are shown in red, other wells in light yellow.
........................................................................................................................................................... 162 xv Figure 9.8 Difference (Repeat – Baseline) in Pwave RMS amplitudes, window of +/2 msec around Marly
trough reflection, with horizontal wells superimposed. CO2 injectors are shown in red, other wells
are shown in light yellow. Areas with CO2 fingering are circled in black. Wells of special interest
are denoted by number. ..................................................................................................................... 164
Figure 9.9 Cumulative CO2 injection in horizontal wells as of 9/2001. Area of circle is proportional to
volume of injected CO2. CO2 injection wells are shown in bold...................................................... 165
Figure 9.10 Possible future process flow for reservoir model inversion from timelapse seismic surveys,
reservoir simulation, and rock and fluid physics models. ................................................................. 166
Figure A.1 Generalized phase diagram for miscible hydrocarbon mixtures.............................................. 174
Figure A.2 Generalized phase diagram for miscible hydrocarbon mixtures.............................................. 175
Figure A.3 Comparison of ultrasonic velocity test results with estimates from STRAPP......................... 180
Figure A.4 Temperature and pressure range for data modeled in STRAPP. Bounds are defined by
combination of maximum and minimum geothermal and geobaric gradients (Reginald Beardsley,
Unocal, personal communication, 2002)........................................................................................... 182
Figure A.5 Density range as a function of average molecular weight. Temperature and pressure vary. .. 185
Figure A.6 Density range as a function of pore pressure. Temperature and molecular weight vary. ....... 185
Figure A.7 Density range as a function of temperature. Molecular weight and pressure vary. ................ 186
Figure A.8 Bulk modulus range as a function of average molecular weight. Temperature and pressure
vary.................................................................................................................................................... 186
Figure A.9 Bulk modulus range as a function of pore pressure. Molecular weight and temperature vary.
........................................................................................................................................................... 187
Figure A.10 Bulk modulus range as a function of temperature. Molecular weight and pressure vary. .... 187
Figure A.11 Density results from cubic EOS for mw = 25........................................................................ 191
Figure A.12 Density results from cubic EOS for mw = 80........................................................................ 191
Figure A.13 Density results from cubic EOS for mw = 150...................................................................... 192
Figure A.14 Density results from cubic EOS for mw = 230...................................................................... 192
Figure A.15 Bulk modulus results from cubic EOS for mw = 25.............................................................. 193
Figure A.16 Bulk modulus results from cubic EOS for mw = 80.............................................................. 193
Figure A.17 Bulk modulus results from cubic EOS for mw = 150............................................................ 194
Figure A.18 Bulk modulus results from cubic EOS for mw = 230............................................................ 194
Figure A.19 Density results from cubic EOS for pressure = 15 MPa........................................................ 195
Figure A.20 Density results from cubic EOS for pressure = 60 MPa........................................................ 195
Figure A.21 Bulk Modulus results from cubic EOS for pressure = 15 MPa. ............................................ 196 xvi Figure A.22 Bulk Modulus results from cubic EOS for pressure = 60 MPa. ............................................ 196
Figure A.23 Density residuals versus density............................................................................................ 197
Figure A.24 Density residuals versus average molecular weight. ............................................................. 198
Figure A.25 Density residuals versus pressure. ......................................................................................... 198
Figure A.26 Density residuals versus temperature. ................................................................................... 199
Figure A.27 Bulk modulus residuals versus bulk modulus........................................................................ 199
Figure A.28 Bulk modulus residuals versus average molecular weight. ................................................... 200
Figure A.29 Bulk modulus residuals versus pressure. ............................................................................... 200
Figure A.30 Bulk modulus residuals versus temperature. ......................................................................... 201
Figure A.31 Critical temperature versus average molecular weight for hydrocarbon samples. ................ 204
Figure A.32 Critical pressure versus average molecular weight for hydrocarbon samples. ...................... 204
Figure A.33 Density versus pressure relationships for two pseudosamples with the same average
molecular weight. Data modeled in STRAPP. ................................................................................. 205
Figure A.34 Bulk modulus versus pressure relationships for two pseudosamples with the same average
molecular weight. Data modeled in STRAPP. ................................................................................. 205
Figure A.35 Relationship between density and molecular weight............................................................. 207 xvii LIST OF TABLES Table 3.1 Description of Marly unit samples............................................................................................... 28
Table 3.2 Description of Vuggy unit samples.............................................................................................. 31
Table 4.1 Comparison of sensitivity of seismic attributes to changes in fluid composition and pressure for
Marly and Vuggy zones. Based on isotropic rock physics model from ultrasonic core testing. ........ 52
Table 4.2 Comparison of average values of KDRY, µDRY from borehole sonic log and ultrasonic core
measurements. ..................................................................................................................................... 57
Table 4.3 Comparison of sensitivity of seismic attributes to changes in fluid composition and pressure for
Marly and Vuggy zones. Based on geophysical logs for well 02042364........................................... 63
Table 4.4 Comparison of sensitivity of seismic attributes to changes in fluid composition and pressure for
Marly and Vuggy zones. Based on isotropic rock physics model from ultrasonic core testing and
geophysical logs. Porosity values are based on well 02042364. ........................................................ 65
Table 4.5 Comparison of sensitivity of seismic attributes to changes in fluid composition and pressure for
Marly and Vuggy zones. Based on isotropic rock physics model from ultrasonic core testing and
geophysical logs. Porosity values are averages for RCP study area................................................... 69
Table 5.1 Description of fractures in Marly and Vuggy zones (from Bunge, 2000). .................................. 72
Table 5.2 Comparison of sensitivity of seismic attributes to changes in fluid composition and pressure for
Marly and Vuggy zones. Based on saturated, anisotropic rock physics model. Porosity values are
averages for RCP study area. .............................................................................................................. 87
Table 7.1 Description of components in Weyburn reservoir simulation. .................................................. 103
Table 7.2 Composition of live Weyburn oil modeled in STRAPP. ........................................................... 105
Table 7.3 Optimum a’i, b’i coefficients for Eq. (6.5) to use in cubic EOS for Weyburn oilCO2 system.
Coefficients correspond to units of L, MPa, K for molar volume, pressure and temperature (R=
0.008314 L MPa/(Kelvin gmole)).................................................................................................... 113
Table 8.1 Summary of reservoir properties in the reservoir simulation model.......................................... 117
Table 8.2 Summary of output pressure and saturation data from the reservoir simulator. Baseline refers to
the time of the first seismic survey before CO2 injection and repeat means the time of the monitor
seismic survey, approximately one year after the commencement of CO2 injection. Gas denotes the
presence of a CO2based phase.......................................................................................................... 122
Table 9.1 Comparison of baseline and monitoring seismic surveys. ......................................................... 152 xviii Table A.1 Composition of gas samples from Fluids Project. Values are mole fraction. .......................... 179
Table A.2 Composition of gas and condensate pseudosamples 1 to 10 modeled in STRAPP. Values are
mole fraction. .................................................................................................................................... 179
Table A.3 Composition of gas and condensate pseudosamples 11 to 18 modeled in STRAPP. Values are
mole fraction. .................................................................................................................................... 181
Table A.4 Composition of oil pseudosamples 1 to 10 modeled in STRAPP. Values are mole fraction. 183
Table A.5 Composition of oil pseudosamples 11 to 20 modeled in STRAPP. Values are mole fraction.
........................................................................................................................................................... 184
Table A.6 Optimal ai, bi, ai’, and bi’ coefficients for Eq. (A.10). Coefficients correspond to units of L,
MPa, K for molar volume, pressure and temperature........................................................................ 188
Table A.7 Sample values for calibration of cubic EOS code..................................................................... 190 xix ACKNOWLEDGMENTS This research was funded by the Reservoir Characterization Project and Rock Physics
Lab at Colorado School of Mines. I am grateful for their support and for the sponsors of these
consortia. Dr. Tom Davis, Dr. Mike Batzle, and Dr. Ilya Tsvankin have been excellent advisors
and I appreciate their guidance and confidence in my ability to succeed as a student.
This work is the result of collaboration with many fellow students, including Hiro
Yamamoto (reservoir simulation), Ida Herawati (seismic inversion), Reynaldo Cardona (rock
physics), and Hans Ecke (Perl programming). Their stimulating conversations and help in
learning have made graduate school a more enjoyable process.
I would also like to acknowledge Phillips Petroleum Corporation and Vaughn Ball for
their financial support in developing the cubic equation of state model for hydrocarbon fluids and
granting permission to use it in this research. Sandy Graham at PanCanadian Petroleum was
particularly helpful in providing necessary data for Weyburn Field. Reginald Beardsley of
Unocal offered helpful advice and data for fluids modeling. Marcia Huber at the National
Institute of Standards went the extra mile in providing technical support for the STRAPP
program. xx For Kristi, Hannah and Katie xxi 1 CHAPTER 1
INTRODUCTION 1.1 Introduction
This research is sponsored by the Reservoir Characterization Project (RCP) and the Rock
Physics Lab, industrysponsored research consortia at the Colorado School of Mines. The main
goal of RCP is to integrate dynamic data from time lapse (4D) multicomponent (9C) seismic
surveys into reservoir simulation. Phases VIII and IX of RCP focus on Weyburn Field,
Saskatchewan, a carbonate reservoir undergoing a CO2 flood. The seismic data acquired during
injection and production are integrated with geological and engineering data to dynamically
characterize the reservoir.
During enhanced oil recovery operations at Weyburn Field, many questions need to be
answered. Where are the injected fluids going? What is the sweep efficiency? Are there areas
with bypassed oil? Are some areas more heavily fractured and more permeable than others?
How can the injection scheme be optimized? To help answer these questions, a timelapse multicomponent seismic study is being conducted in a section of the CO2 flood area at Weyburn Field.
Rock physics theory is the link between the seismic data and the reservoir processes. The
purpose of this work is to establish a theoretical and experimental basis for identifying fluid
composition and fluid pressure changes in the reservoir from changes in the seismic data over
time. First, the necessary data to model these effects are acquired from ultrasonic testing of
reservoir fluids and rock core samples, and analysis of geophysical logs. Second, these data are
integrated to develop a pressuredependent, anisotropic rock physics model for the seismic
response of a saturated, fractured, rockfluid system. Third, this rock physics model is used to
evaluate the sensitivity of the seismic properties of the reservoir to changes in fluid composition
and pressure and to estimate the expected changes in the 4D seismic data based on the reservoir
simulation results. Fourth, the timelapse seismic data are compared with the timelapse seismic
attributes predicted from reservoir simulation. The data are interpreted for changes in fluid 2 pressure and composition due to the enhanced oil recovery (EOR) operations. The information
gained can be used to update the geological and reservoir models for Weyburn Field.
In this introductory section, Weyburn Field is described from both an oil production
(Section 1.2) and geologic (Section 1.3) perspective. The EOR scheme is discussed in Section
1.4. The RCP Phase IX study is introduced in Section 1.5. 1.2 Weyburn Field
Weyburn Field is located in the northern part of the Williston Basin in Saskatchewan,
Canada (Figure 1.1). It is operated by PanCanadian Petroleum. Production is from Mississippian
age carbonates at depths of 1300 to 1500 meters. The Midale reservoir beds are divided into two
main zones, an upper Marly dolostone zone and the lower Vuggy limestone zone. These are
described in Section 1.3. Figure 1.1 Location of Weyburn Field in Saskatchewan, Canada. 3 The following development and production data are taken from PanCanadian (1997).
Weyburn Field was discovered in 1954 and produced on primary production until 1964 when
waterflood was implemented. Production peaked in 1965 at 46,000 barrels/day. Because of the
fractured nature of the Vuggy unit, it was preferentially swept in the waterflood. Horizontal infill
drilling began in 1991 to target bypassed oil in the Marly unit. In 2000, a CO2 injection project
began, which is discussed in more detail in Section 1.4.
The objective of this enhanced oil recovery project is to increase the recovery rate for
Weyburn Field. The original oil in place at Weyburn Field is estimated at 1.4 billion barrels
(PanCanadian Petroleum, personal communication, 2002). As of 2000, approximately 24% of
this oil has been recovered. 1.3 Weyburn Field Geology and Reservoir Properties
The Williston Basin contains sediments of shallow marine origin and Cambrian to
Tertiary age. The Midale beds of the Mississippian Charles Formation were formed during a
transgression – regression sequence. The Bakken Shale is a possible source rock for the medium
gravity crude oil at Weyburn Field. The Midale reservoir beds were deposited in a shallow
carbonate shelf environment. The petroleum trap is both hydrodynamic and stratigraphic
(Churcher and Edmunds, 1994).
A detailed geologic description of Weyburn Field is found in Churcher and Edmunds
(1994), Bunge (2000), and Reasnor (2001). A stratigraphic column representative of Weyburn
Field is shown in Figure 1.2. The Marly beds consist of evaporite and carbonate units. The
carbonate reservoir has been subdivided into the Marly dolostone and the Vuggy limestone.
The following reservoir description and data are from PanCanadian (1997). The Marly
zone consists of chalky, microcrystalline dolostone and dolomitic limestone. These are
commonly separated by tighter, fractured limestone interbeds. In some areas, the Marly
carbonate beds are partially replaced by mudstones, grainstones, and dolomitic muds, formed as a
tidalchannel fill sequence. Net pay thickness in the Marly zone is 0.1 m to 9.8 m, with an
average of 4.3 m. Porosity ranges from 16% to 38%, and averages 26%. Permeability varies
from 1 md to over 100 md, with an average of 10 md. 4 Formation Lithology Age Mississippian Strata Watrous
Charles
Mission Canyon Jurassic Jurassic
Big Group Mississippian Dawson Bay
Elk Point Middle Souris River
Prairie
Evaporite
Winnipegosis Silurian Ashern
Interlake
Stonewall
Stony Mtn. Ordovician Red River Cambrian Charles
Formation Duperow Mission Canyon
Formation Bakken
Three Forks
Birdbear Madison Group Devonian Upper Lodgepole Winnipeg
Deadwood Lodgepole
Formation Watrous Formation
Kibbey Formation
Poplar Beds
Ratcliffe
Quengre Evap.
Beds
Midale Evap.
Midale
Frobisher Evap.
Beds
Hastings Evap.
Frobisher
Beds
Winlaw Evap.
Kisbey Sandstone
Gainsborough Evap. Alida Beds Tilston Beds
Souris Valley Beds
Bakken Formation PreCambrian Figure 1.2 Stratigraphic column for Weyburn Field. Left side is after Dietrich and Magnusson
(1998). Right side is after Wegelin (1984). The Vuggy zone is subdivided into two rock types formed in the shoal and inter shoal
environments. The shoals were higher energy environments where coarse grained carbonate
(limestone) sands accumulated. These sediments, more prevalent in the west end of Weyburn
Field, have porosity from 12% to 20%, with an average of 15%. Permeability is very high, from
10 md to over 500 md, with an average of 50 md.
Between the shoals, muddy carbonate sediments accumulated in a low energy shallow
marine environment. Porosity and permeability in the intershoal sediments are lower than in the
shoal sediments. Porosity ranges from a few % to 12%, with an average of 10%. Matrix
permeability varies from 0.1 md to 25 md, averaging 3 md.
Overall, net pay in the Vuggy zone ranges from 0.1 m to 18.6 m, with an average of 6 m.
Net porosity ranges from 8% to 20%, averaging 11%. The average porosity in the RCP Study 5 area is 10% (Section 8.2.1). Matrix air permeability values range from 0.3 md to over 500 md,
and average 15 md.
An evaporitic dolomite and shale sequence overlies the Midale reservoir and forms a top
seal. These beds are capped by the Midale Evaporite. The Ratcliffe and Poplar Beds, a series of
thin carbonate and evaporitic carbonate sequences, overlie the Midale Evaporite and are truncated
towards the north by the Mississippian unconformity. Beneath the Midale reservoir lie the
Frobisher Beds of the Mission Canyon Formation. The lithology and depositional environment
are similar to the Midale Beds. The original oilwater contact for the Weyburn reservoir is in the
upper part of the Frobisher Vuggy zone. 1.4 Enhanced Oil Recovery Operations
The following information is taken from PanCanadian (1997). The CO2 miscible flood
operation is expected to enhance oil recovery for several reasons. First, CO2 dissolves in and
significantly increases the volume of Weyburn oil. Waterflooding the swollen oil increases the
recovery by a maximum of 23%. Second, dissolved CO2 lowers the viscosity of oil and increases
its mobility in the reservoir.
In most other carbonate fields, CO2 breakthrough occurs before 5% pore volume has been
injected. This is due to large scale reservoir heterogeneities, gravity override of the CO2 at the
top of the formation, and lack of total miscibility of the CO2 and oil. A good way to overcome
early breakthrough problems is to alternate water and CO2 injection (WAG). This minimizes the
formation of the highly mobile CO2rich phase and maintains flow in the oilrich phase.
Simulation results suggest that four different approaches may be successful at Weyburn. These
are: simultaneous but separate injection of water and gas (SSWG), Vuggy water alternating gas
(VWAG), MarlyVuggy water alternating gas (MVWAG), and straight gas injection (SGI).
These methods will be used in different parts of the EOR project area and are expected to
increase oil recovery by more than 16% and extend the life of Weyburn Field by about 25 years.
In the RCP study area, the enhanced oil recovery strategy is SSWG.
In the simultaneous but separate water and gas injection scheme, CO2 is injected into the
Marly zone through horizontal wells. Water is injected into the Vuggy zone with vertical wells. 6 The density contrast should keep the CO2 in the Marly zone. Production is through both
horizontal and vertical wells. A schematic cross section of the EOR scheme is shown in Figure
1.3. The fluids in the reservoir are CO2, mixtures of CO2 and oil, original oil, and brine. The
seismic properties of these fluids are discussed in Chapter 2. Horizontal Producer Horizontal CO2 Injector Marly
CO2
Oil CO2 & Oil Vuggy
Water Vertical
Producer Vertical Water
Injector Figure 1.3 Schematic cross section of reservoir showing EOR strategy. Injection commenced in October 2000 for 19 patterns at Weyburn Field. The volume of
CO2 injected has ranged 3 to 7 million cubic feet per well per day. As of the date of the first
monitor survey (September 2001), 4% to 6% CO2 by pore volume had been injected into
reservoir in the RCP study area (Fall 2001 Sponsors Meeting presentation by David Stachniak). 1.5 RCP Phases VIII and IX Study
The RCP phases VII and IX study area is shown in Figure 1.4. It is approximately 9 km2.
A plan view of the EOR operations is shown in Figure 1.5. 7 Weyburn Field
R .14 R.13 R .12W 2
T.7 Edmo nton Regin a CO2 flood area Calgary T.6 RCP study area
T.5 Figure 1.4 Location of RCP 4D 9C seismic survey area. Modified from PanCanadian
Petroleum. 1 km Figure 1.5 EOR infrastructure at RCP 4D 9C seismic survey area. 8 During the CO2 and water injection, the reservoir is monitored with several timelapse
multicomponent (4D 9C) seismic reflection and vertical seismic profile surveys. The baseline
9C seismic reflection survey was shot immediately prior to CO2 injection in late September and
early October 2000. The first 9C monitoring survey was shot in September 2001. The seismic
surveys are described in Chapter 9. An additional monitoring survey is planned for 2003.
A main focus of RCP research is to integrate timelapse seismic data into reservoir
simulation. This thesis research contributes to that goal by extending reservoir simulation from
flow and compositional modeling to seismic attribute modeling. This is done by integrating the
reservoir model and compositional reservoir simulation with rock and fluid physics modeling.
This aids in interpretation of the timelapse seismic data and is useful for other RCP research in
fracture analysis, seismic inversion, and reservoir simulation. 9 CHAPTER 2
ACOUSTIC PROPERTIES OF RESERVOIR FLUIDS 2.1 Introduction
Estimates of the fluid saturation for the Marly and Vuggy zones at the times of the
baseline and monitor seismic surveys are used together with the acoustic properties of the fluid
phases to obtain the effective bulk moduli and density of the reservoir fluids. These data are
necessary to predict the changes in the seismic properties of the reservoir units due to changes in
fluid pressure, composition, and saturation. This chapter presents models for the acoustic
properties of reservoir fluids. The results are used in rock physics modeling (Chapters 4 and 5),
which will aid in the interpretation of the time lapse seismic data at Weyburn Field.
Fluid properties are used in calculating the elastic stiffness tensor of the rockfluid
system, as in Gassmann’s (1951) or Brown and Korringa’s (1975) equations. Fluids affect the
compressionalwave velocity, VP, of the rockfluid system by influencing the density, ρ, and bulk
modulus, K: VP = K+ 4
µ
3, (2.1) ρ where µ is shear modulus. Gassmann’s equation is a simple model for determining the seismic
velocities of a rock under different fluid conditions, a problem known as “fluid substitution.”
Gassmann theory is based on isostress conditions for an isotropic, homogeneous, monominerallic
rock at the low frequency limit. A common form is
2 K SAT = K DRY K
1 − DRY KM ,
+
φ
1 − φ K DRY
+
−2
K FL K M
KM where KSAT is saturated rock bulk modulus to be used in Eq. (2.1), KDRY is “dry” rock bulk (2.2) 10 modulus, KM is mineral modulus, φ is porosity, and KFL is fluid bulk modulus. “Dry” means
room dry. According to Gassmann’s theory, fluid has no effect on shear modulus, µ, and has a
minor influence on shearwave velocity, VS, through density VS = µ
.
ρ (2.3) To perform fluid substitutions using Gassmann’s or other equations, the density and bulk
modulus of the fluid must be modeled as a function of composition and pressure. 2.2 Weyburn Reservoir Fluids
At Weyburn Field, the enhanced oil recovery operations will change the composition and
saturation of the fluids in the Marly and Vuggy units. The fluids in these units are brine, original
oil, CO2, and mixtures of CO2 and oil. The following historical data on reservoir fluids are from
PanCanadian (1997). The average reservoir temperature is 63 °C. The reservoir pore pressure at
time of discovery is estimated as 14 MPa, the hydraulic gradient pressure. During primary
production, the pore pressure dropped locally to as low as 2 MPa and averaged over 6 MPa.
Under waterflood, the pore pressures vary from 8 to 19 MPa.
Recent pore pressure data are available from reservoir simulation (Section 8.2) and from
a well pressure survey performed by PanCanadian during 2001. The pressure survey shows that
the pore pressures in the RCP study vary from approximately 12.5 MPa to 18 MPa and average
15 to 16 MPa. For modeling purposes the expected pore pressure range is assumed to be 8 to 20
MPa under CO2 flood, with an average of 15 MPa. Higher pressures may occur close to the
injector wells, but pressure should decrease a short distance away from injectors, due to the
fractured and permeable nature of the reservoir. 2.2.1 Brine
According to PanCanadian reports (1997), the Midale reservoir units had an original
salinity of 229,000 ppm total dissolved solids (TDS). After nearly approximately 25 years of
waterflooding with water from the Blairmore Fm. (20,000 ppm TDS), the salinity has been 11 lowered to 85,000 ppm TDS in the reservoir. Given the nonuniformity of the water flood, it is
likely that salinity varies within the Marly and Vuggy units.
Brine is the simplest reservoir fluid to model, with a composition specified by one
variable, salinity. Fluid Acoustics for Geophysics (FLAG 4), a modeling program developed in
the Fluids Project at the Colorado School of Mines and Houston Advanced Research Center,
models the bulk modulus and density of brine as a function of salinity, temperature and pressure.
NaCl is the only dissolved solid considered. The FLAG relations for acoustic properties of brine,
oil, and gas are described in Batzle and Wang (1992) and Han and Batzle (2000b). Figures 2.1
and 2.2 show the bulk modulus and density, respectively, of brine as a function of salinity and
(pore) pressure, at an average reservoir temperature of 63 °C.
Brine is the most dense and least compressible fluid in the reservoir. Its bulk modulus
and density increase with increasing pressure and salinity. For expected reservoir conditions of
85,000 ppm TDS and pore pressures from 8 to 20 MPa, the bulk modulus of reservoir brine is
2.83 to 2.92 GPa. The density ranges from 1.044 to 1.049 g/cc. o Brine, T = 63 C
3.4
45k ppm NaCl
85k ppm NaCl
125k ppm NaCl Bulk Modulus, GPa 3.2 3.0 2.8 2.6 2.4
0 5 10 15 20 25 30 Pressure, MPa Figure 2.1 Bulk modulus of brine as a function of salinity and pressure, at reservoir temperature.
Modeled with FLAG 4.0. 12 o Brine, T = 63 C
1.10
45k ppm NaCl
85k ppm NaCl
125k ppm NaCl Density, g/cc 1.08 1.06 1.04 1.02 1.00
0 5 10 15 20 25 30 Pressure, MPa Figure 2.2 Density of brine as a function of salinity and pressure, at reservoir temperature.
Modeled with FLAG 4.0. For expected reservoir conditions, the acoustic properties of brine can be represented as
linear functions of pressure: K = 7.424 x10 −3 P + 2.773 (2.4) ρ = 4.032 x10 −4 P + 1.041 , (2.5) where K is GPa, P in MPa and ρ in g/cc. 2.2.2 Oil
Weyburn crude is a light to medium oil. PanCanadian (1997) reports the average oil
properties as 29 API gravity, 30 L/L gas oil ratio (GOR), and bubble point pressure of 6 MPa.
The range is 25 to 34 API gravity, 17 to 32 L/L GOR, and 2.2 to 6.7 MPa bubble point pressure.
Reported gas gravities (Gg) range from 0.935 to 1.3. For Weyburn oil, the average Gg is assumed
to be 1.22, the value of a sample tested in the Rock Physics Lab at Colorado School of Mines. 13 FLAG specifies the composition of live oil (oil with dissolved gas) by API gravity, GOR,
and Gg. FLAG Hydrocarbon Model 1 was used for bulk modulus and density calculations.
Figures 2.3 and 2.4 present the bulk modulus and density, respectively, of live Weyburn crude oil.
The bulk modulus and density of oil increase with decreasing API gravity, decreasing
GOR, and increasing pressure. Gas gravity affects the amount of gas that can dissolve in oil, but
has little influence on the properties of live oil. A dramatic decrease in bulk modulus and density
occurs below the “bubble point” pressure, the pressure at which gas begins to come out of
solution in the oil. Since gas is several orders of magnitude more compressible than oil at low
pressures, a small amount of gas will control the bulk modulus of the twophase mixture. As
shown in Figure 2.3, this occurs below 4 to 5 MPa pressure in the model for Weyburn oil. The
density also decreases significantly as shown in Figure 2.4, but not as much as bulk modulus
does, because density is the arithmetic average of the density of the oil and gas phases.
The enhanced oil recovery operations are designed to keep the reservoir pressure above
the bubble point. For the average original Weyburn oil and an expected pore pressure range of 8
to 20 MPa, the bulk modulus is 1.25 to 1.41 GPa, and the density is 0.826 to 0.834 g/cc. o Weyburn Crude Oil, T = 63 C Bulk Modulus, GPa 2.0 1.5 1.0 API 25, GOR 17 L/L, Gg 1.22
API 29, GOR 30 L/L, Gg 1.22
API 34, GOR 32 L/L, Gg 1.22 0.5 0.0
0 5 10 15 20 25 30 Pressure, MPa Figure 2.3 Bulk modulus of Weyburn crude oil as a function of composition and pressure, at
reservoir temperature. Modeled with FLAG 4.0. 14 o Weyburn Crude Oil, T = 63 C
0.9 Density, g/cc 0.8 0.7 0.6
API 25, GOR 17 L/L, Gg 1.22
API 29, GOR 30 L/L, Gg 1.22
API 34, GOR 32 L/L, Gg 1.22 0.5 0.4
0 5 10 15 20 25 30 Pressure, MPa Figure 2.4 Density of Weyburn crude oil as a function of composition and pressure, at reservoir
temperature. Modeled with FLAG 4.0. 2.2.3 CO2
At reservoir conditions CO2 is above its critical point of 7.4 MPa and 31 °C. It is a dense
fluid without a distinct phase transition between liquid and gas. The optimum bottomhole
injection pressure for CO2 at Weyburn Field is 24 MPa (PanCanadian, 1997), although this
pressure may not be realizable throughout the enhanced oil recovery operations. Because of the
small amount of CO2 injected compared to the reservoir volume, the expected temperature
changes in the reservoir are negligible.
The acoustic properties of pure CO2 can be modeled using STRAPP (NIST, 1988, 2002).
STRAPP is a program developed by the National Institute of Standards and Technology to model
thermophysical properties of hydrocarbon mixtures, using an equation of state. Figures 2.5 and
2.6 show the bulk modulus and density, respectively, of CO2 at reservoir conditions, modeled
with STRAPP. For bottomhole pressures of 18 to 24 MPa, the bulk modulus of pure CO2 is 0.05
to 0.17 GPa, and density ranges from 0.66 to 0.76 g/cc. 15 o CO2, T = 63 C 0.30 Bulk Modulus, GPa 0.25
0.20
0.15
0.10
0.05
0.00
0 5 10 15 20 25 30 Pressure, MPa Figure 2.5 Bulk modulus of CO2 as a function of pressure, at reservoir temperature. Modeled
with STRAPP. o CO2, T = 63 C Density, g/cc 0.8 0.6 0.4 0.2 0.0
0 5 10 15 20 25 30 Pressure, MPa Figure 2.6 Density of CO2 as a function of pressure, at reservoir temperature. Modeled with
STRAPP. 16 2.2.4 Oil + CO2 Mixtures
The injected CO2 will not form a pure, separate phase at reservoir conditions. It will
extract hydrocarbons from the oil until it attains a composition that is miscible with the oil above
a minimum miscibility pressure of 14 to 17 MPa (PanCanadian, 1997). Up to approximately 0.66
mole fraction CO2 can dissolve in the oil above the minimum miscibility pressure. If more CO2 is
present, a CO2 rich phase with dissolved light hydrocarbons will be formed. However, since this
phase will be more mobile than the oil, it will flow faster, contact fresh oil, and dissolve to
saturation levels in the oil.
The acoustic properties of oil and CO2 mixtures can be modeled using STRAPP and
FLAG. The addition of CO2 to live oil has the effect of increasing the GOR and the gas gravity.
The mole percent CO2 is calculated based on an estimated average molecular weight of 320 for
the dead oil samples tested for CSM. STRAPP is used to calculate the volume (for GOR) and gas
gravity of the CO2 and hydrocarbon gas. These are then input into FLAG and the bulk modulus
and density can be calculated.
Because the empirical equations in FLAG were optimized for oils with dissolved
hydrocarbon gases, it is necessary to verify FLAG’s accuracy for modeling dissolved CO2. A
sample of live Weyburn oil with 38% mole CO2 was tested in the lab and compared to modeling
results using STRAPP and FLAG. The measurement procedure is described in Han and Batzle
(2000a) and the Fluid Project Report (2000). Hydrocarbon Model 1 fit the data best, as shown in
Figure 2.7. For practical reasons, the test was carried out at room temperature and Pwave
velocity was measured instead of dynamic bulk modulus. The misfit is generally less than 2%.
Based on the laboratory test results, STRAPP and FLAG are adequate for modeling CO2
dissolved in oil. For modeling, the average properties of 29 API gravity, GOR of 30 L/L, gas
gravity of 1.22 were used as the base Weyburn oil. The modeled bulk modulus and density for
mixtures of oil and CO2 are shown in Figures 2.8 and 2.9 respectively. The end members of
original live oil and pure CO2 are included for comparison. In Figure 2.9, the apparent higher
density of pure CO2 compared to oil + 66% CO2 is not real. CO2 properties were calculated in
STRAPP, a highly accurate equation of state based program, whereas oil + CO2 mixtures were
modeled with the empirical FLAG relations, and may be less accurate. 17 o Live Weyburn Oil + 38% CO 2, T ~ 22.5 C PWave Velocity, m/s 1600 Lab Measurements
Flag Hydrocarbon Model 1 1500 1400 1300 1200
0 10 20 30 40 50 60 70 Pressure, MPa Figure 2.7 Comparison of ultrasonic measurement results of Pwave velocity with FLAG model. o Weyburn Oil + CO2 mixtures, T = 63 C
1.6 Bulk Modulus, GPa 1.4
1.2
1.0 Oil
Oil + 10% CO2
Oil + 25% CO2
Oil + 40% CO2
Oil + 66% CO2
CO2 0.8
0.6
0.4
0.2
0.0
0 5 10 15 20 25 30 Pressure, MPa Figure 2.8 Bulk Modulus of oil and CO2 mixtures as a function of composition and pressure, at
reservoir temperature. Modeled with STRAPP and FLAG. 18 o Weyburn Oil + CO2 mixtures, T = 63 C Density, g/cc 0.8 0.6
Oil
Oil + 10% CO2
Oil + 25% CO2
Oil + 40% CO2
Oil + 66% CO2
CO2 0.4 0.2 0.0
0 5 10 15 20 25 30 Pressure, MPa Figure 2.9 Density of oil and CO2 mixtures as a function of composition and pressure, at
reservoir temperature. Modeled with STRAPP and FLAG. The addition of CO2 to live Weyburn oil lowers the bulk modulus and density and raises
the bubble point pressure. However, the results below the modeled bubble point pressure may be
inaccurate because CO2 would not exsolve as a gas, as assumed in FLAG, but would form a CO2
phase rich with dissolved light hydrocarbons, as discussed previously. The behavior of the bulk
modulus below the bubble point pressure or minimum miscibility pressure would likely be
smoother than the FLAG results.
The acoustic properties of the phase consisting of CO2 with dissolved light hydrocarbons
could be modeled using STRAPP if the composition were known. Since data on the range in
composition of this phase were not available, it was not modeled. The properties would likely
fall between those of pure CO2 and Oil + 66% CO2. Based on reservoir simulation (Chapter 8),
this phase is expected to be present close to the CO2 injection well. Properties of and options for
modeling the CO2based phase are discussed in Chapter 7.
The expected bulk modulus for oil with dissolved CO2 is 0.6 to 1.4 GPa at pore pressures
of 15 to 24 MPa. The expected density is 0.675 to 0.83 g/cc. A significant amount of dissolved 19 CO2 (25% +) is necessary for the acoustic properties of oil and CO2 mixtures to be different than
those for the range in original oil compositions (compare Figures 2.3 and 2.8, 2.4 and 2.9). 2.3 Summary
The generalized bulk modulus and density relations for Weyburn reservoir fluids under
expected reservoir pressure and temperature conditions are presented in Figures 2.10 and 2.11.
The properties of these homogeneous fluids can be used for fluid substitutions to determine the
maximum fluid effects on the seismic properties of a saturated rock. In the reservoir, however,
the pores may contain brine, oil, and CO2rich phases. To determine the effective properties of
this heterogeneous fluid mixture, a mixing law must be applied, such as Wood’s (1955) relation.
Wood’s relation is based on isostress behavior of an isotropic, linear, elastic material.
The effective bulk modulus, KEFF, is expressed as K EFF = 1
n f
∑ Ki
i =1
i , (2.6) and the average density, ρ, is
n ρ = ∑ fi ρi , (2.7) i =1 where fi, Ki, and ρi are the volume fraction, bulk moduli, and densities of the phases, respectively.
By combining the relations for acoustic properties of single fluid phases, and the fluid
saturation and pressure estimates from reservoir simulation, the effective fluid density and bulk
modulus can be calculated for the Marly and Vuggy zones. These calculations can be done for
the reservoir conditions at the baseline and repeat seismic surveys. The data are used in rock
physics modeling to aid in interpretation of the timelapse seismic data for changes in reservoir
fluid saturation and pressure. 20 o Weyburn Reservoir Fluids, T = 63 C
3.0
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 25% CO2
Oil + 40% CO2
Oil + 66% CO2
CO2 Bulk Modulus, GPa 2.5
2.0
1.5
1.0
0.5
0.0
6 8 10 12 14 16 18 20 22 24 26 Pressure, MPa Figure 2.10 Bulk modulus of reservoir fluids as functions of composition and pressure, at
reservoir temperature. Modeled with STRAPP and FLAG. o Weyburn Reservoir Fluids, T = 63 C
1.0 Density, g/cc 0.8
0.6
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 25% CO2
Oil + 40% CO2
Oil + 66% CO2
CO2 0.4
0.2
0.0
6 8 10 12 14 16 18 20 22 24 Pressure, MPa Figure 2.11 Density of reservoir fluids as functions of composition and pressure, at reservoir
temperature. Modeled with STRAPP and FLAG. 26 21 CHAPTER 3
ULTRASONIC MEASUREMENTS ON CORE SAMPLES AND DRY ROCK
ACOUSTIC PROPERTIES 3.1 Introduction
Ultrasonic measurements were performed on core samples from the Marly and Vuggy
units for two main purposes. First, the dry rock elastic moduli are necessary for fluid
substitutions, i.e. Gassmann’s theory (Eq. 2.2), or anisotropic rock physics modeling (Chapter 5).
Dry rock moduli can be estimated from geophysical logs but there is some uncertainty in fluid
saturation and pressure. Through lab testing, the pressure dependence of dry rock moduli can be
determined. This is necessary to model the effects of changing pore pressure during the enhanced
oil recovery operations.
The second purpose for laboratory measurements in this work is to measure velocity
changes during a simulation of the CO2 flood. A core sample from the Vuggy unit was tested
with brine, oil with dissolved CO2, and CO2 saturation. From these test results, the magnitude of
the expected changes in seismic properties of this reservoir zone can be estimated. The results
can also be compared to modeling results to test the applicability of models such as Gassmann’s
theory. From test results, seismic attributes are identified that have potential for discriminating
seismic effects of pore fluid composition and pressure change.
Although the elastic properties of a rock can be measured very precisely through
ultrasonic measurements, there are several drawbacks. The core sample does not contain the
largescale reservoir features, such as fractures, that influence both fluid flow and seismic
response. This limitation is shared to some extent with geophysical logs. The information gained
from core testing is still useful as representative of the background rock matrix in the reservoir, to
which fractures can be added for modeling. The frequency range of ultrasonic testing (~ 500
kHz) does not correspond to the seismic bandwidth (5 to 100 Hz). In saturated rocks, seismic
velocities may exhibit significant dispersion over this frequency range, depending on the
viscosity of the fluid and the permeability of the rock. In the ultrasonic testing system used in the 22 CSM rock physics lab, the confining stress is isotropic. The confining and differential stresses in
the reservoir are not isotropic, but the uncertainty in the magnitude of the horizontal stresses in
the reservoir justifies lab testing at an estimated equivalent isotropic stress state. 3.2 General Fluid Substitution Theory
The effect of fluids on the seismic properties of a rockfluid system is generally known as
fluid substitution. Fluids affect seismic velocities in several ways. Increasing the fluid bulk
modulus increases the bulk modulus of the rockfluid system, tending to increase the Pwave
velocity. Less compressible fluids are usually also denser. Increasing fluid density increases the
density of the porous fluidrock system and tends to lower P and Swave velocities. Changes in
pressure cause additional complications. Fluid bulk modulus and density both increase with
increasing pore pressure and cause changes mentioned previously. If the pressure change causes
a phase change in the pore fluid, the effects can be large. At constant confining (overburden)
pressure, increasing the pore pressure decreases the differential pressure on the rock skeleton (Eq.
3.1) and lowers its dry rock moduli, tending to lower P and Swave velocities.
Differential pressure, Pd, is the difference between confining pressure, Pc, and pore
pressure, PP, Pd = Pc − Pp . (3.1) Confining pressure is often called ‘total’ pressure or stress in engineering disciplines. Pressure
refers to an isotropic or hydraulic stress. Stresses can be arbitrarily anisotropic. The differential
pressure is usually equivalent to the ‘effective’ pressure or stress.
The magnitude of the changes due to fluid saturation and pressure are determined by the
exact rock and fluid properties and reservoir process. As an example, in a waterflood of an oil
reservoir, the Pwave velocity will tend to increase because the bulk modulus of water is usually
higher than that for oil. If pore pressures are maintained or increased, Pwave velocity may
decrease due to decreasing differential pressure. The total change is the sum of these opposing
effects. Figure 3.1 shows these processes. 23 5100 PWave Velocity, m/s Brine
Oil
5000 Net change
Increase SW
4900 Increase pore pressure
4800 Increasing differential pressure
4700
6 8 10 12 14 16 18 20 Pore Pressure, MPa Figure 3.1 Example of change in Pwave velocity of reservoir rock under waterflood due to
changes in fluid saturation and pressure. 3.3 Ultrasonic Testing Procedure
The measurement procedure is described in Castagna et al. (1993) and summarized here. 3.3.1 Sample Preparation
The sample is prepared from precleaned 1 ½” core supplied by PanCanadian Petroleum.
The samples are identified by well number and depth, which can be compared to geophysical logs
to identify the geologic unit. The top and bottom of the sample are cut and ground flat and
perpendicular to the plug axis. Sample dimensions are measured, from which the volume is
computed. Both the dry sample mass and the watersaturated sample mass are measured. From
the sample volume and mass, dry density is calculated. Porosity fraction, φ, is calculated from
dry density, ρdry, saturated density, ρsat, and fluid density, ρfl φ= ρ sat − ρ dry
,
ρ fl (3.2) 24 where the units for ρsat, ρdry, ρfl are the same, i.e. g/cc. Grain density, ρgrain, is calculated from dry
density and porosity, ρ grain = ρ dry
,
1−φ (3.3) and is helpful in determining grain composition. Permeability is measured with an air
permeameter. 3.3.2 Ultrasonic Testing Equipment
A schematic of the ultrasonic testing equipment setup is shown in Figure 3.2. The
sample and lower part of the displacement transducers are enclosed in shrink wrap or Tygon
tubing, and clamped with metal wire. The transducer casing also has pore fluid lines. The
transducer wires and pore fluid lines are connected to the ultrasonic testing apparatus. The
sample assembly is enclosed in a metal confining vessel, so that high confining pressures can be
applied across the sample membrane. Figure 3.2 Schematic of ultrasonic testing setup for saturated cores. 25 Pore pressure is controlled independently by a digital pump. A transfer vessel with
floating piston allows water to be used in the pump but the fluid in the sample to be easily
changed, i.e. from brine to oil or CO2. The volume of fluid flowing into the sample is measured
with the digital pump and the output volume is measured in a graduated cylinder. Both pore and
confining pressures can range from 0 psi to 10,000 psi (69 MPa).
A pulse generator is used to send a signal to the transmitting piezometric transducer and
to trigger the digital oscilloscope used for recording the output signal from the second receiving
transducer (Figure 3.3). The frequency of the resulting signal is in the range of 500,000 Hz. The
transducer wires are connected to a junction box so that the input signal and measurement axis
can be varied from Pwave (particle motion parallel to long axis of sample) and two orthogonal
shearwave directions (particle motion perpendicular to long axis of sample). The wave
velocities are calculated from the travel time of the pulse through the sample, as measured by the
corrected arrival time at the second transducer (tC), and the length of the sample (L), V= L
.
tC (3.4) The arrival time must be corrected for by the zerotime calibration of the transducer pair. The
transducer waveforms can be recorded by a PC running LabView and analyzed later if there is
uncertainty in arrival times. Figure 3.3 Recording setup for ultrasonic testing 26 3.3.3 Ultrasonic Measurements
Once the sample has been prepared and placed in the measurement apparatus, pressure
lines, transducer wires, and signal cables are connected. For dry measurements, the pore fluid
lines are left open to the air. The confining pressure is raised incrementally for each
measurement. After the highest confining pressure is reached, the measurements are repeated,
decreasing the pressure incrementally to the starting pressure. This shows the stress hysteresis of
the sample.
At each confining pressure, the sample is allowed to equilibrate for approximately five
minutes. The confining pressure is recorded from a digital readout. Pwave and two Swave
measurements are made sequentially. The two shearwave measurements, VS1 and VS2, are
orthogonal but not necessarily oriented in the principle directions related to core fabric. On the
received waveforms, the times of the first break, first trough and first peak are recorded.
Depending on the sample velocity, the arrival times are recorded to a precision of 0.02 or 0.05 µs.
For typical Weyburn samples with lengths of 4 to 5 cm, this precision corresponds to an error of
less than 0.5% for Pwave velocity (VP) and 0.3 % for Swave velocity (VS). Optionally, the
waveforms are also recorded using LabView. Because of attenuation, the first break may be the
most accurate measurement but is often obscured by noise. Heterogeneity is usually the largest
source of error. The average of the calibrationcorrected first break, trough, and peak times gives
a stable estimate and reduces the random error in VP and VS.
For saturated core measurements, the measurement procedure is largely the same, with
the additional variable of pore pressure. The pore pressure input line is connected to a fluid
pressure vessel. The confining and pore pressures are usually kept above the bubble point of the
test fluid. Five to ten pore volumes of fluid are flowed through the sample, as measured in the
graduated cylinder at the pore fluid output line. The pore output line is then closed so that any
remaining gas will dissolve in the test fluid and the sample will be fully saturated.
Differential pressure controls the stiffness of the rock frame in a fluidrock system. To
insure that the core sample remains intact, differential pressure is always kept positive. Several
options exist for saturated core measurements, including constant differential pressure, in which
confining and pore pressures are changed by the same increment, and constant confining
pressure, in which only the pore pressure is varied. The usefulness of each test is explained in 27 sections 3.5.1 and 3.5.2. For saturated core measurements, the sample is allowed to sit for 10 to
15 measurements between measurements so that pore pressures can equilibrate. 3.4 Dry Rock Ultrasonic Measurements
Dry rock measurements are much simpler and quicker to perform than saturated rock
measurements. The pore pressure is zero (atmospheric), so the confining pressure is equal to the
differential pressure. From measurements of VP and VS for a range of confining pressures, the
dependence of dry rock bulk modulus, KDRY, on differential pressure can be determined by 4 K = ρ VP2 − VS2 .
3 (3.5) This information is necessary for pressuredependent fluid substitutions using Gassmann’s
equation (Chapter 4). The variation of dry rock density is very small for the Weyburn samples
tested, and was therefore not measured or applied in Eq. (3.5). In calculating isotropic elastic
moduli, VS is taken to be the average of VS1 and VS2.
For the rock physics modeling scheme used in this work, the variation of dry rock Lame’s
parameter, λDRY, and shear modulus, µDRY, with differential pressure are required. These can be
calculated by λ = ρ (VP2 − 2VS2 ),
µ = ρVS2 . (3.6)
(3.7) 3.4.1 Marly Unit
Two samples from the Marly unit were tested. They are described in Table 3.1. Figures
3.4 and 3.5 present the variation of VP and VS with differential pressure for samples PC 1213
4626 and PC 1411 4659.5, respectively. Dry rock bulk modulus as a function of pressure for
both Marly samples is shown in Figure 3.6. Figures 3.7 and 3.8 present the variation of λ and µ
with differential pressure for samples PC 1213 4626 and PC 1411 4659.5, respectively. 28 Table 3.1 Description of Marly unit samples.
Sample Name
(well, depth in ft) Porosity
(fraction) Permeability
(md) Dry Density
(g/cc) Grain Density
(g/cc) PC 1213 4626’ 0.29 21.7 1.93 2.88 PC 1411 4659.5’ 0.333 34.8 1.815 2.72 Dry Marly Sample PC 1213, 4626 ft
3600
2100 2000
3200
1900
3000
VP
VS1
VS2 2800 1800 1700 2600
0 5 10 15 20 25 30 35 Confining = Differential Pressure, MPa Figure 3.4 Variation of ultrasonic P and Swave velocities with pressure for dry Marly sample
PC 1213 4626. Hysteresis can be seen in the variation between the increasing versus
decreasingpressure measurements. SWave Velocity, m/s PWave Velocity, m/s 3400 29 Dry Marly Sample PC 1411, 4659.5 ft 3400
PWave Velocity, m/s 2200 2100 3200 2000 3000 SWave Velocity, m/s 3600 1900
VP
VS1
VS2 2800 1800 2600
0 5 10 15 20 25 1700
35 30 Confining = Differential Pressure, MPa Figure 3.5 Variation of ultrasonic P and Swave velocities with pressure for dry Marly sample
PC 1411 4659.5. Dry Marly Samples
17
16
Bulk Modulus, GPa 15
14
13
12
11
PC 1213, 4626 ft, φ = 0.29
PC 1411, 4659.5 ft, φ = 0.33 10
9
8
0 5 10 15 20 25 30 35 Confining = Differential Pressure, MPa Figure 3.6 Variation of dry rock bulk modulus with pressure for Marly samples PC 1213 4626
and PC 1411 4659.5. 30 Dry Marly Sample PC 1213, 4626 ft Lame's Parameter, GPa 11 10 10 9 9 8 8 7 λ
µ 6
0 5 10 15 20 25 30 Shear Modulus, GPa 11 7 6
35 Confining = Differential Pressure, MPa Figure 3.7 Variation of Lame’s parameter and shear modulus with pressure for dry Marly sample
PC 1213 4626. Dry Marly Sample PC 1411, 4659.5 ft 7.5
Lame's Parameter, GPa 8.0
7.5 7.0 7.0 6.5 6.5 6.0 6.0
λ
µ 5.5
5.0
0 5 10 15 20 25 30 Shear Modulus, GPa 8.0 5.5
5.0
35 Confining = Differential Pressure, MPa Figure 3.8 Variation of Lame’s parameter and shear modulus with pressure for dry Marly sample
PC 1411 4659.5. 31 Because of the relationship between seismic velocities and moduli (Eq. 3.5, 3.6, and 3.7),
errors are amplified from velocities to moduli. The dependence of moduli on porosity is shown
in the sample data (Fig. 3.6). The data on the pressure dependence of elastic moduli are
smoothed for modeling purposes. The estimated range in differential pressure in the reservoir is 2
to 14 MPa, as discussed in Section 4.2.2. From 2 to 14 MPa differential pressure, the change in
Lame’s parameter is +11% and +17% and the change in shear modulus is +18% for the two
Marly samples tested. 3.4.2 Vuggy Unit
The two samples from the Vuggy unit that were tested are described in Table 3.2.
Figures 3.9 and 3.10 show the variation of VP and VS with differential pressure for samples PC
1411 4698 and PC 1213 4643, respectively. Dry rock bulk modulus as a function of pressure
for both Vuggy samples is shown in Figure 3.11. Figures 3.12 and 3.13 present the variation of λ
and µ with differential pressure for samples PC 1411 4698 and PC 1213 4643, respectively. Table 3.2 Description of Vuggy unit samples.
Sample Name
(well, depth in ft) Porosity
(fraction) Permeability
(md) Dry Density
(g/cc) Grain Density
(g/cc) PC 1213 4643’ 0.115 1.3 2.305 2.60 PC 1411 4698’ 0.133 103 2.34 2.70 Due to testing difficulties, only VS1 could be measured on sample PC 1213 4643.
Sample PC 1411 4698 is more representative of the Vuggy reservoir units than is PC 1213
4643. PC 1213 4643 is the rock referred to as “Vuggy tight” in PanCanadian’s core testing
report (Core Laboratories, 1998). It is much stiffer, less permeable, and its elastic moduli have
little pressure dependence in comparison with sample PC 1411 4698. The properties of PC 1411 4698 are used for the Vuggy unit in rock physics modeling described in Chapter 4. 32 Dry Vuggy Sample PC 1411, 4698 ft
5100 2950
VP
VS1
VS2 2900 4900 2850 4800 2800 4700 SWave Velocity, m/s PWave Velocity, m/s 5000 2750 4600
0 5 10 15 20 25 30 2700
35 Confining = Differential Pressure, MPa Figure 3.9 Variation of ultrasonic P and Swave velocities with pressure for dry Vuggy sample
PC 1411 4698. Dry Vuggy Sample PC 1213, 4643 ft
5900 2800
VP
VS1 2750 5700 2700 5600 2650 5500 SWave Velocity, m/s PWave Velocity, m/s 5800 2600 5400
0 5 10 15 20 25 30 2550
35 Confining = Differential Pressure, MPa Figure 3.10 Variation of ultrasonic P and Swave velocities with pressure for dry Vuggy sample
PC 1213 4643. 33 Dry Vuggy Samples
60 Bulk Modulus, GPa 55
50
45 PC 1213, 4643 ft, φ = 0.11
PC 1411, 4698 ft, φ = 0.13 40
35
30
0 5 10 15 20 25 30 35 Confining = Differential Pressure, MPa Figure 3.11 Variation of dry rock bulk modulus with pressure for Vuggy samples PC 1213 4643
and PC 1411 4698. Dry Vuggy Sample PC 1411, 4698 ft
19 λ
µ 22 18 20
17
18 16
0 5 10 15 20 25 30 16
35 Confining = Differential Pressure, MPa Figure 3.12 Variation of Lame’s parameter and shear modulus with pressure for dry Vuggy
sample PC 1411 4698. Shear Modulus, GPa Lame's Parameter, GPa 24 34 Dry Vuggy Sample PC 1213, 4643 ft Lame's Parameter, GPa 19 44 18 42 17 40 16 Shear Modulus, GPa 46 λ
µ
38
0 5 10 15 20 25 30 15
35 Confining = Differential Pressure, MPa Figure 3.13 Variation of Lame’s parameter and shear modulus with pressure for dry Vuggy
sample PC 1213 4643. Over a differential pressure range of 2 to 14 MPa, the change in Lame’s parameter is 2%
and +10% and the change in shear modulus is +4% for the two Vuggy samples tested. Although
a decrease in Lame’s parameter with increasing differential pressure is unusual, it occurs when
the change in VS with pressure is more than one half of the change in VP (Eq. 3.6). 3.5 Saturated Rock Ultrasonic Measurements
Ultrasonic measurements were performed on Vuggy sample PC 1411 4698 under
different saturation conditions. This sample was selected because the Vuggy unit is the thicker of
the two reservoir zones and for that reason, changes within it may have the most effect on the
seismic data. The Marly zone shows greater fluid and pressure effects (Chapter 4) but is thinner
and changes within it may be more difficult to resolve seismically.
The ultrasonic measurements were performed first with the sample saturated with brine,
then 27 API Weyburn crude oil with 39% dissolved CO2, and finally pure CO2. Five to ten pore
volumes of each fluid were flowed through the sample before each testing phase, but there was 35 still residual fluid saturation (i.e. brine) during the oil and CO2 measurement phases. The
measurement of residual fluid saturation was imprecise and is therefore omitted here.
Measurements were made at confining pressures up to 34.5 MPa (5000 psi) and pore
pressures up to 27.6 MPa (4000 psi). The data were sorted into two groups for analysis: constant
differential pressure and constant confining pressure. 3.5.1 Constant Differential Pressure
The advantage of making ultrasonic measurements at a constant differential pressure is
that the “fluid effect” is readily seen. In this measurement scheme, confining pressure and pore
pressure are varied by equal amounts. Because the dry rock bulk and shear moduli depend on
differential pressure, which is held constant, differences in velocities measured in the sample
under different saturation are due to differences in bulk modulus and density of the pore fluid.
For a common fluid, differences in seismic velocities with changing pore pressure are due to
changes in the fluid bulk modulus and density with pressure. These measurements provide a
good estimate of the magnitude of the “fluid effect” on the rock.
Figures 3.14 and 3.15 show the variation of P and Swave velocity, respectively, for
saturated Vuggy sample PC 1411 4698 at a constant differential pressure of 6.9 MPa (1000 psi).
This pressure is close to the estimated average differential pressure for the reservoir of 7 MPa.
Variations in seismic velocities are similar at a differential pressure of 13.8 MPa, so those results
are omitted here.
As shown in Figure 3.14, the difference in Pwave velocity for the different reservoir
fluids is small, 1% to 2%. As expected the Pwave velocities are highest for the least
compressible reservoir fluid, brine, and lower for oil + CO2 and CO2. Below the bubble point
pressure of the oil + CO2, mixture, Pwave velocity does not distinguish CO2 from oil with CO2.
CO2 and oil + CO2 show greater variation in bulk modulus with pressure than brine does
and this is reflected in the greater change in Pwave velocity with pore pressure for the sample
saturated with these fluids. However, the changes are only slightly greater than the measurement
precision of 0.5%. 36 Saturated Vuggy Sample PC 1411 4698, Differential Pressure = 6.9 MPa
5100
Brine: 50k ppm NaCl
Oil + 39% CO2
CO2 PWave Velocity, m/s 5050
5000
4950 below
bubble point 4900
4850
4800
0 5 10 15 20 25 30 Pore Pressure, MPa Figure 3.14 Variation of ultrasonic Pwave velocity with pore pressure at constant differential
pressure for saturated Vuggy sample PC 1411 4698. Saturated Vuggy Sample PC 1411 4698, Differential Pressure = 6.9 MPa
Brine: 50k ppm NaCl
Oil + 39% CO2
CO2 SWave Velocity, m/s 2750 2700 2650 0 5 10 15 20 25 Pore Pressure, MPa Figure 3.15 Variation of ultrasonic Swave velocity with pore pressure at constant differential
pressure for saturated Vuggy sample PC 1411 4698. 30 37 For the Vuggy sample tested, there is little variation in shear wave velocity with fluid
type or pore pressure (Fig. 3.15). This is expected from Gassmann’s theory, which predicts that
fluid has no effect on shear modulus. Any change in Swave velocity would be due to difference
in density between the fluids and would be minor. At high pore pressure (27.5 MPa) this is seen
in the brinesaturated Swave velocity, which is about 1% lower than the CO2 or oil + CO2
saturated Swave velocity. Over the expected pore pressure range of 8 to 20 MPa, there is no
variation in Swave velocity within the precision (0.3%) of the ultrasonic test data. 3.5.2 Constant Confining Pressure
Ultrasonic measurements at constant confining pressure simulate the enhanced oil
recovery operations. In this measurement scheme, confining pressure is held constant and pore
pressure is varied. The confining pressure in ultrasonic testing is isotropic, and represents the
average of the vertical and horizontal “total” or confining stresses. In the reservoir, the confining
stress can be thought of as the vertical stress due to the overburden material and horizontal
stresses due to lateral support from adjacent material.
In the enhanced oil recovery operations, the vertical confining (overburden) stress in the
reservoir should remain approximately constant. Injecting fluids at high pressure will raise the
reservoir pore pressure and producing fluids may decrease the pore pressure. Raising the pore
pressure will also decrease the differential pressure on the rock, both of which will affect seismic
velocities, as explained in Section 3.2.
At a single confining pressure and pore pressure, the “fluid effect” is seen by the
difference in velocities for different pore fluids. For a common pore fluid and constant confining
pressure, measurements at different pore pressures show both the “fluid effect” and the rock
skeleton “pressure effect.” Figures 3.16 and 3.17 show the variation of P and Swave velocity,
respectively, for saturated Vuggy sample PC 1411 4698 at a constant confining pressure of 34.5
MPa (5000 psi).
A confining pressure of 34.5 MPa is probably higher than the average of horizontal and
vertical confining stresses in the reservoir (see Sections 4.2.2, 8.5). However, since
measurements were made for a wider range of pore pressures than expected in the reservoir, the
range of differential pressure in the reservoir is covered by the test data. 38 Saturated Vuggy Sample PC 1411 4698, Confining Pressure = 34.5 MPa PWave Velocity, m/s 5100 5050 5000 4950 4900
10 Brine: 50k ppm NaCl
Oil + 39% CO2
CO2
15 20 25 30 Pore Pressure, MPa Figure 3.16 Variation of ultrasonic Pwave velocity with pore pressure at constant confining
pressure for saturated Vuggy sample PC 1411 4698. Saturated Vuggy Sample PC 1411 4698, Confining Pressure = 34.5 MPa
2760 SWave Velocity, m/s 2740
2720
2700
2680
2660
10 Brine: 50k ppm NaCl
Oil + 39% CO2
CO2
15 20 25 Pore Pressure, MPa Figure 3.17 Variation of ultrasonic Swave velocity with pore pressure at constant confining
pressure for saturated Vuggy sample PC 1411 4698. 30 39 The difference in Pwave velocity at constant confining and pore pressures is 1% to 2%
for the Vuggy sample saturated with brine, oil + CO2, and CO2 (Fig. 3.16). Over a pore pressure
range of 10 MPa, the variation in Pwave velocity for the sample saturated with a common fluid
is approximately 1%. For the Vuggy sample tested, Pwave velocity does not uniquely
distinguish fluid saturation from differential pressure effects. For example, the Pwave velocity
for brine saturated rock at high pore pressure could be the same as oil + CO2 saturated sample at
low pore pressure.
As expected from Gassmann’s theory, there is little difference in Swave velocity for
different fluid saturations in the sample at the same pore and confining pressures (Fig. 3.17). The
Swave velocities for the brine saturated sample may be slightly lower because of the fluid
density effect (Eq. 2.3). Swave velocity varies by approximately 1% over a pore pressure range
of 10 MPa for the Vuggy sample tested. For an isotropic rock, Swave velocity does not
discriminate fluid saturation well, but is sensitive to pore pressure. 3.5.3 Seismic Attributes
One of the purposes of the ultrasonic velocity measurements is to identify seismic
attributes that may be useful in interpreting the timelapse seismic data for fluid and pore pressure
changes. Although there are major differences between the lab core samples and measurements
and the seismic data collected over the reservoir, the seismic attributes from the lab can at least be
qualitatively interpreted in terms of reservoir processes. More realistic modeling of the seismic
response of the reservoir can be done for these attributes, so that they can be more quantitatively
interpreted.
Because the P and Swave velocities of the Vuggy sample exhibit approximately the
same pore pressure dependence (Figs. 3.16 and 3.17), the VP/VS ratio removes the pore pressure
variation in seismic velocities. Changes in the VP/VS ratio should then show only differences in
velocity due to changes in fluid saturation. This is shown in Figure 3.18. The magnitude of the
difference between brine, oil + CO2, and CO2 saturation is 1% to 2%. These values are for the
Vuggy matrix only and do not take into account large scale fractures present in the reservoir. 40 Saturated Vuggy Sample PC 1411 4698, Confining Pressure = 34.5 MPa
1.88 VP/VS 1.86 1.84 1.82
Brine: 50k ppm NaCl
Oil + 39% CO2
CO2 1.80
10 15 20 25 30 Pore Pressure, MPa Figure 3.18 Variation of ultrasonic VP/VS with pore pressure at constant confining pressure for
saturated Vuggy sample PC 1411 4698. Velocities were measured for shear waves polarized in orthogonal directions, and are
denoted VS1 and VS2. If these directions correspond to the polarizations of the fast and slow
shearwaves, the shear splitting parameter for an anisotropic rock with transversely isotropic (TI)
symmetry can be calculated. As defined by Thomsen (1986) for TI media, VS21 − VS22
γ≡
,
2VS22 (3.8) where γ is the shearwave splitting parameter, and VS1 and VS2 are the fast and slow shearwave
velocities.
In the ultrasonic velocity measurements, the orientation of the core samples in the
borehole was not known and no factors that may create anisotropy, such as layering or fractures,
were visible. The core samples were oriented arbitrarily for Swave testing, so that the measured
shearwave splitting does not necessarily correspond to the true shearwave splitting for the
principal polarization directions. As shown in Figures 3.19 and 3.20, splitting was observed in
the core samples. 41 Saturated Vuggy Sample PC 1411 4698, Confining Pressure = 34.5 MPa Shear Splitting Parameter, γ, % 2.0 1.5 Brine: 50k ppm NaCl
Oil + 39% CO2
CO2 1.0 0.5 0.0
10 15 20 25 30 Pore Pressure, MPa Figure 3.19 Variation of shear splitting parameter with pore pressure at constant confining
pressure for saturated Vuggy sample PC 1411 4698. Dry Marly Samples PC 1213 4626 and PC 1411 4659.5 Shear Splitting Parameter, γ, % 7
6
5
4
3
2
PC 1213 4626
PC 1411 4659.5 1
0
0 5 10 15 20 25 30 Confining = Differential Pressure, MPa Figure 3.20 Variation of shear splitting parameter with differential pressure for dry Marly
samples PC 1213 4626 and 1411 4659.5. 35 42 Figure 3.19 shows that shearwave splitting increases with increasing pore pressure at a
constant confining pressure. This behavior is interpreted as being caused by the opening up of
preferentiallyoriented microcracks in the sample at high pore pressure (low differential pressure).
Similar rock behavior is discussed in Steensma (1999) and Walsh (1965). Qualitatively, this
shows what should happen around a highpressure injection well in the reservoir. Because the
microcracks or reservoir fractures are preferentially oriented, shearwave velocity varies with
azimuth. As pore pressure increases, increased shearwave splitting should be observed in the
seismic data.
In the Vuggy test sample, shearwave splitting does not vary with fluid saturation in a
meaningful way, and hence cannot be used to discriminate between fluids. Changes in shearwave splitting should indicate changes in pore (or differential) pressures in the reservoir.
Figure 3.20 shows similar shearwave splitting for the Marly core samples. Because
shearwave splitting is more sensitive to pressure than fluid saturation, it can also be measured on
dry core samples. At high differential pressure, most of the microcracks within the samples are
closed and splitting is low. As differential pressure decreases, splitting increases due to the
opening up of preferentially aligned microcracks in the rock.
Because the measurement orientations do not necessarily correspond to the principal
shear directions, it is not possible to quantitatively correlate changes in shearwave splitting from
these lab data to changes observed in the seismic data. The lab data show that shearwave
splitting may be more dependent on pressure in the Marly unit than in the Vuggy unit. However,
the greater fracture density in the Vuggy zone (Section 5.2) suggests that the Vuggy zone may
exhibit greater shearwave splitting at the seismic scale. 3.6 Summary of PanCanadian Rock Testing Results
Core Laboratories Canada Ltd. carried out an acoustic velocity study of Weyburn core
samples for PanCanadian Resources (Core Laboratories, 1998). Samples of the ‘Marly Porous’,
‘Vuggy Porous’, and ‘Vuggy Tight’ were tested under varying saturation, confining pressures,
and temperatures. The results are generally compatible with the testing results of this research.
Compressional wave velocity was found to be sensitive to both differential pressure and fluid 43 saturation, whereas shearwave velocity was sensitive to differential pressure only. Compared to
the samples tested in this research, the seismic velocities of samples tested by Core Laboratories
differed according to porosity. If the porosity of the Core Laboratories samples was higher,
velocities were lower and vice versa. This emphasizes the limitations of characterizing a
heterogeneous reservoir by testing of only a few small core samples.
For the ‘Marly Porous’ the variability of Pwave velocity from air to brinesaturation
and differential pressure from 10 to 25 MPa is 5% to 12%, consistent with sample PCP 1411
4659.5 and PCP 1213 4626. The variation of Pwave velocity of the ‘Vuggy Porous’ for the
same conditions is 5% to 9%, higher than results from sample PCP 1411 4698. The Pwave
velocity of the ‘Vuggy Tight’ varies by less than 1%, consistent with results from sample PCP
1213 4643. The effects of temperature variation on the seismic properties were found to be
insignificant. 3.7 Summary
Ultrasonic measurements were performed on core samples from the Vuggy and Marly
reservoir units for two main purposes. First, from measurements of P and Swave velocities on
dry samples, the pressure variation of elastic moduli is defined. These dry rock moduli are used
in isotropic rock physics modeling (Chapter 4). They are representative of the dry background
rock matrix, to which the addition of fractures and fluids will be modeled.
Second, ultrasonic measurements were made on a Vuggy core sample under different
pressure and fluid saturation conditions. From a simulation of the CO2 flood, the magnitude of
changes in seismic velocities due to changes in pressure and fluid saturation are observed. From
data sets at both constant differential pressure and constant confining pressure, the pressure
effects on the rock skeleton can be separated from the fluid saturation effects.
Seismic attributes are identified that can be used in the interpretation of the timelapse
seismic data for changes in fluid pressure and fluid composition. Changes in VP/VS may be used
to monitor changes in fluid composition. Changes in shearwave splitting may indicate changes
in fluid pressure. 44 Seismic attributes such as VP/VS and shearwave splitting are not uniquely defined by
fluid composition and pressure. For example, VP/VS is also sensitive to mineralogy, porosity,
pore shape, and fracture density. Shearwave splitting is also sensitive to fracture density,
direction, and aspect ratio. In timelapse monitoring, using seismic difference volumes allows the
use of changes in these attributes for interpretation. Whereas VP/VS or shearwave splitting from
a single survey may not be uniquely interpretable, changes in these attributes should be
constrained to changes in fluid composition and pore pressure, because the other variables remain
constant in the reservoir over time.
The magnitude of the changes in seismic velocities for the Vuggy sample tested is 1% to
2% for fluid composition changes and approximately 1% for expected pore pressure changes.
These changes are small and may not be resolved by the timelapse seismic data. However, the
pressure effects for the Marly unit are greater and the fluid effects are also expected to be greater
because of its higher porosity and lower elastic moduli. The expected changes in the seismic data
may be larger once the effect of fractures is included in the modeling. The data collected by
ultrasonic measurements on cores will be used in rock physics modeling to determine the
sensitivity of the timelapse seismic data to changes in fluid composition and pore pressure in the
Weyburn reservoir. 45 CHAPTER 4
ISOTROPIC ROCK PHYSICS MODELING 4.1 Introduction
Rock physics modeling is used to quantitatively relate changes in seismic response to
changes in reservoir fluid saturation and pressure. Modeling is used to calculate the sensitivity of
the reservoir zones to fluid saturation and pressure effects (Chapters 4 and 5), obtain the
parameters necessary for modeling synthetic seismograms (Chapter 6), and predict differences in
the time lapse seismic data from the reservoir simulation (Chapter 8). An anisotropic rock
physics model is required to realistically model the reservoir changes over time. The model used
must balance simplicity with accuracy and available input data. This chapter focuses on the
development of an isotropic rock physics model to predict elastic moduli based on porosity and
pressure. It is used to describe the matrix material for the anisotropic rock physics model
(Chapter 5).
Isotropic rock physics modeling is done based on ultrasonic core data (Chapter 3) and
geophysical logs (Section 4.3). These data are combined to obtain a pressure and porositydependent, isotropic, dry rock model for the elastic constants for the Marly and Vuggy reservoir
zones. The sensitivity of this model to fluid and pressure changes is calculated for comparison
with the more detailed anisotropic rock model. 4.2 Gassmann’s Theory
Fluid substitutions for isotropic rock models are commonly done with Gassmann’s
equation, discussed in Section 2.1 and repeated here
2 K SAT = K DRY K
1 − DRY KM +
,
φ
1 − φ K DRY
+
−2
K FL K M
KM (2.2) 46 where KSAT is saturated rock bulk modulus, KDRY is dry rock bulk modulus, KM is mineral
modulus, φ is porosity, and KFL is fluid bulk modulus.
Gassmann’s equation may not be directly applicable to seismic modeling of Weyburn
Field because the reservoir is fractured and anisotropic. However, the Gassmann model is useful
for comparison with ultrasonic core measurements. Gassmann’s equation is the low frequency
limit (relaxed fluidrock state) for wave propagation in saturated media; ultrasonic measurements
may be in the high frequency (unrelaxed fluidrock state) or transition zone (Batzle et al., 1999).
Comparison is useful to test the validity of assumptions in Gassmann’s equation and to estimate
error due to fluid viscosity effects and uncertain input parameters.
Gassmann’s equation is also commonly used for initial modeling, i.e. fluid substitutions
on geophysical logs. The results can then be compared to anisotropic rock physics modeling to
evaluate the effects of fractures on the seismic properties of the reservoir. 4.2.1 Comparison of Gassmann’s Equation and Lab Results
Gassmann’s equation is used to model saturated rock bulk modulus. Combined with rock
density estimation, Pwave velocity can be computed for direct comparison with ultrasonic Pwave measurements on core samples. This is shown for Vuggy sample 1411 4698 in Figure 4.1.
For Gassmann’s equation, the KDRY vs. pressure relationship was estimated from dry rock
ultrasonic measurements (Fig. 3.11), KFL was calculated using STRAPP and FLAG 4.0, porosity
was measured on the sample, and KM was estimated as 72 GPa (Reasnor, 2001). For Pwave
velocity (Eq. 2.1), variation of shear and bulk moduli with pressure and dry rock density were
measured on the core sample and fluid density values were calculated with STRAPP and FLAG
4.0.
As shown in Figure 4.1, Pwave velocity estimated with Gassmann’s equation is
consistently lower than the measured values by approximately 2.5% to 4%. Part of this
difference could be due to the rockfluid system being in the unrelaxed (high frequency) state in
the ultrasonic measurements. More importantly for timelapse seismic monitoring purposes, the
magnitude of the changes in Pwave velocity with different fluid saturation and pressure are
similar for the laboratory measurements as for Gassmann’s equation. 47 Saturated Vuggy Sample PC 1411 4698, Confining Pressure = 34.5 MPa
5100 PWave Velocity, m/s 5000 4900 4800
Brine: 50k ppm NaCl
Oil + 39% CO2
CO2 4700
10 Gassmann Brine
Gassmann Oil + 39% CO 2
Gassmann CO 2 15 20 25 30 Pore Pressure, MPa Figure 4.1 Comparison of Pwave velocities calculated with Gassmann’s equation and ultrasonic
core measurements. Confining pressure was fixed at 34.5 MPa. With increasing pore pressure,
differential pressure decreases. Gassmann’s equation predicts a greater fluid effect (3% change between brine and CO2) than do
the laboratory measurements (2% change between brine and CO2). This is probably because the
core had residual water and oil saturation during testing with CO2 as the primary pore fluid, and
modeling was done assuming pure phase fluids.
Modeling with Gassmann’s equation predicts a slightly greater pore pressure effect, 1.5%
to 2% over the pressure range tested, versus 1% to 1.5% for laboratory measurements. This could
be due to hysteresis due to different stress paths for the dry and saturated rock ultrasonic
measurements. Several other assumptions in Gassmann’s equation, such as homogeneity and 48 isotropy, are not strictly true for the core testing and probably contribute to the discrepancy
between modeling and laboratory results.
The laboratory measurements and comparable modeling with Gassmann’s equation are
extreme cases. In the reservoir, there will be mixed fluid saturation, not pure phase fluids. The
change in pore pressure due to EOR operations is also expected to be much less than the 13 MPa
range presented in Figure 4.1. However, the salinity of reservoir brine is 85,000 ppm and has a
higher bulk modulus than the 50,000 ppm brine used in lab testing (Section 2.2.1). The bulk
modulus and density of CO2 is also considerably lower at reservoir temperatures (63 °C) than at
room temperature. These effects will cause a greater difference in seismic properties at the endmember fluid saturations at reservoir conditions compared to lab testing. Good correspondence
between the limited lab results and modeling with Gassmann’s equation shows that it will also be
accurate in modeling differences at reservoir conditions. 4.2.2 Modeling with Gassmann’s Equation and Lab Data
Gassmann’s equation can be used to estimate the values of seismic attributes at reservoir
conditions (63 °C) so that the magnitude of differences due to fluid composition and pressure
effects can be estimated. As discussed in Section 2.2.1, the average value of brine salinity is
85,000 ppm, which will be used throughout the modeling in this chapter. The equivalent
isotropic confining pressure, PC, is calculated by PC = (σ V + σ H 1 + σ H 2 )
,
3 (4.1) where σV and σH1,2 are the vertical and orthogonal horizontal total stresses, respectively. The
total vertical stress due to the overburden material is estimated at approximately 1 psi/ft (22.6
MPa/km) depth. The total horizontal stresses or horizontal confining pressures are dependent on
material and deformation history. Horizontal stresses may range from values corresponding to
zero effective stress in areas with open vertical fractures to magnitudes greater than the vertical
stress.
There are no data on the magnitude of horizontal stresses at Weyburn Field. If an
effective Poisson’s ratio, ν, were known for the reservoir zone, the ratio of average horizontal to
vertical stress could be calculated by 49 σH
ν
,
=
σV 1 −ν (4.2) assuming the reservoir is elastic, isotropic, with no lateral strain. A common assumption when
data on horizontal stresses are lacking is that the horizontal stresses are on average one half of the
vertical stress. This corresponds to ν = 0.33. Average Poisson’s ratio calculated from the dipole
sonic log is 0.3, so this assumption is reasonable for Weyburn reservoir. For average horizontal
stresses equal to one half of vertical stress, the equivalent isotropic confining pressure for
Weyburn reservoir would be approximately 22 MPa. For an average of 15 MPa pore pressure
(Section 2.2), the average differential pressure in the reservoir is approximately 7 MPa. For a
range in pore pressure from 8 to 20 MPa, the range in differential pressure is 2 to 14 MPa.
Ultrasonic testing on saturated core samples (Section 3.5.1) showed little difference in the
sensitivity of the rock moduli to fluid and pressure changes at differential pressures of 6.9 MPa
and 13.8 MPa. However at very low differential pressures (i.e. 0 to 3 MPa), the sensitivity
differences may be significant.
Figures 4.2 to 4.5 show the variation of calculated Pwave velocity and VP/VS with fluid
composition and pressure for the reservoir temperature and pressure range. For the Vuggy zone,
the properties are based on sample PC 1411 4698 and for the Marly zone, they are based on
sample PC 1213 4626. Grain density values are 2.88 for the Marly zone and 2.70 for the Vuggy
zone, as measured on the samples (Tables 3.1 and 3.2). The variation of VP and VP/VS with pore
pressure is generally smooth. The kinks in the calculated curves for oil + CO2 mixtures are due to
phase changes in the fluid.
In general, brinesaturated rock has the highest Pwave velocity and VP/VS, followed by
rock saturated by oil, and then CO2. The Pwave velocities for mixtures of oil and CO2 are
complicated due to phase changes with changing pore pressure. Below the bubble point pressure,
a rock saturated with an oil + CO2 mixture may have a lower Pwave velocity than one saturated
with CO2. This is because the fluid bulk moduli are comparable but the density of the oil + CO2
mixture is higher than for pure CO2. However, below the modeled bubble point pressure of oil +
CO2 mixtures, the new phase coming out of solution may be a liquidlike fluid of CO2 enriched
with light hydrocarbons, not a gas, so the modeled results may be unreliable below the bubble
point pressures. 50 Marly Sample PCP 1213 4626, φ = 0.29, PC = 22 MPa
3800 PWave Velocity, m/s 3700
3600
3500
3400
3300
Brine: 85k ppm NaCl
Oil: API 29 GOR 30 L/L
CO2 3200 Oil + 25% CO2
Oil + 40% CO2
Oil + 66% CO2 3100
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 4.2 Variation of Marly dolostone Pwave velocity with fluid composition and pressure,
calculated with Gassmann’s equation. Marly Sample PCP 1213 4626, φ = 0.29, PC = 22 MPa
2.10 Brine: 85k ppm NaCl
Oil: API 29 GOR 30 L/L
CO2 2.05 Oil + 25% CO 2
Oil + 40% CO 2
Oil + 66% CO 2 VP/VS 2.00
1.95
1.90
1.85
1.80
8 10 12 14 16 18 Pore Pressure, MPa Figure 4.3 Variation of Marly dolostone VP/VS with fluid composition and pressure, calculated
with Gassmann’s equation. 20 51 Vuggy Sample PCP 1411 4698, φ = 0.11, PC = 22 MPa PWave Velocity, m/s 5000 4900 4800
Brine: 85k ppm NaCl
Oil: API 29 GOR 30 L/L
CO2 4700
8 Oil + 25% CO2
Oil + 40% CO2
Oil + 66% CO2 10 12 14 16 18 20 Pore Pressure, MPa Figure 4.4 Variation of Vuggy limestone Pwave velocity with fluid composition and pressure,
calculated with Gassmann’s equation. Vuggy Sample PCP 1411 4698, φ = 0.11, PC = 22 MPa
1.86
Brine: 85k ppm NaCl
Oil: API 29 GOR 30 L/L VP/VS Oil + 25% CO 2
Oil + 40% CO 2 CO2 1.84 Oil + 66% CO 2 1.82
1.80
1.78
1.76
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 4.5 Variation of Vuggy limestone VP/VS with fluid composition and pressure, calculated
with Gassmann’s equation. 52 Because VP/VS normalizes the effects of changing fluid density with pore pressure, it
shows separation based on fluid composition (bulk modulus). Because it also partially
normalizes the effect of changing shear and bulk moduli with differential pressure, it shows less
pressure dependence than Pwave velocity does.
Table 4.1 shows the sensitivity calculations made assuming a confining pressure of 22
MPa and a pore pressure of 15 MPa. In this model, the seismic velocities in the Marly zone are
roughly twice as sensitive to changes in fluid composition than the Vuggy zone is. The Marly is
more sensitive because it has lower elastic moduli and higher porosity than the Vuggy. The
seismic properties of the Marly zone are also more sensitive to pressure than are those of the
Vuggy zone. For both the Marly and the Vuggy zones, VP/VS is more sensitive to changes in
fluid composition than to changes in pore pressure. Table 4.1 Comparison of sensitivity of seismic attributes to changes in fluid composition and
pressure for Marly and Vuggy zones. Based on isotropic rock physics model from ultrasonic core
testing.
Change Marly Zone, φ = 0.29 Vuggy Zone, φ = 0.13 ∆VP ∆VS ∆VP/VS ∆VP ∆VS ∆VP/VS Brine to Oil 4.2 % 1.4 % 5.5 % 2.0 % 0.6 % 2.6 % Brine to CO2 7.8 % 3.2 % 10.7 % 3.8 % 1.3 % 5.1 % 0.60 % 0.90 % 0.30 % 0.19 % 0.15 % 0.04% Increase Pore
Pressure 1 MPa
(Oil @ 15 MPa) 4.3 Fluid Substitution on Geophysical Logs
Through ultrasonic lab testing, the variation of rock properties with pressure and fluid
saturation can be determined for a few samples with fixed lithology and porosity. With sonic
logging, the variation of seismic properties with lithology, porosity, and fluid saturation can be
measured, at approximately fixed pressure. An advantage of geophysical logging is that it 53 samples a larger volume of the reservoir than laboratory measurements on cores do. Geophysical
logs are affected to some extent by reservoir fractures. From geophysical logs, fluid substitutions
can be performed using Gassmann’s equation and the magnitude of the expected changes
compared to predictions from laboratory data.
Most of the necessary information for fluid substitutions using Gassmann’s equation can
be derived from geophysical logs. The input logs are resistivity, neutron porosity, bulk density,
Pwave velocity, and Swave velocity. Although it is outside the RCP study area, well 02042364
was chosen for modeling because it is the closest well with the required geophysical logs through
the reservoir interval. Figure 4.6 shows the location of well 02042364 in relation to the RCP
study area and other wells and core samples used in this research. 02042364 N
RCP
STUDY AREA Samples
1213
Samples
1411 410416613 1 km Figure 4.6 Location of wells 02042364 and 41041663 in relation to RCP study area. The
locations of core samples are also shown. The logs are shown in Figure 4.7. The high reservoir quality Marly zone is from
approximately 1400 m to 1405 m, with a transition zone upward to the Midale Evaporite at 1396 54 m. For sensitivity calculations, the Marly is assumed to extend from 1398 m to 1405 m. The
Vuggy zone is from 1405 to 1423 m. From the resistivity and porosity logs, water saturation, Sw,
can be estimated using Archie’s equation aR
Sw = m w
φ R
t 1
n , (4.3) where a is a tortuosity factor, Rw is the resistivity of formation water at formation temperature, Rt
is the true resistivity formation (corrected for invasion), m is the cementation exponent, and n is
the saturation exponent (Archie, 1942). The factors are usually calibrated to the area based on
core measurements. Using water saturation, average Weyburn brine (85,000 ppm NaCl) and oil
(API 29, GOR 30 L/L) properties, and FLAG 4.0 software, the bulk modulus (KFL) and density
(ρFL) of the pore fluids are calculated. A pore pressure of 15 MPa is assumed. 1395 Marly 1400 1410 Vuggy Depth, m 1405 1415 1420 1425
0 70 0.0 Resistivity, Ωm 0.2 Porosity 0.4 2.0 2.5 3.0 4000 6000 2000 3000 Density, g/cc PVelocity, m/s SVelocity, m/s Figure 4.7 Geophysical logs from well 02042364. Marly zone is from approximately 1398 m to
1405 m, Vuggy zone is from approximately 1405 to 1423 m. 55 From P and Swave velocity and bulk density, saturated rock bulk modulus is
determined, using Eq. 3.5 4 K = ρ VP2 − VS2 .
3 (3.5) Shear modulus is determined from the shear wave log and bulk density log, using Eq. 3.7,
repeated here µ = ρVS2 . (3.7) The two Swave velocities from the dipole sonic log were averaged for the calculations. The dry
rock bulk modulus can be calculated with an alternative form of Gassmann’s equation K DRY φK K SAT M + 1 − φ − K M
K FL .
=
φK M K SAT
+
−1−φ
K FL
KM (4.4) Because of the mixed lithology, mineral moduli are estimated by extrapolating the dry rock
moduli versus porosity relationships to zero porosity. For the Vuggy zone, this yields KM of 72
GPa and µM of 33.5 GPa; for the Marly zone KM is 83GPa and µM is inconclusive (Reasnor,
2001). For the Marly zone, µM is estimated to be 48 GPa, based on a range of values for dolomite
provided in Mavko et al. (1998). The values of the mineral moduli are used later as the maximum
bounds for the dry rock moduli (Section 4.4.2). In Eq. 4.4, porosity is taken from the neutron
porosity log, and KFL is from FLAG modeling of fluid saturation (Sw and So). The calculated
KDRY log is shown in Figure 4.8. 56 1395
1400
Marly Dolostone Depth, m 1405
1410
Vuggy Limestone 1415
1420
1425
0 10 20 30 40 50 60 KDRY, GPa Figure 4.8 Dry rock bulk modulus versus depth for well 02042364. Marly zone is approximately
1398 m to 1405 m; Vuggy zone is approximately 1405 m to 1423 m. For the Marly zone, KDRY ranges from 11 GPa to 22 GPa. Most of the variation is due to
the transition zone with interbedded anhydrite and dolomite. In the Vuggy zone, KDRY varies
from 20 GPa to 59 GPa. Why is KDRY so variable within the Vuggy zone? Possible factors
include mineralogy, fabric, porosity, fracturing. Figure 4.9 shows that porosity influences KDRY
but is not the only factor.
Although the geophysical logs are generally compatible with laboratory measurements, it
may be misleading to compare these values directly to core data because of the effect of reservoir
fractures, heterogeneities, and differing local stress distribution on the geophysical logs. A sonic
log often cannot accurately measure velocity in a fractured, vuggy limestone. Some energy goes
around the pore and fracture network rather than through it, dependent on scale. For a sonic log
frequency of 10,000 Hz and Pwave velocity of 5000 m/s, the wavelengths are 0.5 m.
Heterogeneities close to this size or larger will not be measured accurately. 57 Vuggy Zone, 1405 to 1423 m, Well 02042364
60 KDRY, GPa 50 40 30 20
0.00 0.05 0.10 0.15 0.20 Porosity Figure 4.9 Dry rock bulk modulus versus porosity for Vuggy Zone, well 02042364. The average values from sonic logs and ultrasonic testing are presented in Table 4.2. For
the sonic logs, the Backus average was used for the moduli (Backus, 1962), which for the vertical
incidence case corresponds to the (1/moduli) average. For the ultrasonic measurements, a
differential pressure of 7 MPa was used. Table 4.2 Comparison of average values of KDRY, µDRY from borehole sonic log and ultrasonic
core measurements.
Marly Zone
Moduli
(GPa) Vuggy Zone Sonic Log Ultrasonic Sonic Log Ultrasonic (well 02042364) (PC 1213 4626)
φ = 0.29 (well 02042364) (PC 1411 4698)
φ = 0.13 φ = 0.29 φ = 0.10 KDRY 15.0 14.5 33.9 31.2 µDRY 10.4 7.4 19.5 17.7 58 Except for µDRY, the moduli for the Marly zone calculated from core testing and the sonic
log compare well. For the Vuggy zone, the moduli from core testing are lower than those from
the sonic log, which is expected because the core sample has a higher porosity than the average
porosity in the geophysical log.
For fluid substitutions, Gassmann’s equation (Eq. 2.2) is used to calculate saturated rock
bulk modulus for a variety of fluid compositions. To calculate seismic velocities, saturated rock
density is required (Eq. 2.1, 2.3). This is estimated by using the bulk density log, porosity log,
and fluid density information derived from fluid saturation to calculate dry rock density with an
alternative form of Eq. 3.2 ρ sat = ρ dry + φρ fl . (4.5) Dry rock density is combined with porosity and new fluid density to get new saturated rock
density for velocity calculations.
Figure 4.10 shows the calculated Pwave velocity for the reservoir zone with brine, oil,
and CO2 saturation. Figure 4.11 presents the percent difference in Pwave velocity for oil and
CO2 saturation versus brine saturation. The variation of VP/VS with fluid saturation is shown in
Figure 4.12. Figure 4.13 shows the percent change in VP/VS from brinesaturated conditions for
oil and CO2 saturation. Fluid substitutions for density are shown in Figure 4.14. Figure 4.15
shows the percent change in density from brinesaturated conditions for oil and CO2 saturation.
The averaged or “blocked” variation of seismic impedance with fluid saturation for the Marly and
Vuggy zones is presented in Figures 4.16 and 4.17. 59 1395
Brine
Oil
CO2 1400
Marly Dolostone Depth, m 1405
1410
1415 Vuggy Limestone 1420
1425
3000 3500 4000 4500 5000 5500 6000 PWave Velocity, m/s Figure 4.10 VP versus depth, as a function of fluid saturation for well 02042364. Marly zone is
1398 m to 1405 m; Vuggy zone is 1405 to 1423 m. 1400
Marly Dolostone Depth, m 1405 1410
Vuggy Limestone 1415 1420
10 Oil
CO2 8 6 4 2 0 % Change in PWave Velocity from Brine Saturated Conditions Figure 4.11 Percent change in VP from brine saturation, for well 02042364. Marly zone is 1398
m to 1405 m; Vuggy zone is 1405 to 1423 m. 60 1400
Marly Dolostone Depth, m 1405 1410
Vuggy Limestone 1415
Brine
Oil
CO2 1420
1.6 1.8 2.0 2.2 VP/VS Figure 4.12 VP/VS versus depth, as a function of fluid saturation, for well 02042364. Marly zone
is 1398 m to 1405 m; Vuggy zone is 1405 to 1423 m. 1400
Marly Dolostone Depth, m 1405 1410
Vuggy Limestone 1415
Oil
CO2 1420
12 10 8 6 4 2 % Change in VP/VS from Brine Saturated Conditions Figure 4.13 Percent change in VP/VS from brine saturation, for well 02042364. Marly zone is
1398 m to 1405 m; Vuggy zone is 1405 to 1423 m. 0 61 1395
1400
Marly Dolostone Depth, m 1405
1410
1415
1420
1425
2.1 Vuggy Limestone
Brine
Oil
CO2 2.2 2.3 2.4 2.5 2.6 2.7 Density, g/cc Figure 4.14 Density versus depth, as a function of fluid saturation, for well 02042364. Marly
zone is 1398 m to 1405 m; Vuggy zone is 1405 to 1423 m. 1400
Marly Dolostone Depth, m 1405 1410
Vuggy Limestone 1415 1420
8 Oil
CO2 6 4 2 0 % Change in Density from Brine Saturated Conditions Figure 4.15 Percent change in density from brine saturation, for well 02042364. Marly zone is
1398 m to 1405 m; Vuggy zone is 1405 to 1423 m. 62 1400
Marly Dolostone, φ = 0.29 Depth, m 1405 1410
Vuggy Limestone, , φ = 0.10
Brine
Oil
Oil + 40% CO2
CO2 1415 1420
7 8 9 10 11 12
6 13 14 2 PWave Impedance, 10 kg/m /s Figure 4.16 Blocked Pwave impedance versus depth, as a function of fluid saturation, for well
02042364. 1400 Marly Dolostone, φ = 0.29 Depth, m 1405 1410 1415 1420
4.5 Vuggy Limestone, , φ = 0.10 Brine
Oil
Oil + 40% CO2
CO2 5.0 5.5 6.0 6.5
6 7.0 7.5 2 SWave Impedance, 10 kg/m /s Figure 4.17 Blocked Swave impedance versus depth, as a function of fluid saturation, for well
02042364. 63 In the Vuggy zone, sensitivity to fluid saturation varies considerably, i.e. from 0.5% to
8% between brine and CO2saturation. The variation in porosity and elastic moduli lead to the
differences in fluid sensitivity. The average sensitivity to fluids is shown in Table 4.3. While
100% saturation of any pore fluid is not realistic for the reservoir, these curves are useful for
showing the maximum sensitivity of seismic properties to fluid composition. In this model, the
seismic velocities in the Marly zone are roughly 50% more sensitive to changes in fluid
composition than the Vuggy zone is. Table 4.3 Comparison of sensitivity of seismic attributes to changes in fluid composition and
pressure for Marly and Vuggy zones. Based on geophysical logs for well 02042364.
Change Marly Zone, φ = 0.29 Vuggy Zone, φ = 0.10 ∆VP ∆VS ∆VP/VS ∆VP ∆VS ∆VP/VS Brine to Oil 3.3% 1.3% 4.6% 2.0% 0.4% 2.4% Brine to CO2 6.2% 2.8% 9.0% 3.9% 0.9% 4.8% 4.4 Integration of Ultrasonic Core and Geophysical Log Data
As shown in Table 4.2, the data from the ultrasonic core measurements and geophysical
logging are generally compatible. The core data sample a much smaller part of the reservoir than
the logging data do. As shown in Figure 4.8, there is significant variation in the seismic
properties of the Marly and Vuggy zones. The geophysical logs can be averaged so they are
representative of entire reservoir zones. While the absolute values of the ultrasonic core data may
not be representative of the reservoir zones due to sampling limitations, they can be used to
establish relative changes and pressure dependence of seismic parameters in the zones. The
geophysical logs provide information on the relationship between elastic moduli and porosity. A
porosity and pressuredependent model is achieved by integrating the ultrasonic core and the
geophysical log data. 64 4.4.1 Pressure Dependence of Elastic Moduli
The pressure dependence curves from ultrasonic testing are scaled so that the moduli
values at the estimated average differential pressure of 7 MPa are equal to the average moduli
calculated from geophysical logs. This combines the advantages of each data set. In calculating
the dry rock moduli from the geophysical logs, a pore pressure of 15 MPa and a differential
pressure of 7 MPa are assumed. For the Marly zone, the pressure dependence is described by
sample PC 1213 4626 and for the Vuggy zone it is described by 1411 4698. The pressure
dependence curves are scaled so that the dry rock moduli from lab testing and geophysical logs
are equal at a differential pressure of 7 MPa.
The equations for the variation of dry rock moduli for the total Marly zone are K DRY = 1.731x10 −4 P 3 − 1.325x10 −2 P 2 + 0.3542 P + 13.11 (4.6a) µ DRY = 1.157 x10 −4 P 3 − 9.828 x10 −3 P 2 + 0.2989 P + 8.78 , (4.6b) where the units of KDRY and µDRY are GPa and the units of pressure, P, are MPa. The pressuredependent moduli equations for the total Vuggy zone are K DRY = −2.560 x10 −3 P 2 + 0.2616 P + 32.23 (4.7a) µ DRY = −8.437 x10 −4 P 2 + 6.777 x10 −2 P + 19.05 . (4.7b) The dry rock moduli relations in Eqs. 4.6 and 4.7 are appropriate for the reservoir zones
in well 02042364, which is outside the RCP study area. They correspond to the average porosity
values of 0.29 for the Marly zone and 0.10 for the Vuggy zone. If these equations are to be used
for Weyburn reservoir rock with significantly different porosity, they must be adjusted so the
moduli and porosity correspond. Table 4.4 shows the sensitivity of seismic attributes to fluid
saturation and pressure. For velocity calculations, grain density values of 2.87 for the Marly zone
and 2.71 for the Vuggy zone were used. As expected, the sensitivity to fluid saturation is similar
to that of the geophysical logs only (Table 4.3), with the Marly zone being roughly 50% more
sensitive to changes in fluid composition than the Vuggy zone. 65 Table 4.4 Comparison of sensitivity of seismic attributes to changes in fluid composition and
pressure for Marly and Vuggy zones. Based on isotropic rock physics model from ultrasonic core
testing and geophysical logs. Porosity values are based on well 02042364.
Marly Zone, φ = 0.29 Change Vuggy Zone, φ = 0.10 ∆VP ∆VS ∆VP/VS ∆VP ∆VS ∆VP/VS Brine to Oil 3.4 % 1.4 % 4.7% 2.2 % 0.4 % 2.6 % Brine to CO2 6.2 % 3.2 % 9.1 % 4.2 % 1.0 % 5.2 % Increase Pore
Pressure 1 MPa
(Oil @ 15 MPa) 0.64 % 0.90 % 0.26 % 0.18 % 0.15 % 0.04 % 4.4.2 Porosity Dependence of Elastic Moduli
As Figure 4.9 shows, dry rock moduli are related to porosity. With Gassmann’s equation
(Eq. 2.2), the chosen values of porosity must correspond to KDRY for realistic sensitivity of KSAT
to fluids. As an example, if KDRY is held constant and porosity is reduced by one half, the rock
becomes roughly twice as sensitive to fluids. This seems counterintuitive, since rocks with higher
porosity are typically more sensitive to fluids.
KDRY is lower than KM because of the compressible pore space, as shown in Eq. 4.8: 1 K DRY = 1
φ ∂v p
+
K M v p ∂σ dry , (4.8) where KDRY is dry rock bulk modulus, KM is mineral bulk modulus, φ is porosity, vp is pore
volume, and ∂v p
∂σ dry is the derivative of pore volume with respect to hydrostatic confining stress under dry conditions (Mavko et al., 1998). Pore space stiffness is defined as v p ∂σ
∂v p dry . If a small value of porosity causes a large reduction of KDRY from KM, the pore space stiffness is low,
making it very sensitive to fluids. It is equivalent to assuming the pores have a very small aspect
ratio in crack models (Hudson, 1980, 1981; Kuster and Toksoz, 1974). 66 In the reservoir simulation model (Section 8.2.1), porosity varies by a factor of about four
in each layer. A porositydependent model for elastic moduli is therefore necessary to
realistically perform fluid substitutions before upscaling the layers to the reservoir zones. This
porositydependent model is achieved by incorporating a porosity scaling factor obtained from
the geophysical logs. Figures 4.18 though 4.21 show the variation of KDRY and µDRY with neutron
porosity for the Marly and Vuggy zones, with linear trend lines fit to the data. For the Marly
zone in well 02042364, values were taken from a depth of 1396 m to 1405 m, so that a wide
range of porosity was included.
The variability in the relations between dry rock moduli and porosity shows that Eq. (4.8)
only approximates Weyburn reservoir rock behavior. This is probably due to changes in
lithology within the reservoir zones, so that KM and ∂v p
∂σ dry vary. The trend lines could be changed to polynomials for a better fit to the data for well 02042364. However the uncertainty
and variability in the reservoirwide relations probably do not warrant this refinement. Care
should be taken not to extrapolate the trend lines for the Marly zone moduli to low porosity
values where unreasonably high or negative moduli would be predicted. The maximum dry rock
moduli should be less than or equal to the mineral moduli.
The porosity scaling factors are derived from the linear trend lines fit to 1/moduli versus
porosity. The porosityscaling factors for the Marly zone are K DRY (φ )
1
=
K DRY (φ = 0.29) 3.50φ − 0.015 (4.9a) 1
µ DRY (φ )
=
.
µ DRY (φ = 0.29) 3.17φ + 0.0807 (4.9b) The porosityscaling factors for the Vuggy zone are K DRY (φ )
1
=
K DRY (φ = 0.1) 5.60φ + 0.44 (4.10a) 1
µ DRY (φ )
=
.
µ DRY (φ = 0.1) 5.06φ + 0.494 (4.10b) 67 Marly Zone, 1396 to 1405 m, Well 02042364
0.10 0.08 1/Kdry = 0.233φ  0.00151
2 1/KDRY, 1/GPa R = 0.72
0.06 0.04 0.02 0.00
0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.35 0.40 Porosity Figure 4.18 Variation of 1/KDRY with porosity for the Marly zone. Marly Zone, 1396 to 1405 m, Well 02042364
0.14
0.12
1/µDRY, 1/GPa 1/µdry = 0.305φ + 0.00737
0.10 2 R = 0.76 0.08
0.06
0.04
0.02
0.00
0.05 0.10 0.15 0.20 0.25
Porosity Figure 4.19 Variation of 1/µDRY with porosity for the Marly zone. 0.30 68 Vuggy Zone, 1405 to 1423 m, Well 02042364
0.06
1/Kdry = 0.165φ + 0.013
2 1/KDRY, 1/GPa R = 0.46
0.04 0.02 0.00
0.00 0.05 0.10 0.15 0.20 0.25 0.20 0.25 Porosity Figure 4.20 Variation of 1/KDRY with porosity for the Vuggy zone. Vuggy Zone, 1405 to 1423 m, Well 02042364
0.10 0.08 1/µdry = 0.260φ + 0.0252
2 1/µDRY, 1/GPa R = 0.46
0.06 0.04 0.02 0.00
0.00 0.05 0.10 0.15
Porosity Figure 4.21 Variation of 1/µDRY with porosity for the Vuggy zone. 69 The porosity scaling factors are independent of pressure and pressure scaling is
independent of porosity. In reality, low porosity rocks may be less pressure dependent. These
assumptions are made for simplicity and data availability and should only cause minor deviations
over the pore pressure changes in the RCP study area. Equations (4.9) and (4.10) should be used
with caution outside of the porosity range from which they were developed. The values of dry
rock moduli should be less than the values of mineral moduli.
Table 4.5 shows the sensitivity of the seismic attributes to fluid saturation and pressure in
the porositydependent model. The main difference is that porosity of the Marly zone is reduced
to 0.24, the average value for the RCP study area. The sensitivity to fluids and pressure is almost
the same as in the model with porosity of 0.29 (Table 4.4), showing that porosity variation has
only a second order effect on fluid and pressure sensitivity at Weyburn Field. Seismic velocities
in the Marly zone are approximately 50% more sensitive to fluid changes than are the velocities
in the Vuggy zone. Table 4.5 Comparison of sensitivity of seismic attributes to changes in fluid composition and
pressure for Marly and Vuggy zones. Based on isotropic rock physics model from ultrasonic core
testing and geophysical logs. Porosity values are averages for RCP study area.
Change Marly Zone, φ = 0.24 Vuggy Zone, φ = 0.10 ∆VP ∆VS ∆VP/VS ∆VP ∆VS ∆VP/VS Brine to Oil 3.3% 1.1% 4.3% 2.2 % 0.4 % 2.6 % Brine to CO2 6.0% 2.5% 8.3% 4.2 % 1.0 % 5.2 % Increase Pore
Pressure 1 MPa
(Oil @ 15 MPa) 0.64% 0.90% 0.26% 0.18 % 0.15 % 0.04 % 4.5 Summary
A pressure and porositydependent isotropic model for elastic moduli of the reservoir
zones is developed from ultrasonic core data and geophysical logs. The pressure dependence is 70 taken from ultrasonic measurements, which are then scaled to average moduli for the reservoir
zones. The variation of elastic moduli with porosity is added from geophysical log data. This
model is used to calculate the sensitivity of Weyburn reservoir rock to changes in fluid saturation
and pressure using Gassmann’s equation, and as input to the anisotropic rock physics model
described in Chapter 5.
Gassmann’s equation does not accurately reproduce the seismic velocities measured on
core velocities, perhaps due to frequency effects. However it does reproduce the differences in
seismic velocities due to fluid and pressure changes, so it is acceptable for fluid substitution on
Weyburn reservoir rock for timelapse monitoring purposes. The sensitivity of models from core
data alone, geophysical logs alone, and integrated core and log data is approximately the same for
reasonable ranges in porosity for the reservoir zone. This shows that changes in seismic
velocities due to changes in fluid saturation and pressure in the reservoir can be calculated using
simplified models representative of the Marly and Vuggy zones. Spatial variability of lithology
within these zones causes only a second order effect on sensitivity of the reservoir to fluid
saturation and pressure. 71 CHAPTER 5
ANISOTROPIC ROCK PHYSICS MODELING 5.1 Introduction
The dry, isotropic rock physics model from the combined ultrasonic test data and
geophysical logs should give a reasonable estimate of the elastic moduli of the reservoir matrix
material. To more accurately simulate the seismic properties of the reservoir, the effect of
fractures should be included in an anisotropic rock physics model. Hudson’s crack model is used
to calculate the elastic stiffness tensor of the dry, fractured rock (Section 5.2). Brown and
Korringa’s equation is then applied to calculate the elastic stiffness tensor of the saturated,
fractured rock (Section 5.3). The elastic coefficients and density data are combined to compute
seismic velocities (Section 5.4). The end result is an anisotropic rock physics model that can be
used to model the variation of seismic velocities with pressure and fluid saturation. 5.2 Fractures
There are several indications of fractures in Weyburn reservoir. Oil production shows
definite trend and offtrend patterns. For this reason horizontal wells are oriented parallel to the
apparent dominant fracture direction. Shearwave splitting is observed in both the dipole sonic
logs and multicomponent seismic data. Reservoir core samples have fractures.
Core samples show varying degree of fracturing, with the Marly zone less fractured than
the Vuggy zone. There are vertical anhydritefilled fractures less than 1 mm thick, vertical
extensional fractures off stylolites, and large rough fractures 0.1 m to 0.5 m long. Most of the
large fractures do not have mineral growth along the fracture surfaces, making discrimination
between natural and drillinginduced fractures ambiguous. Thin section analysis would provide
more definitive classification of fracture origin. 72 5.2.1 Previous Fracture Characterization Studies
Bunge (2000) obtained information on fracture characteristics from core and image logs.
Fracture height, aperture, density, and orientation were quantified for three sets of open fractures.
Table 5.1 lists properties of the open fracture sets at Weyburn. Table 5.1 Description of fractures in Marly and Vuggy zones (from Bunge, 2000).
Fracture Set Strike Azimuth
(degrees) Fracture Density Marly
(fractures/m) Fracture Density Vuggy
(fractures/m) A 40° ± 5° 2.3 3.8 B 285° ± 7° 1.5 2.5 C 328° ± 11° 1.0 1.6 The average dip of all fracture sets is 80°, with fractures dipping from 80° to 90°
accounting for 69% of fractures (Bunge, 2000). This supports the assumption that the fractures
can be modeled as vertical. The average fracture aperture is 0.1 mm, with a range of 0.004 to 0.2
mm. Fracture height averaged 30 cm for all sets and for both the Marly and Vuggy zones, with
fractures 10 cm to 30 cm high accounting for 70% of fractures.
Fracture porosity may be calculated by φf = afdf , (5.1) where af is fracture aperture (length) and df is fracture density (1/length). This equation is for
linear porosity, which assumes that the fractures extend infinitely in the plane normal to the line.
For Vuggy zone set A, this gives an estimate of 0.038% average fracture porosity.
An alternative method is to assume that the linear fracture density is the same in all
directions, so that the volumetric fracture density in fractures/m3 is three times the linear fracture
density in fractures/m. The geometry of the fractures must be chosen. Assuming the fractures are
disks, with fracture length equal to fracture height, and width given by aperture, φf = πh 2 a f d f
f
4 , (5.2) 73 where hf is the fracture height. For Vuggy zone set A, this gives an estimate of 0.0081% for φf, a
factor of about five lower than the fracture porosity from Eq. (5.1). Without additional
information on fracture shape, it is not possible to know the correct fracture porosity.
Ideally we would either have data on fracture shape and porosity or have independent
seismic velocities of the fractured Weyburn reservoir with which we could calibrate Hudson’s
crack model. Because the Vuggy and Marly zones are thin, accurate seismic velocities over these
separate intervals are difficult to obtain.
For a single fracture set, crack density, ηc, is approximately equal to Thomsen’s shear
splitting parameter, γ, for transversely isotropic media. Shearwave splitting values for the
reservoir level are available from the dipole sonic log. These values are representative of seismic
data if the well penetrates representative fractures sets, which may be unlikely for a vertical well
sampling vertical fractures. For well log 02042364, the average shear splitting is 2.2% for the
Marly zone and 10% for the Vuggy zone. With a fracture aspect ratio of 0.00033 (0.1 mm
average aperture ÷ 30 cm average height) and density of fracture set A, this corresponds to a
fracture porosity of 0.0031% for the Marly zone and 0.014% for the Vuggy zone. For the Marly
and Vuggy zones, these are closer to the volumetrically derived fracture porosity (Eq. 5.2) than
the linearly derived porosity (Eq. 5.1). The core and image logs show approximately 40% less
fracture density for the Marly zone versus the Vuggy zone, but the dipole sonic logs show 75%
less fracturing in the Marly zone compared to the Vuggy zone. The ultrasonic testing of the
Marly core shows that there may be 3% to 5% shearwave splitting (Figure 3.20). 5.2.2 Hudson’s Model for Cracked Media
Hudson’s model for cracked media (1980, 1981) can be used to compute the effective
moduli of an isotropic background matrix with fractures. The model is based on scattering theory
for a wave field in an elastic solid with idealized ellipsoidal cracks. The cracks may be either dry,
fluidfilled, or filled with “weak” solids. Hudson’s model assumes that the cracks are isolated
with respect to fluid flow, which corresponds to highfrequency fluidrock behavior. Because
seismic data are probably in the low frequency range at Weyburn, Hudson’s model is used to
calculate the elastic constants for a dry, fractured rock and then Brown and Korringa’s equation (a
generalization of Gassmann’s equation) is used to calculate the saturated properties. 74 In Hudson’s crack model, fractures are described by the crack density, ηc. Crack density
is formally defined as ηc = 3φ f
4πα , (5.3) Where φf is fracture (crack) porosity and α is crack aspect ratio. There is inherent nonuniqueness in the interpretation of crack density. For example, a fracture set with high fracture
porosity and large aspect ratio could have the same crack density as a set with low fracture
porosity and small aspect ratio. 5.2.3 Single Fracture Set – Hudson’s Model
Because the main focus of this work is on characterizing pore pressure and fluid
saturation changes, not fractures, rock physics modeling is done only for the average parameters
for a single set of fractures. The following summary of Hudson’s (1980,1981) model for a single
set of fractures with normals parallel to the x1axis (vertical fractures) is adapted from Mavko et
al. (1998) and Bakulin et al. (2000). The effective medium exhibits HTI symmetry (transversely
isotropic with a horizontal symmetry axis). The effective moduli, cijeff are described by
eff
0
1
2
cij = cij + cij + cij , (5.4) where cij0, cij1, cij2 are the isotropic background moduli and the first and secondorder corrections
for cracks, respectively. For vertical wave propagation, the moduli of interest are c33 (Pwave),
c44 (S1wave), c66 (S2wave). The isotropic background moduli are determined by Lame’s
parameters λ and µ. The corrections are
1
c11 = − (λ + 2µ ) 2 µ 2
η cU 1 , c11 = q
( λ + 2 µ )(η cU 1 ) 2
15 (5.5a) 1
c13 = − q
λ (λ + 2 µ ) 2
2
η cU 1 , c13 = λ (η cU1 ) 2
15
µ (5.5b) 1
c33 = − q
λ2
λ2
2
η cU 1 , c33 =
(η cU 1 ) 2
µ
15 ( λ + 2 µ ) (5.5c) 75 2
c1 = 0 , c44 = 0
44 (5.5d) 2 µ (3λ + 8µ )
(η cU 3 ) 2
15 (λ + 2 µ ) (5.5e) U1 = 4( λ + 2 µ )
,
3( λ + µ ) (5.5f) U3 = 1
2
c66 = − µη cU 3 , c66 = 16( λ + 2 µ )
,
3(3λ + 4 µ ) (5.5g) where ηc is crack density and λ2
λ
q = 15 2 + 28 + 28 .
µ
µ (5.5h) The other stiffness constants are obtained using the structure of the stiffness tensor for HTI
symmetry with the x1symmetry axis, as shown in Eq. 5.6: c11
c 13
c13
cij = 0
0 0 c13
c33
c23 c13
c23
c33 0 0 0 0 0 0 0 0 c44 0 0 0 0 c66 0 0 0 0 0 0 , c23 = c33 − 2c44 .
0 0 c66 0 (5.6) Figure 5.1 shows the variation of c33 with crack density for the Marly and Vuggy zones. Figure
5.2 presents the variation of c44, and c66 with crack density for the Marly and Vuggy zones. In the
Marly and Vuggy zones, the crack density is estimated to vary from 0% to 10%. Seismic
velocities are proportional to the square root of elastic coefficients. For vertically traveling waves
in an HTI medium, crack density has no effect on c44 (VS1), and approximately linear effects on
c33 (VP) and c66 (VS2). 76 Differential Pressure = 7 MPa
65
60
Dry Rock C33, GPa 55
50
45
40
35
30
25
20
0.00 Vuggy, φ = 0.1
Marly, φ = 0.24
0.01 0.02 0.03 0.04
0.05
0.06
Crack Density 0.07 0.08 0.09 0.10 Figure 5.1 Variation of c33 (Pwave modulus) with crack density at constant differential pressure,
calculated using Hudson’s crack model. Differential Pressure = 7 MPa
22
Dry Rock C44 and C66, GPa 20
18
16
14
12
10
8
6
0.00 C44 Vuggy, φ = 0.1
C66 Vuggy, φ = 0.1 0.01 0.02 0.03 C44 Marly, φ = 0.24
C66 Marly, φ = 0.24 0.04
0.05
0.06
Crack Density 0.07 0.08 0.09 0.10 Figure 5.2 Variation of c44 (S1wave modulus) and c66 (S2wave modulus) with crack density at
constant differential pressure, calculated using Hudson’s crack model. 77 Figure 5.3 shows the variation of c33 with differential pressure at constant crack density
for the Marly and Vuggy zones. Figure 5.4 presents the variation of c44, and c66 with differential
pressure at constant crack density for the Marly and Vuggy zones. For modeling purposes, the
crack density is 3% in the Marly zone and 10% in the Vuggy zone. In this constant crack density
model, the relative amount that cij is reduced due to cracks is largely independent of differential
pressure. In other words, cij0 varies but cij1 and cij2 in Eq. (5.4) are approximately constant. Crack Density: Marly 3%, Vuggy 10%
65
60
Dry Rock C33, GPa 55
50
45
40
35
30 Vuggy, φ = 0.1
Marly, φ = 0.24 25
20
2 4 6 8
10
Differential Pressure, MPa 12 14 Figure 5.3 Variation of c33 (Pwave modulus) with differential pressure at constant crack density,
calculated using Hudson’s crack model. 78 Crack Density: Marly 3%, Vuggy 10%
24
Dry Rock C44 and C66, GPa 22
20
18
16
14
12
C44
C66
C44
C66 10
8
6 Vuggy, φ = 0.1
Vuggy, φ = 0.1
Marly, φ = 0.24
Marly, φ = 0.24 4
2 4 6 8
10
Differential Pressure, MPa 12 14 Figure 5.4 Variation of c44 (S1wave modulus) and c66 (S2wave modulus) with differential
pressure at constant crack density, calculated using Hudson’s crack model. 5.2.4 Multiple Fracture Sets
For reservoir development and seismic response, only one set of vertical fractures is
usually considered at Weyburn Field. However, the work by Bunge (2000) shows that there may
be several fracture sets. Several fracture sets would reduce the elastic moduli further from one
set, change the symmetry of the elastic stiffness matrix, and possibly decrease shearwave
splitting. While these effects are important, it is not practical to model them for several reasons.
First, because of the lack of information on fracture dimensions, the uncertainty in
fracture porosity would lead to unacceptable error in crack density and the elastic stiffness matrix.
Although Hudson’s crack model can be used to model the effect of multiple fracture sets, several
of its assumptions would be grossly violated in the case of Weyburn Field. Second, the cracks
are assumed to be uniform ellipses (pennyshaped), which is convenient but not necessarily
realistic. Third, Hudson’s model should only be used for crack density less than 0.1. A crack
density of 0.1 for a single fracture set is necessary to match the average shear splitting in the
Vuggy zone; multiple fracture sets would require a higher crack density to match the shearwave 79 splitting. Fourth, multiple intersecting fracture sets would interact. The limitation on crack
density is imposed because the fractures are isolated in Hudson’s model. There are other crack
models besides Hudson’s, such as those by Kuster and Toksoz (1974), O’Connell and Budianski
(1974), Cheng (1993), but they impose similar geometric restrictions on cracks and are usually
valid only for low crack density. Fifth, the fracture characteristics change with changes in
differential stress. Existing theories to model this effect (Eshelby, 1957; Walsh, 1965) require
information on fracture distributions that is not available. The pressure variation of crack density
could be modeled empirically or theoretically but there are no data with which to calibrate the
model.
How significant are the errors introduced by the simplified crack model used in this
research? For the interpretation of the zerooffset Pwave volumes, the differences are minor.
For the interpretation of Pwave AVO or multicomponent data, the differences between the
simplified crack model and the reservoir are important. The crack model could be improved in
terms of both static fracture parameters and dynamic fracture characterization (dependence of
fracture parameters on pressure). 5.3 Stiffness Tensor of Saturated Anisotropic Rock – Brown and
Korringa’s Equation
Brown and Korringa’s (1975) equation relates the effective moduli of a saturated,
anisotropic rock to a dry, anisotropic rock. Like Gassmann’s equation, it applies in the low
frequency seismic limit. The formulation given here does not include a separate pore space
compressibility, and assumes a single elastic compliance tensor for the mineral(s) and complete
saturation.
( sat
( dry
S ijkl ) = S ijkl ) − ( dry
M
(
M
( S ijαα ) − S ijαα )( S kldry ) − S klαα )
αα
( dry
M
( Sααββ) − Sααββ ) + ( β fl − β M )φ , (5.7) where Sijkl(sat), Sijkl(dry), SijklM are the effective elastic compliance tensors of the saturated rock, dry
rock, and mineral material, respectively; β fl and βM are compressibility of the pore fluid and
mineral material; φ is porosity. The elastic compliance tensor is the matrix inverse of the elastic 80 stiffness tensor. The foursubscript notation is used (Aki and Richards, 1980) and summation is
over α and β. A more intuitive formulation of Eq. 5.7, which uses the 2subscript notation, is
given in Cardona et al. (2001). The symmetry of the saturated rock is the same as that of the dry
rock (Eq. 5.6).
For the case of HTI symmetry and vertically traveling waves, the moduli that control
shear wave velocity, c44 and c66, are not influenced by fluid saturation. C33, which controls
vertical incidence Pwave velocity, increases with increasing fluid bulk modulus, as shown in
Figure 5.5. Therefore fluid saturation should not change shearwave splitting but will affect
VP/VS ratio. As shown in Figure 5.5, the elastic properties of the Marly zone are more sensitive
to fluid saturation than are those of the Vuggy zone. Differential Pressure = 7 MPa, Crack Density: Marly 3%, Vuggy 10%
70 Saturated Rock C33, GPa 65
60
55
50
Vuggy, φ = 0.1
Marly, φ = 0.24 45
40
35
30 oil gas 25
0.0 0.5 1.0 1.5
2.0
Fluid Bulk Modulus, GPa brine
2.5 3.0 Figure 5.5 Modeled variation of c33 (Pwave modulus) with fluid bulk modulus at constant crack
density and differential pressure. 5.4 PressureDependent, Anisotropic, Saturated Rock Physics Model
The stiffness matrix of the saturated, anisotropic rock is combined with bulk density to
calculate the seismic velocities of the rock. Figures 5.6 and 5.7 show the variation of density with 81 pore pressure and fluid composition, for a constant confining pressure and crack density.
Because bulk density is linearly related to fluid density at constant porosity, matrix density, and
grain density (Eqs. 3.3 and 4.5), the bulk density relations are simply scaled versions of the fluid
density relations (Figure 2.11). Grain density values are 2.87 for the Marly zone and 2.71 for the
Vuggy zone.
Figures 5.8 through 5.13 present the variation of VP, VS1, and VS2 with pore pressure and
fluid composition at constant confining pressure and crack density, for the Marly and Vuggy
zones. The seismic velocities change due to variation of fluid density and bulk modulus with
pore pressure and variation of the elastic stiffness matrix (cij) with differential pressure. Marly Dolostone, φ = 0.24, ρgrain = 2.87
2.45 Bulk Density, g/cc 2.40 2.35 2.30
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO2
CO2 2.25 2.20
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 5.6 Marly bulk density as a function of fluid saturation and pore pressure. Fluid density
values are shown in Figure 2.11. 82 Vuggy Limestone, φ = 0.10, ρgrain = 2.71 Bulk Density, g/cc 2.54 2.52 2.50 2.48 Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO 2
CO2 2.46
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 5.7 Vuggy bulk density as a function of fluid saturation and pore pressure. Fluid density
values are shown in Figure 2.11. Marly Dolostone, φ = 0.24, Crack Density = 0.03, PC = 22 MPa
4300 PWave Velocity, m/s 4200
4100
4000
3900
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO2
CO2 3800
3700
3600
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 5.8 Marly Pwave velocity as a function of fluid saturation and pore pressure at constant
crack density and confining pressure. 83 Marly Dolostone, φ = 0.24, Crack Density = 0.03, PC = 22 MPa S1Wave Velocity, m/s 2500 2400 2300
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO2
CO2 2200 2100
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 5.9 Marly S1wave velocity as a function of fluid saturation and pore pressure at constant
crack density and confining pressure. Marly Dolostone, φ = 0.24, Crack Density = 0.03, PC = 22 MPa
2400 S2Wave Velocity, m/s 2350
2300
2250
2200
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO2
CO2 2150
2100
2050
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 5.10 Marly S2wave velocity as a function of fluid saturation and pore pressure at constant
crack density and confining pressure. 84 Vuggy Limestone, φ = 0.10, Crack Density = 0.10, PC = 22 MPa PWave Velocity, m/s 5100 5000 4900 4800 Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO2
CO2 4700
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 5.11 Vuggy Pwave velocity as a function of fluid saturation and pore pressure at constant
crack density and confining pressure. Vuggy Limestone, φ = 0.10, Crack Density = 0.10, PC = 22 MPa S1Wave Velocity, m/s 2840
2820
2800
2780
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO2
CO2 2760
2740
8 10 12 14 16 18 Pore Pressure, MPa Figure 5.12 Vuggy S1wave velocity as a function of fluid saturation and pore pressure at
constant crack density and confining pressure. 20 85 Vuggy Limestone, φ = 0.10, Crack Density = 0.10, PC = 22 MPa S2Wave Velocity, m/s 2540
2520
2500
2480
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO2
CO2 2460
2440
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 5.13 Vuggy S2wave velocity as a function of fluid saturation and pore pressure at
constant crack density and confining pressure. In this constant crack density model, VS2 is a scaled version of VS1. For the Marly zone,
the reduction factor is approximately 3.4% for 3% crack density. For the Vuggy zone, the
reduction factor is approximately 11.9% for 10% crack density. The effect of cracks is also to
reduce VP by approximately 0.4% for the Marly zone and 1.1% for the Vuggy zone compared to
the intact, isotropic model.
The variation of VP/VS with pore pressure and fluid composition is shown in Figures 5.14
and 5.15. VP/VS1 was chosen because it is usually more accurately derived from seismic data
than VP/VS2, due to noise considerations. In this constant crack density model, VP/VS2 is simply a
scaled version of VP/VS1 and contains no new information.
Figures 5.14 and 5.15 show that the seismic attribute VP/VS1 is less sensitive to pressure
than VP alone. Table 5.2 summarizes the sensitivity of these seismic attributes to changes in fluid
saturation and pressure, as calculated with the saturated, anisotropic rock physics model.
Calculations are made assuming a confining pressure of 22 MPa and pore pressure of 15 MPa. 86 Marly Dolostone, φ = 0.24, Crack Density = 0.03, PC = 22 MPa
1.90
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO 2
CO2 1.85 VP/VS1 1.80
1.75
1.70
1.65
1.60
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 5.14 Marly VP/VS1 as a function of fluid saturation and pore pressure at constant crack
density and confining pressure. Vuggy Limestone, φ = 0.10, Crack Density = 0.10, PC = 22 MPa
1.84
1.82 Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO2
CO2 VP/VS1 1.80
1.78
1.76
1.74
1.72
1.70
8 10 12 14 16 18 Pore Pressure, MPa Figure 5.15 Vuggy VP/VS1 as a function of fluid saturation and pore pressure at constant crack
density and confining pressure. 20 87 Table 5.2 Comparison of sensitivity of seismic attributes to changes in fluid composition and
pressure for Marly and Vuggy zones. Based on saturated, anisotropic rock physics model.
Porosity values are averages for RCP study area.
Change Marly Zone, φ = 0.24 Vuggy Zone, φ = 0.10 ∆VP ∆VS1 ∆VP/VS1 ∆VP ∆VS1 ∆VP/VS1 Brine to Oil 3.4% 1.1% 4.4% 2.8% 0.4% 3.2% Brine to CO2 6.3% 2.5% 8.6% 5.6% 1.0% 6.5% Increase Pore
Pressure 1 MPa
(Oil @ 15 MPa) 0.65% 0.90% 0.25% 0.16% 0.15% 0.01% Table 5.2 can be compared to Table 4.5, which shows the sensitivity of the saturated,
isotropic rock physics model. The fractures increase the sensitivity of the rock to fluid saturation
changes. After accounting for fractures, the Marly zone is approximately 25% more sensitive to
changes in fluid saturation than the Vuggy zone.
Compared with ∆VP, ∆VP/VS1 is more sensitive to fluid composition change and less
sensitive to change in pore pressure. Because the current model does not include the variation of
crack density with pressure, shear wave splitting cannot be modeled. The estimation of ∆VP/VS2
with pressure would also be only approximate, since VS2 would vary differently with pressure
than VS1 in a pressuredependent crack density model. If the multiple vertical fracture sets
present at Weyburn were modeled, both VS1 and VS2 would vary differently with pressure than in
the current model.
The sensitivity values are calculated for the average of the reservoir zones. In the high
quality reservoir zone in the Marly, the sensitivity to fluids is at least twice that of the average
values of the Vuggy zone (Figures 4.11, 4.13). The Marly core samples for ultrasonic testing
were taken from the high quality Marly zone and exhibit greater sensitivity to fluids (Table 4.1)
than the anisotropic model (Table 5.2). By including lower quality reservoir rock in the Marly
models, the sensitivity to fluids decreases but the height of the zone increases. This is done to
best match the Marly model in the reservoir simulator (Chapter 8). 88 Although the sensitivity of the reservoir rock is shown for pure fluid phases in Tables 4.1,
4.3, 4.4, 4.5, and 5.2, this maximum sensitivity is not much different than the maximum changes
expected in the reservoir. This is because the effective bulk modulus of a multiple phase fluid
mixture with a highly compressible phase (i.e. gas) is similar to the bulk modulus of the
compressible phase. At Weyburn reservoir, the acoustic properties of a fluid mixture with a light
hydrocarbonenriched CO2 phase are similar to those of pure CO2. The maximum expected
difference in the reservoir is from 2phase (brineoil) saturation to 3phase (brineoil saturated
with CO2CO2 with light hydrocarbons) saturation. The change in seismic velocities should be
between the values for brine to oil and brine to CO2 given in Table 5.2, for Pwave velocity, the
expected changes are approximately 4% to 5% for the reservoir. Synthetic seismic modeling
(Chapter 6) is used to estimate the effect of these changes on the timelapse seismic reflection
data and gauge their detectability. 5.5 Summary
An anisotropic, pressuredependent model for the acoustic properties of Weyburn
reservoir rock is developed from the integrated isotropic model. A single set of vertical fractures
is modeled with Hudson’s crack model and the elastic constants of the saturated, anisotropic rock
are calculated with Brown and Korringa’s equation and the fluid properties for Weyburn
reservoir. The seismic velocities are calculated for a range in pressure and pure phase saturation
conditions, to show the maximum sensitivity of the rock to changes in pressure and fluid
saturation. In general, the Marly zone exhibits greater changes in seismic properties with changes
in pressure and fluid composition than does the Vuggy zone.
As Tables 4.1, 4.3, 4.4, 4.5 and 5.2 show, the sensitivity to fluid and pressure changes in
the isotropic rock physics models from ultrasonic testing is similar to the sensitivity of the
integrated, anisotropic rock physics model. Why is this and are there any advantages of using the
more complicated anisotropic model? The values of porosity, dry density, and elastic moduli are
slightly different between the isotropic and anisotropic models. In addition, the anisotropic
model includes a set of fractures. These differences affect the absolute value of the calculated
seismic velocities, but are only second order effects in calculating the changes in seismic 89 velocities due to changes in fluid pressure and composition. This finding supports the use of a
reservoirzone rock model for timelapse studies at Weyburn Field. Although the porosity,
lithology, pressure, and fracture characteristics vary spatially in the reservoir, the first order
effects in the seismic difference data are caused by changes in pore fluid pressure and saturation.
The integrated, anisotropic rock physics model has several advantages over core or logbased models. By integrating geophysical logs with ultrasonic core data, the model is more
representative of the reservoir zones. This is a major improvement in calculating seismic
velocities, but less important for estimating changes in seismic velocities, as discussed in the
previous paragraph. The incorporation of fractures causes a static decrease in elastic moduli,
which increases the sensitivity of the rock to changes in fluid saturation. This effect is significant
for the Vuggy zone. Although it would improve the modeling of seismic attributes such as
∆VP/VS and shear wave splitting, the pressure dependence of crack density was not incorporated
into the model because the required data are lacking. The simplicity of the fracture model is
justified because predicted pore pressure changes under CO2 flood at Weyburn Field are small
(Section 8.2.2). The seismic interpretation in this research also focuses mainly on the Pwave
data (Chapter 9), which are less sensitive to vertical fractures than are the Swave data.
Another advantage of the integrated, anisotropic rock physics model is that it can be used
in further RCP research. For some reservoir processes, pressure changes may be significant and a
more detailed fracture model may be required. In this case, the framework of the anisotropic rock
physics model is in place and can easily be modified to include this. Timelapse studies should
be able to use changes in NMO ellipses, azimuthal AVO response etc. – signatures that cannot be
modeled isotropically. 90 CHAPTER 6
SYNTHETIC SEISMIC MODELING 6.1 Introduction
Synthetic seismic modeling is helpful in identifying the horizon of interest in seismic
data. In timelapse seismic reflection surveys, synthetics can be used to quantify the change in
the seismic data due to changes in reservoir properties. In its simplest form, modeling synthetics
requires a layered earth model of seismic impedance and a source wavelet. The modeling in this
chapter consists of two main parts: modeling of both P and Swave impedance under varying
pore fluid and pore pressure conditions, and modeling of Pwave reflection synthetics with fluid
substitution in the reservoir zones. The results indicate the vertical resolution and guide the
interpretation method of the timelapse seismic data at Weyburn Field (Chapter 9). 6.2 Seismic Impedance Modeling
Amplitudes of reflection events in seismic data are determined by the changes in seismic
impedance with depth in the earth (Eqs. 6.2 and 6.3). Seismic impedance, Z, is the product of
density, ρ, and seismic velocity, V: Z = ρV . (6.1) The anisotropic rock physics model discussed in Chapter 5 is used to model the variation of Pand Simpedance with varying pore fluid and pore pressure. Figures 6.1 and 6.2 present the Marly
zone P and S1impedance. Figures 6.3 and 6.4 present the Vuggy zone P and S1impedance.
The baseline (before CO2 flood) insitu impedance is between the values for pure brine and pure
oil saturation. The repeat (1 year after CO2 flood) insitu impedance varies spatially over the
reservoir and ranges from the values for pure CO2 to between the values for pure brine and pure
oil saturation. 91 10.0 6 2 PImpedance, 10 kg/m /s Marly Dolostone, Crack Density = 0.03, Confining Pressure = 22 MPa 9.5 9.0
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO 2
CO2 8.5
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 6.1 Variation of Marly Pwave impedance with fluid saturation and pressure at constant
crack density. Marly Dolostone, Crack Density = 0.03, Confining Pressure = 22 MPa 5.6 6 2 S1Impedance, 10 kg/m /s 5.8 5.4 Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO 2
CO2 5.2 5.0
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 6.2 Variation of Marly S1wave impedance with fluid saturation and pressure at constant
crack density. 92 Vuggy Limestone, Crack Density = 0.1, Confining Pressure = 22 MPa 12.8 6 2 PImpedance, 10 kg/m /s 13.0 12.6
12.4
12.2
12.0
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO 2
CO2 11.8
11.6
11.4
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 6.3 Variation of Vuggy Pwave impedance with fluid saturation and pressure at constant
crack density. Vuggy Limestone, Crack Density = 0.1, Confining Pressure = 22 MPa 6 2 S1Impedance, 10 kg/m /s 7.15
Brine: 85k ppm NaCl
Oil: API 29, GOR 30 L/L
Oil + 40% CO 2
CO2 7.10 7.05 7.00 6.95 6.90
8 10 12 14 16 18 20 Pore Pressure, MPa Figure 6.4 Variation of Vuggy S1wave impedance with fluid saturation and pressure at constant
crack density. 93 As shown in Figures 6.1 and 6.3, changes in fluid saturation have a greater effect on Pwave impedance of the reservoir zones, considering that the expected pore pressure change is less
than 5 MPa. The Marly zone S1impedance is more sensitive to pressure changes than to fluid
saturation. In this model, the Vuggy zone S1impedance is much less sensitive to fluid saturation
or pressure changes than is the Marly zone. This is because the elastic moduli of the Vuggy core
samples are higher than those of the Vuggy zone and exhibit less pressure dependence (Chapter
3).
To aid in simplifying impedance contrasts in geophysical logs, the logs can be blocked by
seismic units. Impedance values from the previous charts could be used for the effective
properties of the reservoir units in synthetic modeling. Examples of this are shown in Figures 6.5
and 6.6 for P and S1wave impedance, respectively. Fluid substitutions were performed in the
Marly and Vuggy zones only. As shown in Figures 6.5 and 6.6, the top of the Marly is an
impedance decrease, the MarlyVuggy interface is an impedance increase, and the bottom of the
Vuggy is an impedance decrease.
For thick beds, the seismic data are easily correlated with seismic impedance differences
at layer boundaries. For thin beds, tuning effects make this correlation more difficult, as
discussed in the Section 6.3.1. 94 1395
1400 Marly Dolostone,
φ = 0.24 Depth, m 1405
1410
1415 Vuggy Limestone,
φ = 0.10 Brine
Oil
Oil + 40% CO2
CO2 1420
1425 8 10 12
6 14 16 2 PWave Impedance, 10 kg/m /s Figure 6.5 Variation of Pwave impedance with depth, Weyburn reservoir model. Crack density
of Marly = 0.03, Vuggy = 0.1, PC = 22 MPa, PP = 15 MPa. 1395
1400 Marly Dolostone,
φ = 0.24 Depth, m 1405
1410
Vuggy Limestone,
φ = 0.10 1415
Brine
Oil
Oil + 40% CO2
CO2 1420
1425
4 5 6 7
6 8 9 2 S1Wave Impedance, 10 kg/m /s Figure 6.6 Variation of S1wave impedance with depth, Weyburn reservoir model. Crack density
of Marly = 0.03, Vuggy = 0.1, PC = 22 MPa, PP = 15 MPa. 95 6.3 Synthetic Seismograms
HampsonRussell AVO software is used to compute synthetic seismograms. HR AVO
uses a convolutional model implemented in the frequency domain to determine the offsetdependent seismic trace due to propagation of a source wavelet through a 1D layered earth
model (HampsonRussell Software Services, Ltd., 1994). Multiples and mode converted waves
are not included. The modeling options chosen in this research are simple: vertical incidence
(zerooffset), with no transmission losses or geometric spreading or added noise. This is so that
the synthetics can be best compared to the processed seismic data. The seismic data volume is
poststack, has been migrated to approximate zerooffset conditions, and has had amplitude
balancing to counteract transmission losses and geometric spreading.
The Zoeppritz equation for reflection coefficient, RC, at vertical incidence is: RC = Z 2 − Z1
,
Z 2 + Z1 (6.2) where Z is seismic impedance and the subscripts 1 and 2 represent the upper and lower layers,
respectively (Zoeppritz, 1919). The seismic trace, T(t) is computed by T (t ) = RC (t ) * W (t ) , (6.3) where W(t) is the source wavelet, RC(t) is the reflection coefficient series as a function of two
way travel time for seismic waves in the media, and * denotes convolution. 6.3.1 Correlation of Reflectors and Reflections
The correlation of reflectors (layers in the earth) and reflections (events in the seismic
data) is done in HR AVO through synthetic modeling. This aids in interpreting the reservoir
horizons in the seismic data. By changing the seismic properties of the reservoir and generating
synthetics, the changes in the seismic data can be estimated (Section 6.3.2).
Synthetic modeling requires geophysical logs and an estimate of the source wavelet. Ida
Herawati (personal communication, 2002) of the RCP provided the Pwave source wavelet. It
was estimated by leastsquares inversion of the seismic data in a 130 ms window around the
reservoir and three wells in the survey area. The wavelet is shown in Figure 6.7. The dominant
frequency based on the width of the main peak is approximately 33 Hz. 96 Figure 6.7 Pwave wavelet for RCP seismic survey (from Ida Herawati). The logs for generating the synthetic were made by combining geophysical logs from
wells 410416613 and 02042364. Well 410416613 is outside of the survey area (Figure 4.6) but
has continuous data from the near surface to more than 1000 m below the reservoir. Well
02042364 has the logs necessary for fluid substitution in the reservoir zone. The reservoir zone
data of well 02042364 were substituted for the reservoir zone data of well 410416613.
Figure 6.8 presents the geophysical logs and the Pwave synthetic seismic trace generated
from them. The largest event is the Mississippian unconformity, which is a large peak on the
seismic data. Because the Marly and Vuggy zones are thin beds, they are not represented by
separate events in the seismic data. The Marly zone is approximately 1/16 of the dominant
wavelength and the Vuggy zone is 1/8 of the dominant wavelength. The top of the Marly zone
corresponds to a trough. The middle of the Vuggy zone corresponds to the zero crossing of the
same wavelet. 97 Mississippian Unconformity
Ratcliffe Beds
Midale Evaporite
Marly
Vuggy
Frobisher Beds Figure 6.8 Pwave synthetic seismogram from wells 02042364 and 410416613, wavelet in
Figure 6.7. 98 6.3.2 Sensitivity of Synthetics to Fluid Saturation
Synthetic modeling is used to estimate the changes in the timelapse seismic data due to
fluid saturation changes. For the following Pwave synthetics, fluid substitutions were performed
on the reservoir zones and the seismic traces were generated using the composite log and wavelet
described in Section 6.3.1. The fluid substitutions are described in Section 4.3. Because the
emphasis is on differences in the synthetics, not the individual synthetics, the results from these
well logs are qualitatively representative of the RCP survey area. As discussed in Reasnor (2001)
spatial variation in layer thickness, fracture density, porosity, and lithology may have first order
effects on the seismic data. These variations are expected only to have second order effects on
the timelapse difference, however, and are therefore not modeled. The sensitivity calculations in
Chapters 4 and 5 for various isotropic and anisotropic models support this. Because the Pwave
data are interpreted in this thesis, only Pwave synthetics are generated.
Swave synthetics could be generated to show the sensitivity to fluid saturation or
pressure, using the relations presented in Figures 6.2 and 6.4 and blocked logs. The Swave
synthetics are expected to be similar to the Pwave synthetics for several reasons (Reasnor 2001).
The impedance changes generally correlate between the S and Pwave velocity profiles. The
reservoir zones are thin beds to shear waves. The Swave velocity is approximately 55% of the
Pwave velocity and the dominant frequency is about 40% of that of the Pwave, so the resolution
of the Swave data is less than that of the Pwave data.
The following figures show the individual synthetics after fluid substitution, with the
Marly and Vuggy zones delineated. The original synthetic for insitu fluid saturation is
subtracted from the synthetics for brine, oil, and CO2 saturations. The original insitu fluid
saturation varies with depth. It averages 47% brine and 53% oil in the Marly zone and 72% brine
and 28% oil in the Vuggy zone, based on analysis of geophysical logs for well 02042364. To aid
in identifying changes, the “difference synthetic” is amplified by six. Figure 6.9 presents the
synthetics for identical fluid substitutions in the Marly and Vuggy zones. Figure 6.10 shows the
results for fluid substitutions in the Marly zone only and Figure 6.11 shows the synthetics for
fluid substitutions in the Vuggy zone only. 99 Brine 1060 Oil CO2
Original Change from original synthetic, amplified by 6
Brine Oil CO2 Time, ms 1080 1100 Marly
Vuggy 1120 1140
Amplitude Figure 6.9 Change in Pwave synthetic trace due to fluid change in both Marly and Vuggy zones. Brine 1060 Oil CO2
Original Change from original synthetic, amplified by 6
Brine Oil Time, ms 1080 1100 Marly
Vuggy 1120 1140
Amplitude Figure 6.10 Change in Pwave synthetic trace due to fluid change in Marly zone. CO2 100 Brine 1060 Oil CO2
Original Change from original synthetic, amplified by 6
Brine Oil CO2 Time, ms 1080 1100 Marly
Vuggy 1120 1140
Amplitude Figure 6.11 Change in Pwave synthetic trace due to fluid change in Vuggy zone. As seen in Figures 6.9 to 6.11, the predicted changes in the seismic data are not unique
for changes in the fluid saturation of the individual reservoir zones. Although the maximum
difference changes slightly in time, the fluid changes in the Marly, Vuggy, or Marly and Vuggy
zones cause changes in the reflections in the seismic data corresponding to both the Marly and
Vuggy zones. Because of this, the seismic data cannot resolve changes in the Marly and Vuggy
zones separately. An increase in the bandwidth of the source wavelet is needed to improve the
resolution of the seismic data.
The maximum expected changes in amplitude due to fluid saturation changes are 10% to
20%, with fluid changes in the Vuggy causing roughly half the effect of fluid changes in the
Marly. Impedance decreases in either zone, as would be expected with CO2 injection, produce a
larger trough that correlates with the top of the Marly zone and a larger peak that is close to the
bottom of the Vuggy zone. To capture the whole change in the seismic data, a window should be
used in calculating RMS amplitudes of the reflectors corresponding to the reservoir zones. 101 6.4 Summary
P and Swave impedance of the reservoir zones are modeled as a function of fluid
saturation and pore pressure. The results can be used in synthetic modeling to estimate the
changes in the seismic data due to expected fluid and pressure changes. As an example, Pwave
synthetic modeling is performed to quantify changes in the timelapse seismic data due to fluid
saturation changes in the reservoir zones. The results show that the effects of fluid changes in the
Marly and Vuggy zones are not unique, but produce differences in the reflection amplitudes over
a greater time window than has been considered to correspond to that reflector. This is expected
because both the Marly and Vuggy zones are seismically thin beds. This suggests that to identify
timelapse difference anomalies in the seismic data, RMS amplitude be measured on a window
around the reflections corresponding the reservoir zones.
The results apply qualitatively to all Pimpedance changes in the reservoir zones. For Pwaves, a decrease in pore pressure would cause an increase in impedance and have an effect
similar to changing the fluid saturation to brine. An increase in pore pressure would have an
effect similar to changing the fluid saturation to oil or CO2 in the previous examples. Additional
synthetic modeling could be done to estimate the expected timelapse changes in the Swave data
or to investigate the sensitivity of the change in reflection amplitudes for reservoir zones of
varying thickness. 102 CHAPTER 7
ACOUSTIC PROPERTIES OF WEYBURN OIL – CO2 SYSTEM 7.1 Introduction
In previous modeling of fluid and rock properties (Chapters 2 to 6), FLAG 4.0 has been
used to model the acoustic properties of brine and hydrocarbons. This modeling has served to
estimate the sensitivity of the seismic properties of the reservoir rock to fluid and pressure
changes. End members were used in the fluid substitutions to estimate the maximum sensitivity
to fluids. While this is valuable for feasibility studies and qualitative interpretation of timelapse
seismic data, more accurate fluid modeling is needed for quantitative interpretation.
A different method of fluid modeling is needed for three reasons. First, there will be
multiple phase fluid saturation: brine, oil, and/or a CO2rich fluid phase. Fluid substitutions need
to be done for the estimates of fluid saturation and composition from reservoir simulation, which
are not in a format FLAG can use. Second, previous modeling was done for oil with up to 66%
mole fraction dissolved CO2. Below the minimum miscibility pressure of approximately 14 to 17
MPa (Section 2.2.4) or for higher concentrations of CO2, a CO2 phase enriched with light
hydrocarbons separates from the oil phase. The reservoir simulation (Chapter 8) shows that the
CO2 phase forms near the injectors. This CO2 phase cannot be modeled accurately with FLAG
because it behaves differently than the hydrocarbons in the database upon which the FLAG
relations are built. Third, FLAG cannot be used directly on the detailed compositional output
from the reservoir simulation.
FLAG 4.0 can be used to model the properties of brine from the reservoir simulation
output, but not the oil phase and CO2rich phase. FLAG requires API gravity, gas oil ratio, and
gas gravity to specify composition (Section 2.2.2), but the reservoir simulator provides mole
fraction of seven pseudocomponents (Table 7.1). This chapter discusses the compositional data
from the reservoir simulation and the options for modeling its acoustic properties. 103 7.2 Compositional Reservoir Simulation
Eclipse 300 has been used to model the EOR operations and match the production history
at Weyburn Field. The Eclipse simulator has been used in compositional mode with the PengRobinson equation of state (EOS) (Peng and Robinson, 1965) and three phase fluid treatment.
The PengRobinson cubic EOS is used for the “flash” calculations. Flash calculations
refer to the determination of the composition and relative quantities of the phases present at
specified total composition, temperature, and pressure. The parameters in an EOS model are
fully determined for given pressure and temperature conditions. From an EOS model, all the
thermodynamic properties of a mixture can be determined, including bulk modulus and density.
The oil and gas hydrocarbon phases are modeled with seven pseudocomponents, which
correspond to groups of compounds with similar properties. Their properties are fit to PVT data
of the total mixture, as described in PanCanadian (1994) and Schlumberger Geoquest (1996).
Molecular weight and critical properties are summarized in Table 7.1. Additional properties
include boiling temperature, critical Z factor, acentric factor, static volume shift, LBC viscosity,
parachor, reference density, and binary interaction coefficients. Table 7.1 Description of components in Weyburn reservoir simulation.
Component Molecular
Weight Critical Temp.,
°C Critical Pressure,
MPa Critical Volume,
L LN2C 18.57 100.3 7.14 0.0573 LCO2 44.01 11.8 5.98 0.109 LC2C 39.19 69.9 3.04 0.263 LC4C 79.32 196.2 4.85 0.227 MC8C 168.61 348.5 2.62 0.542 HC10 315.76 463.7 1.29 1.168 H3 647.84 502.9 1.05 1.655 An advantage to compositional simulation is that fluid phase changes can be modeled
accurately. Fluid phase changes may be the dominant cause of changes in timelapse data. 104 Development of a small amount of gas or compressible liquid phase in an otherwise oil or brinesaturated rock lowers the bulk modulus significantly. 7.3 Fluid Modeling Options
Fluid modeling options for oilCO2 systems range from simple to complex, with varying
accuracy. A simple approximation method is to average the end member properties of oil and
CO2 based on volume fraction (Duranti, 2001). This is artificial because the mixture behaves
differently than the separate components. Miscibility and phase behavior change with pressure.
On the other extreme, equation of state models may describe the phase behavior and acoustic
properties very accurately, but require detailed compositional information and are
computationally intensive. For timelapse seismic modeling, the fluid model must accurately
describe the significant changes in acoustic properties and still be practical to implement from
available data. 7.3.1 FLAG 4.0
FLAG 4.0 can be used to model the properties of CO2rich oil, after STRAPP is used to
determine the gas oil ratio (Section 2.2.4). The sevencomponent reservoir simulation data could
be converted empirically to the API gravity, gas oil ratio, and gas gravity for modeling with
FLAG. This process may produce a rough estimate for the oilrich phase but introduce
significant error for the CO2rich phase. Currently FLAG has separate models for oil and gas that
are difficult to apply in the region inbetween liquid oil and gas and the relations are not
optimized for CO2. 7.3.2 STRAPP
In order to model acoustic properties of fluids with an equation of state such as STRAPP,
the composition and properties of the components must be precisely known. Chemical analysis
of Weyburn crude oil samples is available for separator gas, separator liquid, recombined oil and
dead oil (Core Laboratories, 2000). The composition of live (recombined) oil samples is broken
down into nonhydrocarbon gases, alkanes up to C30+, napthenes, and aromatic hydrocarbons. 105 The average molecular weight of one separator liquid sample at 35 °C and 130 psi is estimated at
168. Average molecular weight of dead oil was obtained by freezing point depression, yielding
values of 219 and 364 for two samples.
In STRAPP, Weyburn oil composition was specified by 20 components up to C30H62
based on chemical analysis of recombined oil. The amount of gas components was adjusted until
the live oil average molecular weight matched that of the background oil in the reservoir
simulator, 157. This corresponds to a dead oil average molecular weight of approximately 210.
The composition is listed in Table 7.2. To the live oil composition, CO2 was added in 2.5%
molar increments, up to 100% CO2. Density and bulk modulus were calculated in STRAPP for a
constant temperature of 63 °C and pressures from 3 to 30 MPa in 0.5 MPa increments. Table 7.2 Composition of live Weyburn oil modeled in STRAPP.
Constituent Mole Percent Constituent Mole Percent N2
CO2
H2S
CH4
C2H6
C3H8
C4H10
C5H12
C6H14
C7H16 1.50%
0.91%
0.19%
5.94%
3.54%
7.19%
6.86%
6.28%
8.22%
7.29% C8H18
C9H20
C10H22
C12H26
C14H30
C16H34
C18H38
C20H42
C24H50
C30H62 6.94%
3.54%
3.60%
6.28%
5.61%
4.22%
3.23%
2.66%
4.98%
11.02% The following graphs present the effective properties of the oilCO2 system. In the two
phase region, bulk moduli were averaged with Wood’s (1955) equation (Eq. 2.6) and density
were averaged by volume (Eq. 2.7). Figures 7.1 and 7.2 show the variation of effective bulk
modulus and density, respectively, with pressure and mole fraction CO2. There are some
questionable points around 0.5 mole fraction CO2 and 10 to 15 MPa pressure due to numerical
instability in STRAPP. A few unreasonable density values in this region were removed. 106 Figure 7.1 Variation of effective bulk modulus with composition and pressure for Weyburn oil –
CO2 system. Temperature is 63 °C. Figure 7.2 Variation of average density with composition and pressure for Weyburn oil – CO2
system. Temperature is 63 °C. 107 The bubble point pressure line is visible in Figures 7.1 and 7.2. As expected, the bubble
point pressure increases with increasing CO2 content. The phase transition is steeper in bulk
modulus than density, similar to the FLAG models shown Figures 2.10 and 2.11. At high
pressure and high CO2 content, the phase transition from CO2rich oil to CO2rich oil and light
hydrocarbonenriched CO2 is not clearly visible in either the bulk modulus or density models. At
reservoir conditions, the density of CO2 is similar to that of Weyburn oil. Although the bulk
modulus of pure CO2 is much lower than that of Weyburn oil, the bulk modulus change from
CO2rich oil to lighthydrocarbonenriched CO2 is much smoother.
How would the acoustic properties of oilCO2 mixtures modeled in STRAPP compare to
a Reuss (isostress) or Voigt (isostrain) average of the end member properties of oil and CO2? A
direct comparison is difficult because volume fractions, not mole fractions, are required for Reuss
or Voigt averages. The Reuss average would produce a concave bulk modulus surface as a
function pressure and CO2 content, similar to a low pressure slice of the surface in Figure 7.1.
The Voigt average would produce a linear relationship between bulk modulus and CO2 content.
Neither average would reproduce the complex phase behavior that produces the alternating
concave and convex bulk modulus surface shown in Figure 7.1.
The density of oilCO2 systems could be more accurately described through endmember
averaging than bulk modulus can be. At Weyburn reservoir conditions of approximately 15 MPa
and 63 °C, the variation in hydrocarbon density is less than 20%, compared to two or three order
of magnitude change in bulk modulus from CO2 to oil. Still, endmember averaging would not
reproduce the density changes due to phase transitions. It would also not match the unexpected
result of density first increasing, then decreasing with increasing CO2 content at high pressures
(Figure 7.2).
While modeling a 20component system is useful for describing the behavior of oil and
CO2 mixtures, it is impractical for modeling the hydrocarbons in the reservoir simulation.
Theoretically, the seven pseudocomponents could be added to the STRAPP component library
and the previous calculations repeated for the sevencomponent system. This effort failed
because of inconsistencies in the component attributes required in the Eclipse and STRAPP
equations of state. STRAPP could not take advantage of the many of the component descriptors 108 from the Eclipse implementation of the PengRobinson EOS and required additional component
data that were not available from the simulator. 7.3.3 Reservoir Simulator
An ideal solution would be to use the equation of state in the reservoir simulator to
calculate the fluid density and bulk modulus directly. This is possible with the PengRobinson
equation of state. It would eliminate errors in translating the compositional data from reservoir
simulation into different formats for other empirical or EOS models. It would also ensure that the
phase treatment is consistent between the transport properties used for flow modeling and the
acoustic properties used in rock physics modeling. Eclipse can output phase density, and this is
used in modeling (Chapter 8). Unfortunately, calculating adiabatic bulk modulus would require
additional programming for the Eclipse software, and is outside the scope of this thesis research. 7.3.4 Cubic EOS built on STRAPP Database
The gas model for FLAG 4.0 is based on Van der Waals equation, one of the simplest
cubic equations of state. Gas composition is specified by gas gravity or average molecular
weight. Average molecular weight is an attractive descriptor of hydrocarbon composition
because it can be calculated from any component mixture for which either mass fractions or mole
fractions are known. The form of the cubic EOS has been used to extend the range of this model
to heavy gases, condensates, and light oils with average molecular weights up 140 (Brown, 2001).
Hydrocarbon mixtures, including oils, are specified by their average molecular weight. In this
cubic EOS, the Van der Waals ‘a’ and ‘b’ coefficients are empirical functions of molecular
weight and pressure, optimized for a hydrocarbon fluid database.
In anticipation of modeling the CO2oil system at Weyburn Field, the cubic EOS was
extended to use with medium oils with average molecular weights of up to 250. This work is
presented in the Appendix. Unfortunately, this cubic EOS produced unacceptable error in
acoustic properties of oil with more than 25% mole fraction CO2. This shows that dissolved CO2
behaves differently than the short chain alkanes (methane, ethane, butane…) commonly
associated with oil. This model would be applicable to Weyburn if the injected CO2 dispersed
rapidly and small concentrations of CO2 were present in the reservoir. The reservoir simulation 109 shows that at least in the Marly zone, the injected CO2 builds up to high molar concentration
(90%+). It is likely that any model calibrated to typical oil properties would not perform well for
this CO2 rich phase.
The cubic EOS was used to develop separate coefficient relations for a database
consisting only of the Weyburn oilCO2 system. The database was calculated using STRAPP and
is the same as that shown in Figures 7.1 and 7.2, with the exception that only data points from 5
MPa to 27 MPa were included and the properties of the fluids in the two phase region are left
separate. Figure 7.3 presents the variation of bulk modulus with pressure and average molecular
weight for each fluid phase. In general, lower average molecular weight corresponds to higher
CO2 content. Figure 7.3 Variation of singlephase bulk modulus with composition and pressure for Weyburn
oil – CO2 system. Temperature is 63° C. (Based on STRAPP modeling). Figures 7.4 and 7.5 present the variation of bulk modulus with average molecular weight
and pressure, respectively. 110 Weyburn Oil  CO2 System 1200 Bulk Modulus, MPa 1000
800
600
400
200
0
20 40 60 80 100 120 140 160 Average Molecular Weight, g/mole Figure 7.4 Variation of bulk modulus with average molecular weight for Weyburn oil – CO2
system. Temperature is 63 °C, pressure varies from 3 to 27 MPa. Weyburn Oil  CO2 System 1200 Bulk Modulus, MPa 1000
800
600
400
200
0
5 10 15 20 25 30 Pressure, MPa Figure 7.5 Variation of bulk modulus with pressure for Weyburn oil – CO2 system. Temperature
is 63 °C, average molecular weight varies from 44 to 205. 111 Bulk modulus does not vary continuously with average molecular weight and pressure.
The discontinuity is also shown clearly in the plot of the average molecular weight – pressure
range of the data (Figure 7.6). The gap is the separation between CO2rich oil and the light
hydrocarbonenriched CO2 phase. In other words, the fractionation of separate phases from the
singlephase oilCO2 mixture is not continuous in CO2 content. Weyburn Oil  CO2 System Pressure, MPa 25 20 15 10 5
20 40 60 80 100 120 140 160 Average Molecular Weight, g/mole Figure 7.6 Range of average molecular weights and pressures that can exist for single fluid
phases in the Weyburn oilCO2 system. Calculated with STRAPP. A detailed description of the optimization of the coefficients for the cubic EOS is
provided in the Appendix. The pertinent parts for obtaining fluid bulk modulus are summarized
here. The cubic EOS is based on Van der Waal’s equation of state for a nonideal gas P= RT
a
− 2,
V −b V (6.1) where P is pressure, R is the universal gas constant, T is absolute temperature, V is molar volume,
‘a’ is coefficient to account for intermolecular attraction, and ‘b’ is a coefficient to correct for
molecular volume. This equation is referred to as ‘cubic’ because rewriting the form to solve for 112 volume results in a cubic polynomial. Isothermal bulk modulus is defined as KT ≡ 1 ∂P
V ∂V T . (6.2) Rather than making the conversion between adiabatic (also called isentropic or dynamic)
and isothermal bulk modulus via the heat capacity ratios, an alternative approach is to assume the
same form and let the empirical coefficients make the conversion implicitly. This produces K= VRT
2a '
− 2,
2
V
(V − b' ) (6.3) where a’ and b’ are dynamic coefficients used for the dynamic bulk modulus in fluid acoustics.
The coefficients a’ and b’ are functions of average molecular weight, mw: a ' = a ' 0 + a '1 mw + a ' 2 mw 2
b' = b' 0 + b'1 mw + b' 2 mw 2 . (6.4) Average molecular weight is defined as
n mw = ∑ f i mwi , (6.5) i =1 where fi and mwi are the mole fraction and molecular weight of each component. For Eq. 6.3,
molar volume is obtained from the density, ρ, output from the reservoir simulation and the
relation: V= mw ρ . (6.6) The a’i, b’i coefficients are optimized through generalized nonlinear inversion on the oilCO2 molar volume and bulk modulus database computed in STRAPP. They are listed in Table
7.3. Figure 7.7 shows the bulk modulus predicted from the cubic EOS of state versus the bulk
modulus calculated with STRAPP. 113 Table 7.3 Optimum a’i, b’i coefficients for Eq. (6.5) to use in cubic EOS for Weyburn oilCO2
system. Coefficients correspond to units of L, MPa, K for molar volume, pressure and
temperature (R= 0.008314 L MPa/(Kelvin gmole)).
a' 0 a' 1 a' 2 b'0 b'1 b'2 4.15 x 101 7.87 x 102 1.34 x 103 1.61 x 102 1.69 x 103 1.08 x 105 Weyburn Oil  CO2 System Bulk Modulus Predicted from EOS, MPa 1200 1000 800 600 400 200 0
0 200 400 600 800 1000 1200 Bulk Modulus, MPa Figure 7.7 Comparison of predicted bulk modulus from cubic EOS with bulk modulus calculated
in STRAPP for Weyburn oil – CO2 system. The average error is in bulk modulus 9 MPa, although it is not evenly distributed. The
greatest errors are at the transition between CO2rich oil and light hydrocarbonrich CO2 phases.
For CO2rich oil, bulk modulus is overpredicted and for the light hydrocarbonrich CO2 phase it is 114 underpredicted. The maximum error in bulk modulus in this dataset is 152 MPa. Considering
that the effective fluid bulk modulus in any reservoir simulation cell is the average of the bulk
moduli of two or three fluid phases and that the effective elastic moduli of the Marly or Vuggy
zones in the reservoir simulation are averaged from six layers each, the error in fluid bulk
modulus from the cubic EOS is minor for purposes of this research.
The cubic EOS coefficients developed for the Weyburn oilCO2 system may not be
generally applicable to other oilCO2 systems. The form of Van der Waals equation is
theoretically based and highly nonlinear to be able to simulate fluid phase changes. For liquids,
the equation can become unstable if V and b are very close in value. On the compositional data
from the reservoir simulator, the cubic EOS predicts bulk modulus from –0.01 GPa to 1.2 GPa.
For an ideal gas, isothermal or static bulk modulus is equal to pressure; dynamic bulk modulus
for a real gas is slightly higher than static bulk modulus. The small negative values of bulk
modulus from the cubic EOS re replaced with small positive values equal to the pressure. The
range in bulk modulus estimated from the cubic EOS corresponds well with that predicted for the
oilCO2 models from FLAG (Figure 2.10). 7.4 Summary
In this chapter several modeling options for fluid bulk modulus of the Weyburn oilCO2
system are investigated. The acoustic properties of the oilCO2 system are complicated due to
compositional and phase changes that are not typical of live oil. The fluid modeling option that
works best is to use the density of the hydrocarbon phases as calculated by the EOS in the
reservoir simulator and to use an empirically optimized cubic EOS to predict bulk modulus. 115 CHAPTER 8
RESERVOIR SIMULATION AND SEISMIC ATTRIBUTE MODELING 8.1 Introduction
Reservoir simulations have been run by Sandy Graham of PanCanadian Petroleum and
Hiro Yamamoto of the RCP to predict and optimize the EOR operations. Integration of the
seismic and reservoir simulation data at Weyburn Field is a collaborative effort involving these
individuals and likely many others in future RCP Phase IX work. This thesis research analyzes
only the output of the reservoir simulation.
The uncertainty in the spatial variation of rock properties such as lithology, fracture
density, porosity and pore shape make it difficult to accurately estimate the seismic properties of
the reservoir over the RCP study area. Most of these uncertain rock properties remain constant in
the reservoir over time and therefore are largely cancelled out in the timelapse seismic data.
This makes it possible to identify the firstorder causes of differences in the timelapse seismic
data: changes in fluid pressure, saturation and composition.
Quantifying fluid changes or estimating reservoir properties from seismic data is an
inverse problem. In order to solve the inverse problem, the forward problem of estimating
seismic properties from reservoir properties must first be formulated. This can be done by
integrating reservoir simulation and rock physics modeling, as explained in this chapter.
The reservoir simulation results contain much of the information necessary to predict
changes in the seismic properties of Weyburn reservoir. Specifically, fluid pressure, saturation,
and composition data from simulation are combined with the rock physics models detailed in
Chapters 4 and 5 to estimate the seismic properties of the reservoir before and during CO2
flooding. Predicted timelapse changes in seismic attributes are compared with the seismic
difference anomaly in Chapter 9. In future work, this method for estimating timelapse seismic
differences from reservoir simulation can be used in refining the reservoir model and the
simulation process. 116 8.2 Reservoir Simulation
Reservoir simulation is used to model reservoir processes and predict future production at
Weyburn Field. Simulation has been used to evaluate the technical and economic feasibility of
enhanced oil recovery operations and in designing the CO2 flood process (Section 1.4).
The Weyburn Field CO2 flood area has been modeled, of which the RCP study area is a
subset. The simulator, Eclipse 300, has been used in compositional mode with the PengRobinson equation of state (EOS), three fluid phase, and single porosity options (Schlumberger
Geoquest, 1996). The RCP study area is covered by a square grid of cells (x = 1 to 42 , y = 21 to
60), each 60 m x 60 m. The survey area extends an additional 180 m SW and 60 SE. There are
15 layers in the reservoir model; 1 to 7 are for the Marly zone and 8 to 15 are for the Vuggy zone.
In the RCP study area, layers 1, 14, and 15 are not present. Cell thickness varies from layer to
layer and within layers.
The reservoir model was created from a detailed stratigraphic model built from
geophysical logs and core data. The reservoir model was then modified until the simulated
production closely matched the historical production. Reservoir pressure was not historymatched. For simulation of the CO2 flood, the reservoir simulation is run with the actual CO2
injection data. The model has not been adjusted to match the observed CO2 response in
production wells. 8.2.1 Summary of Reservoir Model
The 12layer reservoir model possesses much finer detail than can be resolved with
seismic data or that can be portrayed easily. For ease in presentation, the reservoir model and
simulation output have been averaged vertically into the Marly and Vuggy zones. Porosity, pore
pressure, fluid saturation, and fluid properties are averaged, weighted by pore volume. Key
properties in the reservoir model are summarized in Table 8.1. The spatial variation of these
properties is shown in Figures 8.1 through 8.8. 117 Table 8.1 Summary of reservoir properties in the reservoir simulation model.
Parameter Mean Minimum Maximum Marly thickness (m) 6.2 3.8 10.1 Vuggy thickness (m) 17.0 11.9 20.6 Marly porosity 0.220 0.097 0.363 Vuggy porosity 0.096 0.059 0.204 Marly xdir + ydir permeability, md 104 2.1 2454 Vuggy xdir + ydir permeability, md 126 15.2 1139 Figures 8.1 and 8.2 show that both the Marly and Vuggy thickness vary by a factor of
about two within the RCP study area. This may lead to strong tuning effects in the seismic data,
since both reservoir zones are seismically “thin” beds (Section 6.3.1)
As shown in Figures 8.3 through 8.8, the original geologic model used to build the
reservoir model was detailed and smoothly varying. The modifications to that model could be
described as “patches” with strongly contrasting properties. Although they may not fit into the
geologic framework at Weyburn Field, these modifications produce an adequate history match to
the production data. 118 m N 1 km Figure 8.1 Marly zone thickness in RCP study area, from reservoir simulation model. m N 1 km Figure 8.2 Vuggy zone thickness in RCP study area, from reservoir simulation model. 119 N 1 km Figure 8.3 Marly zone porosity, from reservoir simulation model. N 1 km Figure 8.4 Vuggy zone porosity, from reservoir simulation model. 120 md
1000 300 N 100 1 km
30 10 3 1 Figure 8.5 Marly zone xdirection permeability, from reservoir simulation model. md
1000 300 N 100 1 km
30 10 3 1 Figure 8.6 Vuggy zone xdirection permeability, from reservoir simulation model. 121 md
1000 300 N 100 1 km
30 10 3 1 Figure 8.7 Marly zone ydirection permeability, from reservoir simulation model. md
1000 300 N 100 1 km
30 10 3 1 Figure 8.8 Vuggy zone ydirection permeability, from reservoir simulation model. 122 8.2.2 Summary of Output Data
The reservoir simulation fluid saturation and pressure output from the time of the
baseline seismic survey (10/1/2000) and the repeat survey (10/2/2001) are used for modeling.
Table 8.2 summarizes the output at these times. To aid in interpreting the reservoir simulation
data, a map of the reservoir model with injection and production wells is shown in Figure 8.9. Table 8.2 Summary of output pressure and saturation data from the reservoir simulator. Baseline
refers to the time of the first seismic survey before CO2 injection and repeat means the time of the
monitor seismic survey, approximately one year after the commencement of CO2 injection. Gas
denotes the presence of a CO2based phase.
Parameter Mean Minimum Maximum Marly baseline pressure (MPa) 20.6 9.2 26.0 Vuggy baseline pressure (MPa) 21.6 15.8 26.1 Marly repeat pressure (MPa) 20.3 14.3 26.5 Vuggy repeat pressure (MPa) 20.9 15.7 26.7 Marly baseline oil saturation 0.49 0.39 0.64 Vuggy baseline oil saturation 0.38 0.35 0.50 Marly repeat oil saturation 0.50 0.39 0.64 Vuggy repeat oil saturation 0.41 0.30 0.59 Marly baseline water saturation 0.51 0.36 0.61 Vuggy baseline water saturation 0.62 0.50 0.65 Marly repeat water saturation 0.49 0.36 0.61 Vuggy repeat water saturation 0.59 0.38 0.67 Marly repeat gas saturation 0.0031 0 0.085 Vuggy repeat gas saturation 0.0024 0 0.081 123 N 1 km Vertical Well
Vertical H2O Injection Well
Horizontal Well
Horizontal CO2
Injection Well Figure 8.9 Reservoir simulation area with injection and production wells shown. Figures 8.10 through 8.17 present the key pressure and fluid saturation changes from
reservoir simulation. The pore pressure drops by approximately 2 MPa in the northwest half of
the survey area. The pore pressure also drops around the horizontal CO2 injection wells in the
southeastern side. In the Marly zone, pore pressure increases the most near the horizontal
production wells in the eastern quadrant. These were the lowest pressure areas at baseline and
were still low pressure areas at the time of the repeat survey (Figure 8.10).
Around the CO2 injectors, oil saturation increases because the CO2 dissolves in and
swells the oil. The regions with oil saturation increase generally correspond to the regions with
water saturation decrease. Although the water must be displaced, no distinct water front is visible
in the saturation change maps. 124 MPa N 1 km Figure 8.10 Marly zone pore pressure at time of repeat survey, from reservoir simulation. MPa N 1 km Figure 8.11 Vuggy zone pore pressure at time of repeat survey, from reservoir simulation. 125 ∆ MPa N 1 km Figure 8.12 Marly zone change in pore pressure, from reservoir simulation model. ∆ MPa N 1 km Figure 8.13 Vuggy zone change in pore pressure, from reservoir simulation model. 126 N 1 km Figure 8.14 Marly zone change in fractional oil saturation, from reservoir simulation. N 1 km Figure 8.15 Vuggy zone change in fractional oil saturation, from reservoir simulation. 127 N 1 km Figure 8.16 Marly zone change in fractional water saturation, from reservoir simulation. N 1 km Figure 8.17 Vuggy zone change in fractional water saturation, from reservoir simulation. 128 In the baseline simulation results there is no gas phase. Figures 8.18 and 8.19 show the
presence of a gas, or ‘vapor’ phase in the repeat results. The term ‘vapor’ is somewhat arbitrary
in the simulator as it refers to the less dense of the two phases from the EOS flash calculations.
Both phases may have densities and bulk moduli typical of liquid hydrocarbons and the ‘vapor’
phase may actually have a higher bulk modulus. The presence of the ‘gas’ phase in the
simulation is important because it shows that the CO2 content is high enough that a separate CO2based phase is present. The extent of the gas phase correlates well with areas of CO2 mole
fraction above 0.5, as shown in Figures 8.20 and 8.21.
Figures 8.22 and 8.23 show the change in fluid density from reservoir simulation. It is a
combination of both saturation and compositional changes. The fluid density changes in the
reservoir zones are small. They produce a maximum bulk density change of 1.5% in the Marly
zone and 0.6% in the Vuggy zone. ...
View
Full
Document
 Spring '10
 blllas

Click to edit the document details