Brown_2002 - INTEGRATION OF ROCK PHYSICS AND RESERVOIR...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: INTEGRATION OF ROCK PHYSICS AND RESERVOIR SIMULATION FOR THE INTERPRETATION OF TIME-LAPSE 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 time-lapse (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 time-lapse difference anomalies in the P-wave seismic data. The fluid physics models are based on existing empirical relations, laboratory measurements and equation of state modeling. The pressure- and porosity-dependent 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 P-wave reflection amplitude due to CO2 injection is 15% to 20%, and should be detected in the time-lapse seismic data. Through interpretation of P-wave seismic data volumes, areas effectively contacted by CO2 are identified. The observed time-lapse anomalies correlate strongly with the modeled CO2 movement and P-impedance 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 Pressure-Dependent, 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 P-Wave 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 4-D 9-C seismic survey area. Modified from PanCanadian Petroleum. .......... 7 Figure 1.5 EOR infrastructure at RCP 4-D 9-C 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 P-wave 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 P-wave 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 S-wave velocities with pressure for dry Marly sample PC 12-13 4626. Hysteresis can be seen in the variation between the increasing- versus decreasing-pressure measurements. ..................................................................................................................................... 28 Figure 3.5 Variation of ultrasonic P- and S-wave velocities with pressure for dry Marly sample PC 14-11 4659.5.................................................................................................................................................. 29 Figure 3.6 Variation of dry rock bulk modulus with pressure for Marly samples PC 12-13 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 12-13 4626..................................................................................................................................................... 30 Figure 3.8 Variation of Lame’s parameter and shear modulus with pressure for dry Marly sample PC 14-11 4659.5.................................................................................................................................................. 30 Figure 3.9 Variation of ultrasonic P- and S-wave velocities with pressure for dry Vuggy sample PC 14-11 4698..................................................................................................................................................... 32 Figure 3.10 Variation of ultrasonic P- and S-wave velocities with pressure for dry Vuggy sample PC 12-13 4643..................................................................................................................................................... 32 Figure 3.11 Variation of dry rock bulk modulus with pressure for Vuggy samples PC 12-13 4643 and PC 14-11 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 P-wave velocity with pore pressure at constant differential pressure for saturated Vuggy sample PC 14-11 4698. ............................................................................................ 36 Figure 3.15 Variation of ultrasonic S-wave velocity with pore pressure at constant differential pressure for saturated Vuggy sample PC 14-11 4698. ............................................................................................ 36 Figure 3.16 Variation of ultrasonic P-wave velocity with pore pressure at constant confining pressure for saturated Vuggy sample PC 14-11 4698. ............................................................................................ 38 Figure 3.17 Variation of ultrasonic S-wave velocity with pore pressure at constant confining pressure for saturated Vuggy sample PC 14-11 4698. ............................................................................................ 38 Figure 3.18 Variation of ultrasonic VP/VS with pore pressure at constant confining pressure for saturated Vuggy sample PC 14-11 4698............................................................................................................. 40 x Figure 3.19 Variation of shear splitting parameter with pore pressure at constant confining pressure for saturated Vuggy sample PC 14-11 4698. ............................................................................................ 41 Figure 3.20 Variation of shear splitting parameter with differential pressure for dry Marly samples PC 1213 4626 and 14-11 4659.5................................................................................................................... 41 Figure 4.1 Comparison of P-wave 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 P-wave 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 P-wave 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 P-wave impedance versus depth, as a function of fluid saturation, for well 02042364. ............................................................................................................................................................. 62 Figure 4.17 Blocked S-wave 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 (P-wave modulus) with crack density at constant differential pressure, calculated using Hudson’s crack model. ............................................................................................. 76 Figure 5.2 Variation of c44 (S1-wave modulus) and c66 (S2-wave modulus) with crack density at constant differential pressure, calculated using Hudson’s crack model. ........................................................... 76 Figure 5.3 Variation of c33 (P-wave modulus) with differential pressure at constant crack density, calculated using Hudson’s crack model. ............................................................................................. 77 Figure 5.4 Variation of c44 (S1-wave modulus) and c66 (S2-wave modulus) with differential pressure at constant crack density, calculated using Hudson’s crack model. ........................................................ 78 Figure 5.5 Modeled variation of c33 (P-wave 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 P-wave velocity as a function of fluid saturation and pore pressure at constant crack density and confining pressure. ........................................................................................................... 82 Figure 5.9 Marly S1-wave velocity as a function of fluid saturation and pore pressure at constant crack density and confining pressure. ........................................................................................................... 83 Figure 5.10 Marly S2-wave velocity as a function of fluid saturation and pore pressure at constant crack density and confining pressure. ........................................................................................................... 83 Figure 5.11 Vuggy P-wave velocity as a function of fluid saturation and pore pressure at constant crack density and confining pressure. ........................................................................................................... 84 xii Figure 5.12 Vuggy S1-wave velocity as a function of fluid saturation and pore pressure at constant crack density and confining pressure. ........................................................................................................... 84 Figure 5.13 Vuggy S2-wave 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 P-wave impedance with fluid saturation and pressure at constant crack density. ................................................................................................................................................ 91 Figure 6.2 Variation of Marly S1-wave impedance with fluid saturation and pressure at constant crack density. ................................................................................................................................................ 91 Figure 6.3 Variation of Vuggy P-wave impedance with fluid saturation and pressure at constant crack density. ................................................................................................................................................ 92 Figure 6.4 Variation of Vuggy S1-wave impedance with fluid saturation and pressure at constant crack density. ................................................................................................................................................ 92 Figure 6.5 Variation of P-wave 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 S1-wave 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 P-wave wavelet for RCP seismic survey (from Ida Herawati). .................................................. 96 Figure 6.8 P-wave synthetic seismogram from wells 02042364 and 410416613, wavelet in Figure 6.7. ... 97 Figure 6.9 Change in P-wave synthetic trace due to fluid change in both Marly and Vuggy zones............ 99 Figure 6.10 Change in P-wave synthetic trace due to fluid change in Marly zone. ..................................... 99 Figure 6.11 Change in P-wave 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 single-phase 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 oil-CO2 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 x-direction permeability, from reservoir simulation model................................... 120 Figure 8.6 Vuggy zone x-direction permeability, from reservoir simulation model.................................. 120 Figure 8.7 Marly zone y-direction permeability, from reservoir simulation model................................... 121 Figure 8.8 Vuggy zone y-direction 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 P-wave velocity, from reservoir simulation. ....................................... 136 xiv Figure 8.28 Vuggy zone change in P-wave velocity, from reservoir simulation. ...................................... 136 Figure 8.29 Marly zone change in P-wave impedance, from reservoir simulation.................................... 137 Figure 8.30 Vuggy zone change in P-wave impedance, from reservoir simulation................................... 137 Figure 8.31 Marly zone change in S1-wave velocity, from reservoir simulation....................................... 138 Figure 8.32 Vuggy zone change in S1-wave velocity, from reservoir simulation...................................... 138 Figure 8.33 Marly zone change in S1-wave impedance, from reservoir simulation. ................................. 139 Figure 8.34 Vuggy zone change in S1-wave 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 P-impedance from baseline to repeat survey, based on reservoir simulation and rock physics modeling. .............................................................................. 144 Figure 8.40 Weyburn reservoir percent change in S1-impedance 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 P-wave RMS amplitude for +/- 2 ms window around Marly trough reflection. ........ 154 Figure 9.2 Difference (Repeat – Baseline) in P-wave RMS amplitudes, window of +/-2 msec around Marly trough reflection. ............................................................................................................................... 155 Figure 9.3 Percent difference (Repeat – Baseline) in P-wave RMS amplitudes, window of +/-2 msec around Marly trough reflection. ........................................................................................................ 156 Figure 9.4 Difference (Repeat – Baseline) in P-wave RMS amplitudes, window of 3 msec above and 9 msec below Marly trough reflection.................................................................................................. 157 Figure 9.5 Percent difference (Repeat – Baseline) in P-wave RMS amplitudes, window of 3 msec above and 9 msec below Marly trough reflection........................................................................................ 158 Figure 9.6 Process for interpretation of time-lapse seismic data, using an integrated reservoir simulation and rock and fluid physics model...................................................................................................... 161 Figure 9.7 Weyburn reservoir percent change in P-impedance, 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 P-wave 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 time-lapse 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 pseudo-samples with the same average molecular weight. Data modeled in STRAPP. ................................................................................. 205 Figure A.34 Bulk modulus versus pressure relationships for two pseudo-samples 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 oil-CO2 system. Coefficients correspond to units of L, MPa, K for molar volume, pressure and temperature (R= 0.008314 L MPa/(Kelvin g-mole)).................................................................................................... 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 CO2-based 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 pseudo-samples 1 to 10 modeled in STRAPP. Values are mole fraction. .................................................................................................................................... 179 Table A.3 Composition of gas and condensate pseudo-samples 11 to 18 modeled in STRAPP. Values are mole fraction. .................................................................................................................................... 181 Table A.4 Composition of oil pseudo-samples 1 to 10 modeled in STRAPP. Values are mole fraction. 183 Table A.5 Composition of oil pseudo-samples 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, industry-sponsored research consortia at the Colorado School of Mines. The main goal of RCP is to integrate dynamic data from time lapse (4-D) multicomponent (9-C) 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 time-lapse 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 pressure-dependent, anisotropic rock physics model for the seismic response of a saturated, fractured, rock-fluid 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 4-D seismic data based on the reservoir simulation results. Fourth, the time-lapse seismic data are compared with the time-lapse 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 tidal-channel 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 oil-water 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 CO2-rich phase and maintains flow in the oil-rich 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), Marly-Vuggy 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 4-D 9-C seismic survey area. Modified from PanCanadian Petroleum. 1 km Figure 1.5 EOR infrastructure at RCP 4-D 9-C seismic survey area. 8 During the CO2 and water injection, the reservoir is monitored with several time-lapse multicomponent (4-D 9-C) seismic reflection and vertical seismic profile surveys. The baseline 9-C seismic reflection survey was shot immediately prior to CO2 injection in late September and early October 2000. The first 9-C 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 time-lapse 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 time-lapse 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 rock-fluid system, as in Gassmann’s (1951) or Brown and Korringa’s (1975) equations. Fluids affect the compressional-wave velocity, VP, of the rock-fluid 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 shear-wave 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 non-uniformity 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 two-phase 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 P-wave 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 P-Wave 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 P-wave 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 CO2-based 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 CO2-rich 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 time-lapse 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 large-scale 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 rock-fluid 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 rock-fluid system, tending to increase the P-wave velocity. Less compressible fluids are usually also denser. Increasing fluid density increases the density of the porous fluid-rock system and tends to lower P- and S-wave 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 S-wave 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 P-wave 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, P-wave 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 P-Wave 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 P-wave 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 pre-cleaned 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 water-saturated 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 P-wave (particle motion parallel to long axis of sample) and two orthogonal shear-wave 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 zero-time 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. P-wave and two S-wave measurements are made sequentially. The two shear-wave 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 P-wave velocity (VP) and 0.3 % for S-wave 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 calibration-corrected 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 fluid-rock 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 pressure-dependent 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 12-13 4626 and PC 14-11 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 12-13 4626 and PC 14-11 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 12-13 4626’ 0.29 21.7 1.93 2.88 PC 14-11 4659.5’ 0.333 34.8 1.815 2.72 Dry Marly Sample PC 12-13, 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 S-wave velocities with pressure for dry Marly sample PC 12-13 4626. Hysteresis can be seen in the variation between the increasing- versus decreasing-pressure measurements. S-Wave Velocity, m/s P-Wave Velocity, m/s 3400 29 Dry Marly Sample PC 14-11, 4659.5 ft 3400 P-Wave Velocity, m/s 2200 2100 3200 2000 3000 S-Wave 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 S-wave velocities with pressure for dry Marly sample PC 14-11 4659.5. Dry Marly Samples 17 16 Bulk Modulus, GPa 15 14 13 12 11 PC 12-13, 4626 ft, φ = 0.29 PC 14-11, 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 12-13 4626 and PC 14-11 4659.5. 30 Dry Marly Sample PC 12-13, 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 12-13 4626. Dry Marly Sample PC 14-11, 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 14-11 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 14-11 4698 and PC 12-13 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 14-11 4698 and PC 12-13 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 12-13 4643’ 0.115 1.3 2.305 2.60 PC 14-11 4698’ 0.133 103 2.34 2.70 Due to testing difficulties, only VS1 could be measured on sample PC 12-13 4643. Sample PC 14-11 4698 is more representative of the Vuggy reservoir units than is PC 12-13 4643. PC 12-13 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 14-11 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 14-11, 4698 ft 5100 2950 VP VS1 VS2 2900 4900 2850 4800 2800 4700 S-Wave Velocity, m/s P-Wave 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 S-wave velocities with pressure for dry Vuggy sample PC 14-11 4698. Dry Vuggy Sample PC 12-13, 4643 ft 5900 2800 VP VS1 2750 5700 2700 5600 2650 5500 S-Wave Velocity, m/s P-Wave 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 S-wave velocities with pressure for dry Vuggy sample PC 12-13 4643. 33 Dry Vuggy Samples 60 Bulk Modulus, GPa 55 50 45 PC 12-13, 4643 ft, φ = 0.11 PC 14-11, 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 12-13 4643 and PC 14-11 4698. Dry Vuggy Sample PC 14-11, 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 14-11 4698. Shear Modulus, GPa Lame's Parameter, GPa 24 34 Dry Vuggy Sample PC 12-13, 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 12-13 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 14-11 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 S-wave velocity, respectively, for saturated Vuggy sample PC 14-11 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 P-wave velocity for the different reservoir fluids is small, 1% to 2%. As expected the P-wave 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, P-wave 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 P-wave 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 14-11 4698, Differential Pressure = 6.9 MPa 5100 Brine: 50k ppm NaCl Oil + 39% CO2 CO2 P-Wave 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 P-wave velocity with pore pressure at constant differential pressure for saturated Vuggy sample PC 14-11 4698. Saturated Vuggy Sample PC 14-11 4698, Differential Pressure = 6.9 MPa Brine: 50k ppm NaCl Oil + 39% CO2 CO2 S-Wave Velocity, m/s 2750 2700 2650 0 5 10 15 20 25 Pore Pressure, MPa Figure 3.15 Variation of ultrasonic S-wave velocity with pore pressure at constant differential pressure for saturated Vuggy sample PC 14-11 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 S-wave 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 brine-saturated S-wave velocity, which is about 1% lower than the CO2 or oil + CO2 saturated S-wave velocity. Over the expected pore pressure range of 8 to 20 MPa, there is no variation in S-wave 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 S-wave velocity, respectively, for saturated Vuggy sample PC 14-11 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 14-11 4698, Confining Pressure = 34.5 MPa P-Wave 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 P-wave velocity with pore pressure at constant confining pressure for saturated Vuggy sample PC 14-11 4698. Saturated Vuggy Sample PC 14-11 4698, Confining Pressure = 34.5 MPa 2760 S-Wave 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 S-wave velocity with pore pressure at constant confining pressure for saturated Vuggy sample PC 14-11 4698. 30 39 The difference in P-wave 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 P-wave velocity for the sample saturated with a common fluid is approximately 1%. For the Vuggy sample tested, P-wave velocity does not uniquely distinguish fluid saturation from differential pressure effects. For example, the P-wave 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 S-wave velocity for different fluid saturations in the sample at the same pore and confining pressures (Fig. 3.17). The S-wave velocities for the brine saturated sample may be slightly lower because of the fluid density effect (Eq. 2.3). S-wave velocity varies by approximately 1% over a pore pressure range of 10 MPa for the Vuggy sample tested. For an isotropic rock, S-wave 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 time-lapse 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 S-wave 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 14-11 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 14-11 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 shear-waves, 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 shear-wave splitting parameter, and VS1 and VS2 are the fast and slow shear-wave 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 S-wave testing, so that the measured shear-wave splitting does not necessarily correspond to the true shear-wave 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 14-11 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 14-11 4698. Dry Marly Samples PC 12-13 4626 and PC 14-11 4659.5 Shear Splitting Parameter, γ, % 7 6 5 4 3 2 PC 12-13 4626 PC 14-11 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 12-13 4626 and 14-11 4659.5. 35 42 Figure 3.19 shows that shear-wave splitting increases with increasing pore pressure at a constant confining pressure. This behavior is interpreted as being caused by the opening up of preferentially-oriented 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 high-pressure injection well in the reservoir. Because the microcracks or reservoir fractures are preferentially oriented, shear-wave velocity varies with azimuth. As pore pressure increases, increased shear-wave splitting should be observed in the seismic data. In the Vuggy test sample, shear-wave 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 shear-wave splitting for the Marly core samples. Because shear-wave 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 shear-wave splitting from these lab data to changes observed in the seismic data. The lab data show that shear-wave 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 shear-wave 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 shear-wave 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 P-wave velocity from air- to brine-saturation and differential pressure from 10 to 25 MPa is 5% to 12%, consistent with sample PCP 14-11 4659.5 and PCP 12-13 4626. The variation of P-wave velocity of the ‘Vuggy Porous’ for the same conditions is 5% to 9%, higher than results from sample PCP 14-11 4698. The P-wave velocity of the ‘Vuggy Tight’ varies by less than 1%, consistent with results from sample PCP 12-13 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 S-wave 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 time-lapse 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 shear-wave splitting may indicate changes in fluid pressure. 44 Seismic attributes such as VP/VS and shear-wave 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. Shear-wave splitting is also sensitive to fracture density, direction, and aspect ratio. In time-lapse monitoring, using seismic difference volumes allows the use of changes in these attributes for interpretation. Whereas VP/VS or shear-wave 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 time-lapse 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 time-lapse 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 fluid-rock state) for wave propagation in saturated media; ultrasonic measurements may be in the high frequency (unrelaxed fluid-rock 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, P-wave velocity can be computed for direct comparison with ultrasonic Pwave measurements on core samples. This is shown for Vuggy sample 14-11 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 P-wave 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, P-wave 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 rock-fluid system being in the unrelaxed (high frequency) state in the ultrasonic measurements. More importantly for time-lapse seismic monitoring purposes, the magnitude of the changes in P-wave velocity with different fluid saturation and pressure are similar for the laboratory measurements as for Gassmann’s equation. 47 Saturated Vuggy Sample PC 14-11 4698, Confining Pressure = 34.5 MPa 5100 P-Wave 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 P-wave 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 P-wave 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 14-11 4698 and for the Marly zone, they are based on sample PC 12-13 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, brine-saturated rock has the highest P-wave velocity and VP/VS, followed by rock saturated by oil, and then CO2. The P-wave 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 P-wave 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 liquid-like 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 12-13 4626, φ = 0.29, PC = 22 MPa 3800 P-Wave 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 P-wave velocity with fluid composition and pressure, calculated with Gassmann’s equation. Marly Sample PCP 12-13 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 14-11 4698, φ = 0.11, PC = 22 MPa P-Wave 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 P-wave velocity with fluid composition and pressure, calculated with Gassmann’s equation. Vuggy Sample PCP 14-11 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 P-wave 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, P-wave velocity, and S-wave 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 12-13 Samples 14-11 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 P-Velocity, m/s S-Velocity, 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 S-wave 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 S-wave 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 P-wave 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 12-13 4626) φ = 0.29 (well 02042364) (PC 14-11 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 P-wave velocity for the reservoir zone with brine, oil, and CO2 saturation. Figure 4.11 presents the percent difference in P-wave 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 brine-saturated 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 brine-saturated 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 P-Wave 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 P-Wave 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 P-Wave Impedance, 10 kg/m /s Figure 4.16 Blocked P-wave 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 S-Wave Impedance, 10 kg/m /s Figure 4.17 Blocked S-wave 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 CO2-saturation. 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 pressure-dependent 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 12-13 4626 and for the Vuggy zone it is described by 14-11 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 porosity-dependent model for elastic moduli is therefore necessary to realistically perform fluid substitutions before upscaling the layers to the reservoir zones. This porosity-dependent 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 reservoir-wide 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 porosity-scaling 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 porosity-scaling 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 porosity-dependent 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 porosity-dependent 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 time-lapse 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 off-trend patterns. For this reason horizontal wells are oriented parallel to the apparent dominant fracture direction. Shear-wave 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 anhydrite-filled 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 drilling-induced 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. Shear-wave 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% shear-wave 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, fluid-filled, or filled with “weak” solids. Hudson’s model assumes that the cracks are isolated with respect to fluid flow, which corresponds to high-frequency fluid-rock 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 x1-axis (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 second-order corrections for cracks, respectively. For vertical wave propagation, the moduli of interest are c33 (P-wave), c44 (S1-wave), c66 (S2-wave). 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 x1-symmetry 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 (P-wave 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 (S1-wave modulus) and c66 (S2-wave 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 (P-wave 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 (S1-wave modulus) and c66 (S2-wave 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 shear-wave 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 (penny-shaped), 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 shear-wave 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 zero-offset P-wave volumes, the differences are minor. For the interpretation of P-wave 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 four-subscript notation is used (Aki and Richards, 1980) and summation is over α and β. A more intuitive formulation of Eq. 5.7, which uses the 2-subscript 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 P-wave velocity, increases with increasing fluid bulk modulus, as shown in Figure 5.5. Therefore fluid saturation should not change shear-wave 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 (P-wave modulus) with fluid bulk modulus at constant crack density and differential pressure. 5.4 Pressure-Dependent, 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 P-Wave 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 P-wave 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 S1-Wave 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 S1-wave 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 S2-Wave 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 S2-wave 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 P-Wave 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 P-wave 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 S1-Wave 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 S1-wave 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 S2-Wave 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 S2-wave 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 pressure-dependent 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 hydrocarbon-enriched CO2 phase are similar to those of pure CO2. The maximum expected difference in the reservoir is from 2-phase (brine-oil) saturation to 3-phase (brine-oil saturated with CO2-CO2 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 P-wave 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 time-lapse seismic reflection data and gauge their detectability. 5.5 Summary An anisotropic, pressure-dependent 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 reservoir-zone rock model for time-lapse 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 P-wave data (Chapter 9), which are less sensitive to vertical fractures than are the S-wave 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. Time-lapse 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 time-lapse 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 S-wave impedance under varying pore fluid and pore pressure conditions, and modeling of P-wave reflection synthetics with fluid substitution in the reservoir zones. The results indicate the vertical resolution and guide the interpretation method of the time-lapse 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 S-impedance with varying pore fluid and pore pressure. Figures 6.1 and 6.2 present the Marly zone P- and S1-impedance. Figures 6.3 and 6.4 present the Vuggy zone P- and S1-impedance. The baseline (before CO2 flood) in-situ impedance is between the values for pure brine and pure oil saturation. The repeat (1 year after CO2 flood) in-situ 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 P-Impedance, 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 P-wave impedance with fluid saturation and pressure at constant crack density. Marly Dolostone, Crack Density = 0.03, Confining Pressure = 22 MPa 5.6 6 2 S1-Impedance, 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 S1-wave 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 P-Impedance, 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 P-wave impedance with fluid saturation and pressure at constant crack density. Vuggy Limestone, Crack Density = 0.1, Confining Pressure = 22 MPa 6 2 S1-Impedance, 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 S1-wave 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 S1-impedance is more sensitive to pressure changes than to fluid saturation. In this model, the Vuggy zone S1-impedance 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 S1-wave 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 Marly-Vuggy 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 P-Wave Impedance, 10 kg/m /s Figure 6.5 Variation of P-wave 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 S1-Wave Impedance, 10 kg/m /s Figure 6.6 Variation of S1-wave 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 Hampson-Russell AVO software is used to compute synthetic seismograms. H-R 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 1-D layered earth model (Hampson-Russell Software Services, Ltd., 1994). Multiples and mode converted waves are not included. The modeling options chosen in this research are simple: vertical incidence (zero-offset), 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 post-stack, has been migrated to approximate zero-offset 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 H-R 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 P-wave source wavelet. It was estimated by least-squares 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 P-wave 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 P-wave 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 P-wave 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 time-lapse seismic data due to fluid saturation changes. For the following P-wave 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 time-lapse 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 P-wave data are interpreted in this thesis, only P-wave synthetics are generated. S-wave 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 S-wave synthetics are expected to be similar to the P-wave synthetics for several reasons (Reasnor 2001). The impedance changes generally correlate between the S- and P-wave velocity profiles. The reservoir zones are thin beds to shear waves. The S-wave velocity is approximately 55% of the P-wave velocity and the dominant frequency is about 40% of that of the P-wave, so the resolution of the S-wave data is less than that of the P-wave data. The following figures show the individual synthetics after fluid substitution, with the Marly and Vuggy zones delineated. The original synthetic for in-situ fluid saturation is subtracted from the synthetics for brine, oil, and CO2 saturations. The original in-situ 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 P-wave 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 P-wave 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 P-wave 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 S-wave 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, P-wave synthetic modeling is performed to quantify changes in the time-lapse 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 time-lapse 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 P-impedance 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 time-lapse changes in the S-wave 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 time-lapse 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 CO2-rich 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 CO2-rich 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 pseudo-components (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 Peng-Robinson 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 pseudo-components, 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 time-lapse 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 oil-CO2 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 time-lapse 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 CO2-rich oil, after STRAPP is used to determine the gas oil ratio (Section 2.2.4). The seven-component 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 oil-rich phase but introduce significant error for the CO2-rich phase. Currently FLAG has separate models for oil and gas that are difficult to apply in the region in-between 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 non-hydrocarbon 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 oil-CO2 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 CO2-rich oil to CO2-rich oil and light hydrocarbon-enriched 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 CO2-rich oil to light-hydrocarbon-enriched CO2 is much smoother. How would the acoustic properties of oil-CO2 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 oil-CO2 systems could be more accurately described through end-member 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, end-member 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 20-component 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 pseudo-components could be added to the STRAPP component library and the previous calculations repeated for the seven-component 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 Peng-Robinson 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 Peng-Robinson 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 CO2-oil 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 oil-CO2 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 single-phase 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 CO2-rich oil and the light hydrocarbon-enriched CO2 phase. In other words, the fractionation of separate phases from the single-phase oil-CO2 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 oil-CO2 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 non-ideal 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 non-linear 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 oil-CO2 system. Coefficients correspond to units of L, MPa, K for molar volume, pressure and temperature (R= 0.008314 L MPa/(Kelvin g-mole)). a' 0 a' 1 a' 2 b'0 b'1 b'2 -4.15 x 10-1 7.87 x 10-2 -1.34 x 10-3 -1.61 x 10-2 1.69 x 10-3 -1.08 x 10-5 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 CO2-rich oil and light hydrocarbon-rich CO2 phases. For CO2-rich oil, bulk modulus is overpredicted and for the light hydrocarbon-rich 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 oil-CO2 system may not be generally applicable to other oil-CO2 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 oil-CO2 models from FLAG (Figure 2.10). 7.4 Summary In this chapter several modeling options for fluid bulk modulus of the Weyburn oil-CO2 system are investigated. The acoustic properties of the oil-CO2 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 time-lapse seismic data. This makes it possible to identify the first-order causes of differences in the time-lapse 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 time-lapse changes in seismic attributes are compared with the seismic difference anomaly in Chapter 9. In future work, this method for estimating time-lapse 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 12-layer 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 x-dir + y-dir permeability, md 104 2.1 2454 Vuggy x-dir + y-dir 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 x-direction permeability, from reservoir simulation model. md 1000 300 N 100 1 km 30 10 3 1 Figure 8.6 Vuggy zone x-direction permeability, from reservoir simulation model. 121 md 1000 300 N 100 1 km 30 10 3 1 Figure 8.7 Marly zone y-direction permeability, from reservoir simulation model. md 1000 300 N 100 1 km 30 10 3 1 Figure 8.8 Vuggy zone y-direction 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 CO2-based 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

Ask a homework question - tutors are online