Multiphase, multiscale simulation of fractured reservoirs

Numerical simulation of the geometrically complex fractured reservoirs has been a major engineering challenge. The deficiencies of continuum models are often addressed using the discrete fracture network (DFN) models which represent the complex fracture geometry explicitly. The primary goal in this dissertation is to explore ways of applying the DFN methodology to solve a variety of multiphase problems in oil reservoir simulation. Three-dimensional, three-phase simulators using the control-volume finiteelement scheme were used. After completing validation and fracture-property sensitivity studies, the limitation of employing the often-used Oda homogenization method was shown followed by the development of a simpler geometric scheme. The important question of oil recovery from basement reservoirs (Type I) composed of fractures of various sizes was examined in detail. Oil recovery and breakthrough behavior of this system comprised of seismic and subseismic features were investigated for different oil distributions, permeability values, levels of heterogeneity and rate. In general having more oil distributed in smaller systems led to lower recovery and quicker breakthrough. Lower permeabilities in the subseismic features also led to lower recovery. The recovery at given pore volume of water injected was rate dependent in all of the scenarios explored, with the lower rate production leading to about 5% higher oil in place recovery. This phenomenon was consistent when viewed from the point of view of gravity number for each displacement. The mechanism of gravity-dominated oil recovery in two-phase applications was explored, and a "critical rate" concept for obtaining higher recoveries in gravity-dominated flow was developed A multiscale upscaling exercise was performed to match the oil recovery performance from a structured fault zone using a single feature with different sets of relative permeability curves. The effectiveness of using DFN simulations for reservoirs containing matrix and fractures (Type II) was shown using two different systems. It was shown that placing wells either in the fault zone or in the matrix can have significant impact on recovery and breakthrough behavior. It was also demonstrated that fracture networks bring apparent anisotropy, and water-flooding from one direction or the other may affect oil recovery. Fractured reservoir simulation is high-performance computing - data and file management, computation, visualization, etc. are integral components of this exercise. A workflow to facilitate creation of fracture networks, gridding and simulation, and visualization was developed. A fully integrated two-dimensional graphical user interface (java-based) was also built.

MULTIPHASE, MULTISCALE SIMULATION OF FRACTURED RESERVOIRS by Huabing Wang A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Chemical Engineering The University of Utah December 2008 Copyright © Huabing Wang 2008 All Rights Reserved THE U N I V E R S I T Y OF UTAH G R A D U A T E SCHOOL SUPERVISORY COMMITTEE APPROVAL of a dissertation submitted by Huabing Wang This dissertation has been read by each member of the following supervisory committee and by majority vote has been found to be satisfactory. 2- i J 2 <►* $ Chair: Mi^ac^C^Deo S' vx ig B. Forster Edward M. Trujillo (Vjulp/j. Magda ZCOQ> --I-t _______ PaufBorgmeier THE U N I V E R S I T Y OF UTAH G R A D U A T E SCHOOL FINAL READING APPROVAL To the Graduate Council of the University of Utah: I have read the dissertation o f ____________ Huabing Wang____________ jn its final form and have found that (1) its format, citations, and bibliographic style are consistent and acceptable; (2) its illustrative materials including figures, tables, and charts are in place; and (3) the final manuscript is satisfactory to the supervisory committee and is ready for submission to The Graduate School. O& / 6$jC% Date Milind D. Deo Chair: Supervisory Committee Approved for the Major Department X / ----- - ^ - w ~ [ / (^ y Jo Ann S. Lighty Q Chair/Dean ^ Approved for the Graduate Council David S. Chapman Dean of The Graduate School ABSTRACT Numerical simulation of the geometrically complex fractured reservoirs has been a major engineering challenge. The deficiencies of continuum models are often addressed using the discrete fracture network (DFN) models which represent the complex fracture geometry explicitly. The primary goal in this dissertation is to explore ways of applying the DFN methodology to solve a variety of multiphase problems in oil reservoir simulation. Three-dimensional, three-phase simulators using the control-volume finite-element scheme were used. After completing validation and fracture-property sensitivity studies, the limitation of employing the often-used Oda homogenization method was shown followed by the development of a simpler geometric scheme. The important question of oil recovery from basement reservoirs (Type I) composed of fractures of various sizes was examined in detail. Oil recovery and breakthrough behavior of this system comprised of seismic and subseismic features were investigated for different oil distributions, permeability values, levels of heterogeneity and rate. In general having more oil distributed in smaller systems led to lower recovery and quicker breakthrough. Lower permeabilities in the subseismic features also led to lower recovery. The recovery at given pore volume of water injected was rate dependent in all of the scenarios explored, with the lower rate production leading to about 5% higher oil in place recovery. This phenomenon was consistent when viewed from the point of view of gravity number for each displacement. The mechanism of gravity-dominated oil recovery in two-phase applications was explored, and a "critical rate" concept for obtaining higher recoveries in gravity-dominated flow was developed A multiscale upscaling exercise was performed to match the oil recovery performance from a structured fault zone using a single feature with different sets of relative permeability curves. The effectiveness of using DFN simulations for reservoirs containing matrix and fractures (Type II) was shown using two different systems. It was shown that placing wells either in the fault zone or in the matrix can have significant impact on recovery and breakthrough behavior. It was also demonstrated that fracture networks bring apparent anisotropy, and water-flooding from one direction or the other may affect oil recovery. Fractured reservoir simulation is high-performance computing - data and file management, computation, visualization, etc. are integral components of this exercise. A workflow to facilitate creation of fracture networks, gridding and simulation, and visualization was developed. A fully integrated two-dimensional graphical user interface (java-based) was also built. '■ v CONTENTS ABSTRACT............................................................................................................................ iv LIST OF TABLES................................................................................................................. ix LIST OF FIGURES.................................................................................................................xi ACKNOWLEDGMENTS.................................................................................................... xx 1. INTRODUCTION .............................................................................................................1 1.1 Discrete Fracture Network Modeling of Basement Reservoirs................................5 1.2 Modeling Basement Reservoirs on a Variety of Scales............................................ 7 1.3 Upscaling Relative Permeabilities in Basement Reservoirs.....................................7 1.4 Simulation of Reservoirs with Fractures and Matrix................................................ 8 1.5 Workflow...................................................................................................................... 9 2. FUNDAMENTALS OF FRACTURED RESERVOIR MODELING........................11 2.1 Fundamentals of Oil Reservoir Modeling................................................................ 11 2.2 Fractured Reservoir Classification............................................................................15 2.3 Modeling Fractures in Reservoir Simulation...........................................................16 2.4 Geological DFN Characterization.............................................................................20 2.5 Background of Control-Volume Finite-Element (CVFE) Simulator.....................22 2.6 Meshing....................................................................................................................... 39 2.7 Chapter Summary.......................................................................................................40 2.8 Nomenclature..............................................................................................................41 2.9 Bibliography...............................................................................................................42 3. SENSITIVITY STUDY OF OIL PRODUCTION FROM FRACTURED/FAULTED BASEMENT RESERVOIRS.......................................... 44 3.1 Abstract.......................................................................................................................45 3.2 Introduction.................................................................................................................46 3.3 Governing Equations and DFN Model.....................................................................47 3.4 Proposed Domain Characteristics.............................................................................49 3.5 CVFE Simulator Verification.................................................................................. 52 3.6 Applications................................................................................................................ 61 3.7 Chapter Summary......................................................................................................93 3.8 Nomenclature............................................................................................................. 95 3.9 Bibliography ............................................................................................................. 95 4. SIMULATING FRACTURED BASEMENT RESERVOIRS: COMPARING DISCRETE FRACTURE NETWORK MODELS TO THE UPSCALED EQUIVALENTS...................................................................... 97 4.1 Abstract....................................................................................................................... 98 4.2 Introduction.................................................................................................................99 4.3 Numerical Methods Overviews of Fracture Reservoir Simulation....................101 4.4 Proposed Domain Characteristics...........................................................................101 4.5 CVFE Simulator Verification.................................................................................103 4.6 Study of Hypothetical "Real" Basement Reservoir Model...................................103 4.7 Chapter Summary....................................................................................................126 4.8 Nomenclature.......................................................................................................... 127 4.9 Acknowledgments................................................................................................... 128 4.10 Bibliography........................................................................................................... 128 5. MULTIPHASE SIMULATIONS OF FRACTURED BASEMENT RESERVOIR COMPOSED WITH SEISMIC AND SUBSEISMIC SCALE FEATURES.....................................................................................................130 5.1 Abstract..................................................................................................................... 131 5.2 Introduction.............................................................................................................. 133 5.3 Governing Equations...............................................................................................135 5.4 Modeling the Fractured Basement Reservoir........................................................135 5.5 Case Studies..............................................................................................................139 5.6 Chapter Summary....................................................................................................164 5.7 Nomenclature............................................................................................................167 5.8 Bibliography ............................................................................................................168 6. MULTIPHASE SIMULATIONS AND UPSCALING OF EMBEDDED FRACTURED/FAULTED ZONE MODEL ON BASEMENT RESERVOIR.................................................................................. 170 6.1 Abstract.................................................................................................................... 171 6.2 Introduction...............................................................................................................172 6.3 Governing Equations, Discrete Fracture Network (DFN) Model and CVFE Simulator................................................................................... 175 6.4 Model Descriptions..................................................................................................175 6.5 Case Studies............................................................................................................. 180 6.6 Chapter Summary.................................................................................................... 206 6.7 Nomenclature............................................................................................................207 6.8 Bibliography ............................................................................................................208 7. MULTIPHASE FLOW STUDIES ON FRACTURED/FAULTED TYPE II RESERVOIRS WITH DISCRETE FRACTURE NETWORK APPROACH............................................................................................210 7.1 Abstract.....................................................................................................................211 7.2 Introduction.............................................................................................................. 212 7.3 Type II Fractured Reservoir Modeling.................................................................. 215 7.4 CVFE Simulator Verification with Pooladi-Davish and Firoozabadi's Fracture Water Level (FWL) Experiment.............................................................................216 7.5 Applications..............................................................................................................220 7.6 Chapter Summary.................................................................................................... 248 7.7 Nomenclature............................................................................................................249 7.8 Bibliography ............................................................................................................250 8. THE RESERVOIR SIMULATOR REMOTE COMPUTING INTERFACE DEVELOPMENT ............................................................................... 252 8.1 Overview...................................................................................................................253 8.2 Design........................................................................................................................253 8.3 Developments..........................................................................................................256 8.4 Definitions............................................................................................................. 261 8.5 Chapter Summary....................................................................................................262 9. DEVELOPING AND INTEGRATING FRACTURED RESERVOIR SIMULATION WORKFLOW........................................................... 263 9.1 Overview...................................................................................................................263 9.2 Fractured Network Characterization Packages......................................................264 9.3 Mesh Software Packages........................................................................................ 265 9.4 Assembling Simulation Input File..........................................................................266 9.5 Simulators.................................................................................................................266 9.6 Data Analysis........................................................................................................... 267 9.7 Visualization.............................................................................................................267 9.8 Chapter Summary....................................................................................................267 10. SUMMARY ,269 LIST OF TABLES Table Page 3.1 Fracture set statistics........................................................................................................ 54 5.1 Summary of seismic and subseismic features............................................................136 5.2 Summary of reservoir data for modeling seismic and subseismic featured basement reservoir.......................................................................................... 138 5.3 Summary of reservoir properties and operating conditions for case 1 set of simulations.....................................................................................................................140 5.4 Case 1 simulation summaries at the end of simulations (day 2,000)......................... 142 5.5 Summary of reservoir properties and operating conditions on case 2 studies.......... 144 5.6 Recovery comparisons from case 2 sets of simulations............................................. 144 5.7 Summary of reservoir properties and operating conditions on case 3 studies.......... 150 5.8 Summary of reservoir properties and operating conditions for case 4 sets of simulations.....................................................................................................................156 5.9 Case 4 simulation summary at the end of simulations (day 2,000)........................... 156 5.10 Summary of reservoir properties and operating conditions for case 5 sets of simulations....................................................................................................................159 5.11 Case 5 simulation summaries..................................................................................... 159 6.1 Summary of constructed embedded zone model.........................................................177 6.2 Summary of critical data for modeling embedded zone and upscaled single fracture basement model...............................................................................................181 7.1 Summary of critical data for modeling Pooladi-Davish and Firoozabadi's fracture water level experiment.................................................................................... 217 7.2 Summary of critical data for modeling type II reservoirs connected by normal fault....................................................................................................................224 7.3 Summary of critical data for modeling Teasdale Fault type II reservoir..................234 7.4 Summary of data for modeling type II reservoirs connected by normal fault.................................................................................................................... x LIST OF FIGURES Figure Page 2.1 A relative permeability curve for oil-water system.......................................................14 2.2 Schematic cross plot of percent reservoir porosity versus percent reservoir permeability (percent due to matrix versus percent due to fractures)......................... 16 2.3 Single porosity model in finite element mesh................................................................ 17 2.4 Dual porosity model illustrations (after Warren & Root, 1963)...................................18 2.5 DFN model illustrations...................................................................................................19 2.6 An example control volume mesh. The solid lines represent the triangulation T and the dashed lines represent the control volume (B.................................................24 2.7 A control volume with its boundaries across several triangular elements...................25 2.8 Decomposition of a control volume into several sub volumes.....................................25 2.9 Unit outward normals of subvolume bi:fn in triangle tm............................................... 26 2.10 Illustration of upstream nodes and flux directions.......................................................27 ( 2.11 Illustrations of modular reservoir simulator framework............................................. 31 2.12 Illustrations of reservoir simulation input and output information............................ 33 3.1 A control volume with its boundaries (dashed line) across several triangle elements and decomposition of a control volume into several subvolumes...............49 3.2 Regularized DFN basement model (left) and equivalent regularized ECLIPSE single porosity model.................................................................................... 50 3.3 ^regularized DFN basement model with the same model properties as regularized DFN basement model but different orientation and connectivity................................51 3.4 Fracture pole orientation (Wulff equal-angle projection, lower hemisphere) between hypothetical "real system" and "ideal system".............................................. 53 3.5 Comparisons of fracture size distribution (Power Law Plot) between the hypothetical "real system" and "ideal system".............................................................53 3.6 Three-dimensional view (left) and two-dimensional view (right) of discretized, regularized DFN modeled "ideal" basement reservoir domain for the CVFE simulator.............................................................................................................. 56 3.7 Three-dimensional view (left) and two-dimensional view (right) of discretized, regularized finite-difference modeled "ideal" basement reservoir domain for the ECLIPSE simulator......................................................................................................... 57 3.8 Comparisons of oil production rate between the DFN modeled CVFE simulator and the single porosity modeled ECLIPSE simulator for the regularized "ideal system"............................................................................................................................ 58 3.9 Comparisons of oil recovery between the DFN modeled CVFE simulator and the single porosity modeled ECLIPSE simulator for the regularized "ideal system"............................................................................................................................ 58 3.10 Snapshot of oil saturations for the regularized "ideal" basement domain from both CVFE simulator (left column) and ECLIPSE simulator (right column)..................59 3.11 Comparisons between feature/fault grids block thickness (left 20ft; right 0.5ft) by applying two constraints: OOIP & kh value..........................................................60 3.12 Irregularized DFN basement model with the same model properties as regularized DFN basement model with different orientation and connectivity............................ 61 3.13 At day 900, oil saturation distribution snapshots from different homogenous permeability models...................................................................................................... 62 3.14 Comparisons of water breakthrough and water cut among different homogeneous permeability models.............................................................................63 3.15 Illustrations of geometric mean based intrafault permeability distributions........... 64 3.16 At day 900, oil saturation distribution snapshots from different homogenous permeability models...................................................................................................... 65 3.17 Oil production rate comparisons among uniform 100 md and geometric mean 100 md cases........................................................................................................ 66 3.18 Cumulative oil production comparisons among uniform 100 md and geometric mean 100 md cases..................................................................................... 67 3.19 Water cut comparisons among uniform 100 md and geometric mean 100 md cases................................................................................................................. 67 3.20 Water cut comparisons for 33 33 33 geometric mean with random seeds.............68 3.21 Residual oil saturation distribution snapshot due to the fracture network connectivity (Case 1000 md uniform k with day 600).............................................. 70 3.22 Illustration of two separate fracture zones (separated by the acclivitous line in the middle of trace map) connected by four fractures (circles).................................71 3.23 Comparisons of oil saturation distributions between the original fracture network and disconnected model at day 900..............................................................72 3.24 Oil production rate comparisons between base case and reduced connections case..............................................................................................................73 3.25 Oil recovery factor comparisons between base case and reduced connections case............................................................................................................ 73 3.26 Illustration of single DFN modeled fracture orientation (pole and dip vector)............................................................................................................................75 3.27 Illustration of fracture strike orientations.................................................................... 75 3.28 Comparisons of fracture orientation distribution between base case and titled case by stereo-plots..............................................................................................77 3.29 Oil production rate comparisons between the base case and the dip angle case..................................................................................................................................78 3.30 Recovery comparisons between the base case and the dip angle case...................... 78 3.31 Oil concentration distribution comparisons between base case (top) and the tilted domain case (bottom).............................................................................79 3.32 Oil production rate comparisons among base case (full fracture height), half-top case and half-bottom case with same oil in place and well properties....... 81 3.33 Oil recovery comparisons among base case (full fracture height), half-top case and half-bottom case with same oil in place and well properties............................. 81 3.34 Oil saturation distribution comparisons at the end of the primary production (day 300).....................................................................................................83 3.35 Pressure distribution comparisons at the end of the primary production (day 300).........................................................................................................................83 3.36 Water cut comparisons among base case (full fracture height), half-top case and half-bottom case with same oil in place and well properties............................ 84 3.37 Oil saturation distributions at the end of simulation (day 900)...................... 85 3.38 Illustration of preferred pathway characterization (other than circled area) and preferred pathway with cut-off zone (circled)............................................................. 87 3.39 Oil saturation distribution comparisons among different fracture permeabilities (base case, preferred pathway, and preferred pathway with "cut-off zone")............ 89 3.40 Comparisons of oil recoveries for preferred pathway study...........................90 3.41 Dimensionless analysis o f oil recovery vs. injected pore volume............91 3.42 Dimensionless analysis of oil recovery vs. produced water ...........................92 4.1 Irregular "real" DFN basement model (bottom) and equivalent regularized "ideal" DFN basement model (top)........................................................................................102 4.2 Three-dimensional view (left) and two-dimensional view (right) of discretized, irregularized DFN modeled "real" basement reservoir domain for the CVFE simulator........................................................................................................................103 4.3 Simplification of "real" basement system (left) to "ideal" basement system (right).................................................................................................................104 4.4 Comparisons of oil production rate between the "ideal" and "real" system..............106 4.5 Dimensionless analysis of cumulative oil production rate between the "ideal" and "real" systems.................................................................................................................106 4.6 Snapshot of oil saturations for the "real" basement system and "ideal" basement system........................................................................................................... 107 4.7 At day 300 (end of set primary production), the gas saturation (left) and the oil pressure distributions (right)..........................................................................................108 4.8 At the end of simulation (day 900), the pressure distributions on "real" (left) basement model and the "ideal" (right) basement model..........................................109 4.9 "Homogenized" or "Equivalent" single porosity modeled "real" basement reservoirs........................................................................................................................I l l 4.10 Illustration of the block-K method for calculating the x-directional permeability tensor...................................................................................................... 114 4.11 Illustrations of homogenized DFN grid with different cell sizes.................117 4.12 Visualization of porosity on homogenized DFN grids................................117 4.13 Visualization of the permeability tensors on homogenized DFN grids (The grid block sizes chosen are the 40-ft cube, 25x25x5 equal blocks on a 1000ft x 1000ft x 200ft domain.).....................................................................118 4.14 Oil production rate comparisons among homogenized ECLIPSE simulations and CVFE simulation...................................................................................................119 4.15 Water cut comparisons among homogenized ECLIPSE simulations and CVFE simulation......................................................................................................... 120 4.16 Limitation of Oda's permeability tensor calculation with fine grid cells....121 4.17 Illustration of simplified permeability constrain calculations.......................122 4.18 Oil production rate comparisons among homogenized (simplified) ECLIPSE simulation and CVFE simulation................................................................................123 4.19 Water cut comparisons among homogenized (porosity only) ECLIPSE simulations and CVFE simulation..............................................................................123 4.20 Oil saturation distributions for simplified fracture homogenization simulation (at day 900) by ECLIPSE........................................................................................... 124 4.21 Fluid flow through fracture/fault represented by the equivalent grid blocks (DFN -* homogenized DFN -> homogenized grid blocks).....................................125 5.1 Illustration of the basement reservoir domain composed by seismic features and subseismic features............................................................................................... 136 5.2 Comparison of oil recovery for case 1 set of simulations........................................141 5.3 Comparison of water cut versus for case 1 set of simulations..................................141 5.4 Comparisons of oil recovery versus simulation time for case 2 studies...................145 5.5 Comparisons of water cut versus simulation time for case 2 studies with different time scales.....................................................................................................147 5.6 Depth dependent porosity and permeability.......................................................149 x v 5.7 Comparisons of oil recovery versus simulation time for case 3 studies...................151 5.8 Comparisons of water cut versus simulation time for case 3 studies.............151 5.9 Oil saturations at the end of simulation......................................................................152 5.10 Comparisons of oil recovery versus simulation time for case 4 sets of simulations......................................................... ..........................................................157 5.11 Comparisons of water cut versus simulation time for case 4 studies...................... 157 5.12 Comparisons of oil recoveries in case 5.................................................................... 160 5.13 Comparisons of water cut in case 5............................................................................160 5.14 Dimensionless analyses of water and oil productions in case 5............................. 162 5.15 Oil saturation comparisons at the day 2,000 (end of simulation)........................... 163 5.16 Oil saturation comparisons for OOIP characterization of 50/50 cases.................... 165 6.1 Illustration of seismic features and the fault zone model.......................................... 178 6.2 Can embedded fault zone fractured network be represented by single seismic feature in flow simulation? ............................................................................182 6.3 Comparisons of oil recovery between embedded zone model and single feature model with straight-line relative permeability.................................................183 6.4 Comparisons of water cut between embedded zone model and single feature model with straight-line relative permeability................................................183 6.5 Can the embedded fault zone model be represented by a single seismic feature model with different relative permeability?....................................................185 6.6 One example of Corey-exponent relative permeability curves...................................186 6.7 Comparisons of water cut between the embedded fault zone model and the Corey-exponent relative permeability represented single feature............................. 187 6.8 Comparisons of oil recovery between the embedded fault zone model and the Corey-exponent relative permeability represented single feature............................. 187 6.9 An embedded fault zone simulation shows geometric trapping of the oil.................188 XVI 6.10 Illustration of depth-dependent porosity and permeability models and their reference surface in the embedded fault zone model.......................................189 6.11 Comparisons of oil recovery among depth-dependent reservoir property models and nondepth-dependent reservoir property model in the embedded fault zone study...........................................................................................190 6.12 Comparisons of water cut among depth-dependent reservoir property models and nondepth-dependent reservoir property model in the embedded fault zone study.......................................................................... ......................................... 191 6.13 Residual oil saturations in the embedded fault zone depth-dependent studies....... 192 6.14 Comparisons of water cut between the embedded fault zone model and Corey-exponent relative permeability represented single feature with depth-dependent properties.........................................................................................193 6.15 Illustrations of ECLIPSE single feature model for rate sensitivity study................195 6.16 Comparisons of oil recovery with different production rates...................................196 6.17 Comparisons of water cut with different production rates........................................197 6.18 Comparisons of produced water correlated with oil recovery for various rates operations............................................................................................................198 6.19 Illustrations of bottom aquifer support and side aquifer support single feature models...............................................................................................................199 6.20 Comparisons of oil recovery between different boundaries with various rates......200 6.21 Comparisons of water cut between different boundaries with various rates........ ..201 6.22 Snapshots of oil saturation distributions from various production rates and boundary conditions..................................................................................................... 202 6.23 Water cut comparisons of rate sensitivity test with side water support...................204 6.24 Oil recovery comparisons of rate sensitivity test with side water support..............205 6.25 Comparison of oil recovery versus flood gravity number by operating rate variations (Rate from left to right are: 600, 400, 300, 200, and 100 STB/day)......205 7.1 A tetrahedral element with associated control volumes............................................. 215 7.2 Rock assemblies of the stacked blocks and the domain used in the modeling........216 7.3 Oil recovery comparisons between experiment results and simulation results on low water flooding rate (0.011 PV/hr) and high water flooding rate (0.160 PV/hr)................................... ....................................................................218 7.4 Oil recovery comparisons on different finite volume refinement with 0.160 PV/hr water injection rate..................................................................................220 7.5 Cross-sectional illustrations of normal and reverse dip-slip faults............................ 221 7.6 Illustration of three-dimensional heterogeneous type II fractured reservoir connected by normal fault........................................................................................... 222 7.7 Cross-section of type II reservoirs connected by normal dip-slip fault.....................222 7.8 Oil recovery comparisons between two case studies (10000 md and 100 md fault permeabilities)...................................................................................... 226 7.9 Oil production rate comparisons between two case studies (10000 md and 100 md fault permeabilities)...................................................................................... 227 7.10 Water cut comparisons between two case studies (10000 md and 100 md fault permeabilities) ....................................................................................................228 7.11 Snapshots of oil saturations in 10000 md fault permeability case with day 600 (upper) and 6000 (lower)..............................................................................230 7.12 Snapshots of oil saturations in 100 md fault permeability case with day 600 (upper) and 6000 (lower)..............................................................................230 7.13 Illustration of modeled area, outcrop of joint zone with high permeability features and heterogeneous type II flow model.........................................................232 7.14 Comparisons of oil recoveries between line drive 1 and line drive 2 ......... 235 7.15 Oil saturation comparisons between line drive 1 and line drive 2 at the end of primary production (day 600)..........................................................................236 7.16 Oil saturation comparisons between line drive 1 and line drive 2 at the end of simulation (day 6000)...................................................................................... 238 7.17 Oil production rate comparisons between two line drive case studies....................238 7.18 (Line drive 2 oil rate - Line drive 1 oil rate)/ Line drive 1 oil rate......................... 240 7.19 Illustration of preferred pathway characterization in type II fractured reservoir and preferred pathway with cut-off zone (circled)...................................243 7.20 Oil saturation distribution comparisons among different fault permeabilities on different time scales..................................................................... 244 7.21 Oil recovery fraction comparisons for preferred pathway type II fractured reservoir studies.......................................................................................... 245 7.22 Oil production rate comparisons for preferred pathway type II fractured reservoir studies................................ .........................................................................246 7.23 Dimensionless analysis of oil recovery vs. injection water......................................247 8.1 A snapshot of client-side remote fractured reservoir computing interface................258 8.2 Modular client-server interfaces for "real time on-line" reservoir simulations........ 259 9.1 An example of multiphase DFN simulation workflow............................................. 264 9.2 Illustration of workflow on fractured reservoir simulation.........................................268 i xix ACKNOWLEDGMENTS First and foremost, I would like to thank my advisor, Dr. Milind Deo, for his support, encouragement and guidance during this work. I want to thank him for being so supportive and for always giving insightful advice. This dissertation would not have been possible without his guidance. I would also like to thank Dr. Craig Forster and the rest of my dissertation committee members, Dr. Edward Trujillo, Dr. Jules Magda and Dr. Paul Borgmeier for their valuable comments and suggestions on my proposal and dissertation. I especially thank Craig for all his geological support during this work. I also wish to thank Dr. Tom Doe, and I appreciate his effort and time in characterizing DFN modeled fracture networks. I would also like to thank Golder Associates, Schlumberger, Altair, Sandia National Lab and Argon National Lab for their academic license of various software i packages being applied in this research. I have received a lot of help and encouragement directly or indirectly from all members of the Petroleum Research Center and Chemical Engineering Department at the University of Utah. I am truly grateful to everyone who helped me in the past five years. This research was generously supported by the United States Department of Energy through contract number DE-FC26-04NT15531 and CuuLong JOC. This support is greatly appreciated. . Finally, special thanks to my wife, Ning Zhuo, for her support and for giving me the joy of life, our son Jayden. I would also like to thank my parents-in-law for helping with the care of Jayden when he was bom. My parents and brother also deserve my gratitude for supporting and encouraging my education from the very beginning. xxi CHAPTER 1 INTRODUCTION Naturally fractured reservoirs occur worldwide and constitute an important reservoir type. About 60% of the world's remaining resources reside in such formations. The main feature that distinguishes naturally fractured reservoirs from conventional reservoirs is the presence of fractures. In naturally fractured reservoirs fluids exist in two interconnected systems: the rock matrix, which usually provides the bulk of the reservoir volume, and the highly permeable rock fractures. These types of reservoirs are well known as being both highly heterogeneous and complex and therefore may affect the oil and gas production in many ways. The fracture can enhance production if utilized properly, as in gas-induced gravity drainage of oil; or they can adversely affect oil production when the channeling paths of water or free gas break through early. A reduction in risk and an improvement in understanding of reservoir behavior will lead to enhanced profitability from under-exploited fractured fields. Improving the recovery from fractured reservoirs is an increasingly important task for many oil companies. The recovery from reservoirs where fractures dominate permeability is often a fraction of the resource recovered from conventional reservoirs in which matrix permeability dominates. The lower recovery and higher risks associated with recovery from a fractured reservoir and a result of the uncertainty in forecasting how various completion placements, water flood patterns, and tertiary recovery processes will actually perform in fractured reservoirs. Fractures do more than simply increase reservoir permeability; they fundamentally alter reservoir connectivity and heterogeneity. To improve oil and gas recovery in naturally fractured reservoirs, the reservoir model must be characterized correctly. Identifying, characterizing, and mapping the fracture network in terms of the fracture's properties will result in optimal reservoir management. The fractured reservoirs are usually classified based on what positive effects the fracture system provides to the reservoir's overall quality. Nelson's (2001) classification is well respected and useful to the exploring community because it presents a classification system which provides geologists/engineers with a fractured model standard. Below are the four types of Nelson's fractured reservoirs classification: Type I: Fractures provide the essential reservoir porosity and permeability. Type II: Fractures provide the essential reservoir permeability. Type III: Fractures enhance permeability in an already producible reservoir. Type IV: Fractures provide no additional porosity or permeability but create significant reservoir anisotropy (barriers). The first three types describe positive reservoir attributes of the fracture system. The fourth type describes those reservoirs in which fractures are important for both reservoir quality that they impart, and the inherent flow anisotropy and reservoir partitioning that they create. Three main methods have been applied to model the fractured reservoirs: single porosity model, dual porosity model and discrete fracture network model. Comparing with the other two models, the Discrete Fracture Network approach (DFN) models the geometry of the fracture network explicitly and provides a realistic way of modeling fractured reservoir performance (Dershowitz, et al., 1998). The approach consists of three general steps: 1 1) Fracture Data Analysis 2) Generation of Discrete Fracture Networks 3) Discrete Fracture Network Analysis Step one involves analyzing the information from a variety of data sources to derive the parameters needed for step two. These parameters include fracture locations (spatial model), size, shape, orientation, flow properties, and number of distinct fracture sets. Step two involves generating multiple discrete fracture networks based on the results of the data analysis. Some geo-statistical models will be applied during this process. In step three, these networks are analyzed to derive engineering information. These include simple geometric analysis such as the computation of fracture densities, and the computation of complex multiwell flow simulations. Typically, geologists use a forward modeling approach to iterate between steps two and three during the model calibration process. For example, the calculated fracture density or the calculated response from a transient well test is compared to field measurements. If there is not an acceptable match, the input parameters for the fracture network are changed and the network analysis task is repeated with the updated model. Based on the analysis of the fracture reservoir data, two types of probability distributions can be applied to generate the DFN model. These two types are scalar data probability distribution and directional data probability distributions. Examples of scalar distributions are uniform distribution, exponential distribution, normal distribution, lognormal distribution, normal of log distribution and power law distribution. There are four popular types of probability distribution directional vaiables: the univariate Fisher, the bivariate Fisher, the bivariate normal and the bivariate Bingham distribution. Questions often being asked by reservoir management regarding DFN modeled fractured reservoirs are: 1) Are DFN model appropriate representations of fractured reservoirs? How does oil recovery and other recovery parameters (water cut, etc.) compare to results from conventional approach? 2) Fractures appear at various scales. There are several methods of representations of these fractures' existence. What are fundamental displacements in these multi scale environments? 3) How does rock matrix interact with fracture networks so created? Some answers to the above questions could be discovered by the objectives of this research program as: 1) Benchmark the CVFE reservoir simulator with commercial finite difference reservoir simulator such as ECLIPSE (A Schlumberger product) by considering fracture networks of varying density; the idea is to highlight the physical mechanisms that can be represented more accurately using the DFN approach. Explore the similarities and differences in modeling fractured reservoirs using DFN, dual porosity and single porosity methods. 2) Study the impact of geometric aspects of fracture networks on the oil production from a basement reservoir. 3) Perform multiphase, multidimensional CVFE simulations on basement DFN models section to assess how best to represent fractures and fracture zones in basement reservoirs. 4) Study the Type II fractured reservoir by CVFE method; explore the impact of oil/gas production due to the presence of fracture networks. , 5) Create a workflow for CVFE reservoir simulation that includes preprocessing and postprocessing utilities. Create an operating system free interface to make the simulation run locally and/or remotely. Details to achieve the above goals are discussed below. 1.1 Discrete Fracture Network Modeling of Basement Reservoirs Different types of fractured reservoirs may need different simulation approaches. Simulation of Type I reservoirs in detail is considered in this project task. Geologically, Type I reservoirs are granitic/basalt formations. Practically all of the oil in these reservoirs resides in fractures or fault zones. Since the geometry of the fracture network is critical in understanding recovery from these reservoirs, these lend themselves well for treatment by using the DFN approach. The primary questions to be answered are: 1) Are DFN models accurate in modeling physics of displacement in these reservoirs? ' - 2) How do geometric aspects of fractured networks and other fracture characteristic affect displacements in these systems? To answer the first question posed, we constructed a regularized fracture network that is amenable to gridding with the ECLIPSE program. The ECLIPSE program is the conventional finite-difference based reservoir simulator that allows a comparison to be drawn with our CVFE simulation results thereby validating the CVFE simulator. In Chapter 3, this indexing verification will be discussed in greater detail. Once the accuracy of the DFN approach is ascertained, and its effectiveness compared with the conventional single-porosity approach, detailed sensitivity studies of the impact of fracture characteristics on production from basement reservoirs can be conducted. In Chapter 3 of this dissertation, this sensitivity is outlined in great detail. Permeability distributions in the fracture zone, dip orientation and aperture variations, etc. are included in this sensitivity study to help quantify fractured networks uncertainties. A typical fracture network is highly heterogeneous. In order to understand the impact of permeability heterogeneity between fractures, and within fractures, presence of preferred high permeability pathways between injectors and producers was examined in Chapter 3 of this dissertation. The extent to which the preferred pathway network leads to early breakthrough will be quantified to ascertain the risk of preferred flow pathway existence. The impact of shutting off the high permeability fractures using gel treatments on improving oil recovery is also studied to compare with the fracture networks without such treatments. For irregular networks, DFN approach will be used directly to obtain simulation results. It will also be possible to "homogenize" the irregular fracture domain and create "equivalent" Eclipse input files. This entails creating porosity and permeability fields, using one of the methods outlined in the review of upscaling methods in Chapter 4. This study helps to evaluate the effectiveness of the upscaling techniques, and also to compare the DFN approach to what can be considered "single-porosity" simulation method. This research program compared the performance of the DFN simulator to the dual porosity representation of the basement domain. These questions will be discussed and answered in Chapter 3 and 4. 1.2 Modeling Basement Reservoirs on a Variety of Scales Seismic measurements are used to map the locations of major features in basement reservoirs. These features can be brought into a DFN model and the multiphase flow process can be simulated. However, it is recognized that the seismic mapping technique is identify features larger than a certain size. The smaller subseismic features are not represented in a model based only on seismic measurements. Smaller subseismic features do exist in basement reservoirs. The important questions are: 1) What is the relative importance of the subseismic features in oil storage and displacement? 2) How does the production rate affect of the final recoveries? To answer the questions posed above, a conceptual multiscale basement model containing seismic and subseismic features is created and simulated. Greater details are discussed in Chapter 5. 1.3 Upscaling Relative Permeabilities in Basement Reservoirs In basement reservoirs, FMI logs and cores may provide additional details regarding individual fracture zones (seismic features). This information can be used in a fracture characterization environment (such as FRED) to create models of fractures within the fault zone. In this project, these types of models are titled "embedded zone" models. To improve oil recovery in fractured network it becomes necessary to ask how oil-water displacements behave in the embedded zone environment. Once that information has been ascertained, it becomes necessary to compare to its single feature representations adopted in full-field seismic and subseismic model. In Chapter 6, the answer to this question will be explored in greater detail by creating and simulating a number of embedded features. The upscaled representations are simulated using regular polynomial relative permeabilities, while the embedded zone itself is simulated using straight line relative permeabilities. The functionalities that give the "best-match" for oil recovery and/or water cut can be observed through this exercise. These relative permeability functionalities are considered "upscaled" relative permeabilities. 1.4 Simulation of Reservoirs with Fractures and Matrix In most practical fractured reservoir simulation studies, fractures act as conduits while the matrix contains most of the fluids. The geometric aspects of fracture networks in the context of these Type II reservoirs are studied in greater detail in Chapter 7. The detailed Type II reservoirs case studies in Chapter 7 include: 1) Understanding the importance of a "dip-slip fault" that separates and connects two lobes of a fractured reservoir. 2) Understanding the oil production impacts from "strike-slip fault zone" type II reservoirs. . 3) Analysis of preferred pathways in fracture networks on production performance. 9 1.5 Workflow Simulating fractured reservoirs and predicting oil recovery requires integration workflows. The integration workflow should include geologic knowledge, petrophysical properties representation, well modeling and well history input. Most times, required data are scattered in geologic models, production histories and petrophysical property generation programs. One of the main project tasks is to bring about the integration of the geologic model and reservoir simulation input file. "FracMan Reservoir Edition" (FRED) is a sophisticated general program for importing a number of geologic inputs into a geologic model. The program can create fracture sets with characteristics of all of the measured data. FRED is used in this project as the primary fracture generation tool. These fracture sets may be operational in a Type I reservoir environment or a Type II reservoir situation. In Type II fractured reservoirs, the fracture sets are embedded into a reservoir matrix. When creating reservoir simulation input files using the sophisticated fracture network, both these aspects need consideration. A workflow tool will be discussed in Chapter 9 with the following features: 1) Treatment of fracture only reservoir or reservoirs with fractures and matrix. 2) Meshing the domain created. 3) Assigning properties to the control volumes and/or elements created during the meshing exercise. 4) Creating a simulation input file. 5) Assigning wells and well operational parameters. 6) Describing well production histories. 7) Executing the simulator with the given input file. 8) Generating output data for production analysis and for visualization. 9) Visualizing the data. It may also be necessary to develop a geologist vision of what a fractured reservoir would look like. A modularized java tool which draws the fractured reservoir and most of its features is discussed in greater detail in Chapter 8. 10 CHAPTER 2 FUNDAMENTALS OF FRACTURED RESERVOIR MODELING In this chapter, some fundamentals of fractured reservoir modeling are introduced. At the beginning, some basic concepts of black-oil reservoir model are presented. Then other concepts of fractured reservoir modeling are described from both geological and fluid flow aspects. At the end of the chapter, some background of control-volume finite-element (CVFE) reservoir simulator are presented. This CVFE formulation was used to study DFN characterized fractured reservoir throughout this research program. 2.1 Fundamentals of Oil Reservoir Modeling 2.1.1 Black-Oil Model Crude oil may contain over thousand components which makes modeling all components neither possible nor meaningful. The black-oil model is the most widely adapted method in the reservoir simulation world. The assumptions of the black-oil model can be summarized as: • • At most three phases: oil (hydrocarbon liquid), water, hydrocarbon gas; • Oil phase consists of only oil and solution gas; • Gas phase consists of only free hydrocarbon gas; • Isothermal system; • Oil and gas phases can reach equilibrium instantaneous; The black-oil model was used throughout this research work 12 2.1.2 Porosity The porosity of a rock is a measure of the storage capacity (pore volume) that is capable of holding fluids. Quantitatively, the porosity is the ratio of the pore volume to the total volume (bulk volume). This important rock property is defined in Equation (2-1) as: Permeability is a property of the porous medium that measures the capacity and ability of the formation to transmit fluids. The rock permeability, k, is a very important rock property because it controls the directional movement and the flow rate of the reservoir fluids. It is used in Darcy's law to calculate fluid flux. For a single-phase flow, Darcy's law is written in Equation (2-2) as where u is the Darcy's velocity, <p is the fluid potential, and p is the fluid viscosity. The symmetric and positive-define permeability tensor k in three-dimensional space is given by Equation (2-3) as 0 = (pore volume) / (bulk volume) (2- 1) where 0 is the porosity. 2.1.3 Permeability k 11 V(p (2-2) (2-3) 13 2.1.4 Phase Saturation and Potential In multiphase flow problems, multiple types of fluids (at most oil, gas and water for black-oil model) can exist in any one pore. Phase saturation Si describes the fracture of pore volume occupied by phase 1 and is defined in Equation (2-4) as where Vi is the pore volume occupied by phase 1. To satisfy volume conservation, the summation of all phases' saturation should be equals to 1. As shown in Equation (2-5), phase potential is defined by where Pi is the phase pressure, p t is the phase density, g is the gravitational constant, gc is a universal gravitational conversion constant, and z is the elevation of the fluid phase in consideration. For multiphase flow problem, Darcy's law can be extended to Equation (2-6) as where the subscript I denotes an individual phases, k,i is the phase relative permeability. (2-5) 2.1.5 Relative Permeability (2-6) 14 As shown in Figure 2.1, the relative permeability is considered as a rock property and is a function of the phase saturation S/. For an oil-water system, the oil phase cannot move when the oil saturation is lower than the residual oil concentration (Sro). The water phase becomes mobile only when water saturation is higher than the connate water saturation (Siw). Capillary pressure Pc is defined as the pressure difference between the nonwetting phase and the wetting phase as a function of the (wetting phase) saturation. For oil/water systems in porous rock, oil is in general considered to be the least wetting phase. Therefore, the capillary pressure is defined in the Equation (2-7) as: 2.1.6 Capillary Pressure (2-7) 0 Water Saturation Figure 2.1 A relative permeability curve for oil-water system. where subscript o and w represent oil and water phases. In reservoir engineering, Pc is an important parameter for simulation studies (in particular for heterogeneous systems). 2.1.7 Formulation Volume Factor As shown in Equation (2-8), the reservoir fluid formation volume factor, 5/, is defined as the ratio of the volume of each reservoir fluid at the prevailing reservoir temperature and pressure to the volume of that fluid at the standard conditions. (2-8) v l £ T C where subscript I denotes an individual phase, RC and STC denote reservoir conditions and stock tank conditions. 2.2 Fractured Reservoir Classification Based on what positive effects the fracture system provides to overall reservoir quality, Nelson's (2001) classified fractured reservoirs as the following four types: Type I: Fractures provide the essential reservoir porosity and permeability. Type II: Fractures provide the essential reservoir permeability. Type III: Fractures enhance permeability in an already producible reservoir. Type IV: Fractures provide no additional porosity or permeability but create significant reservoir anisotropy (barriers). 15 16 Nelson's fractured reservoir classification is demonstrated in Figure 2.2. The first three types describe positive reservoir attributes of the fracture system. The fourth type describes those reservoirs in which fractures are important not only for reservoir quality, but also the inherent flow anisotropy and reservoir partitioning that they create. 2.3 Modeling Fractures in Reservoir Simulation Three types of main models formulations are popular for modeling and simulating flow through naturally fractured reservoirs: The single porosity model, dual porosity model, and DFN model. Fracture dominated 100% Fracture 100% % of Total porosity 100% Fracture ------------------ 100% Matrix Matrix dominated Figure 2.2 Schematic cross plot of percent reservoir porosity versus percent reservoir permeability (percent due to matrix versus percent due to fractures) 2.3.1 Single Porosity Model The single porosity model is a straightforward application, as shown in Figure 2.3, fractures are represented explicitly by a fine-scale mesh. The single porosity approach provides accuracy, but it is not practical due to a very large member of grids when considering even a few fractures. This is because of the big mesh size contrast ratios generated by the fractures and the matrix. This contrast ratio can reach over tens or hundreds of magnitudes which create numerical instabilities in multiphase flow simulation. - 2.3.2 Dual Porosity Model The dual porosity models are most widely used in large-scale fractured reservoir simulation. The dual-porosity system is represented by two different continua, one representing the porous matrix and the other representing the fractures. Fluid flow is primarily through the high permeability, low porosity fractures surrounding individual 17 Figure 2.3 Single porosity model in finite element mesh matrix elements. The matrix blocks contain a majority of the reservoir volume and act as sources or sinks to the fractures. A shape factor represent transfer function describes the mass transfer between the matrix and fractures. The shape factor a is defined in Equation (2-9) as <7= 4 ( i + 1 + ;|) <2-9) where Lx, Ly and Lz are described in the dual porosity model illustration of Figure 2.4. Despite the numerical efficiency, the dual porosity models have some drawbacks. The model is limited to sugar cube representation of fractured media. The shape factor is too simple to describe the fluid flow when gravity and viscous effects are involved. This approach also assumes that the medium to have dense closely connected fractured network and may not be very accurate when treating only a few fractures. 18 surround by fracturcs rcscrvoironc grid block Figure 2.4 Dual porosity model illustrations (after Warren & Root, 1963) 2.3.3 DFN Model Some deficiencies regarding conventional continuum model (dual-porosity) was discussed above. As a remedy, the DFN model provides relatively new mythologies for addressing some important needs. The discrete-fracture model is an alternative to the single porosity model. In the discrete fracture model, the dimensionality of fractures is reduced from n to (n-1). This reduction greatly decreases the computational time. Compared to the continuum models, there are many advantages of the DFN model: it accounts for the heterogeneity in fractures accurately; the performance of the method is not affected by very thin fractures; it can account explicitly for the effect of even a single fracture on fluid flow; there is no need to compute transfer functions. The DFN model is illustrated by Figure 2.5. Comparing with the single porosity model (Figure 2.3), the DFN model could be treated as a simplified single porosity model. 19 Figure 2.5 DFN model illustrations However, there are a number of unresolved challenges with the DFN model based flow simulations. The bottle neck of simulating DFN characterized fractured reservoir is the unstructured domain meshing. Low quality mesh will result in a significant mathematical problem during the flow simulation. 2.4 Geological DFN Characterization Throughout this research work, DFN characterized fracture networks were studied in various aspects. In this section, the methods of generating DFN models are discussed with extent details. 2.4.1 DFN Modeling Approach The DFN (Dershowitz et al., 1984) approach models the geometry of the fracture network explicitly and provides a realistic way of modeling fractured reservoir performance. The DFN approach consists of three general steps: • Fracture data analysis • Generation of discrete fracture networks • Discrete fracture network analysis Step one involves analyzing the information from a variety of data sources to derive the parameters needed for step two. These parameters are fracture locations (spatial model), size, shape, orientation, flow properties, and number of distinct fracture sets. Step two involves generating multiple discrete fracture networks from the results of the data analysis in step one. In step three, these networks are analyzed to derive engineering information. These include simple geometric analysis like the computation of fracture densities, as well as complex multiwell flow simulations. 20 Typically, a forward modeling approach is applied to iterate between steps two and three during the model calibration process. For example, the calculated fracture density or the calculated response from a transient well test is compared to field measurements. If there is not an acceptable match, the input parameters for the fracture network are changed and the network analysis task is repeated with the updated model. Most DFN models being studied in this work were characterized by FracMan Reservoir Edition (FRED) which is a product of Golder Associates. 2.4.2 Methods to Create DFN Model After the fracture reservoir data analysis, two types of probability distributions are usually applied to generate the DFN domain: scalar data probability distribution and directional data probability distributions. Examples of scalar distributions are uniform distribution, exponential distribution, normal distribution, lognormal distribution, normal of log distribution and power law distribution. For the directional data probability distributions, four types of major models are available in FRED: univariate Fisher distribution, bivariate Fisher distribution, bivariate normal distribution and bivariate Bingham distribution. These distributions are stated in terms of their probability density functions, for a variation about the mean direction. The actual direction chosen from a directional distribution is the composite of directional variation and mean directions. All scalar and directional data probability distribution equations discussed here refer to FRED Manual (Golder Associates, 2007). 2 1 2 2 2.5 Background of Control-Volume Finite-Element (CVFE) Simulator 2.5.1 Governing Equations As shown in Eq (2-10) to (2-12), in the three-phase black-oil model, the governing equations for each phase can be derived from the general continuity equations (Bird et al., 1960) by including porosity 0 in the cumulative term (Aziz and Serrai, 1979; Peaceman, 1977) where subscript o, w, and g represent oil, water and gas phases, q is the source term, Rs is the gas-oil ratio, S is the saturation, B is the formation volume factor. The LHS in above equations represent the flux term, and u denotes to Darcy phase velocity (by combining Equation (2-5) and (2-6)). In this three-phase black-oil model, capillary pressures coupling phase pressures are listed in Equation (2-13) and (2-14) as OIL: (2- 10) WATER: GAS: V ■ -1-U5 ) - + 0 ^ ) +fls<7o + <?/5 (2-12) PC( S J = PD - (2-13) PX Se) = ps - pc (2-14) and the volume conservation is defined in Equation in (2-15) as S0 + S w+ Sg = l (2-15) All equations shown above make the three-phase black-oil problem description complete. For nonfractured reservoirs or continuum modeled fractured reservoirs, this type of problem is usually solved by fmite-difference (FD) discretization of the partial differential equations. However, for DFN modeled complex fractured systems, finite-element discretization is required to solve the flow problem. A multiphase, multidimensional, upstream flux-weighted CVFE simulator is introduced in the rest of this chapter. 2.5.2 CVFE Formulations The CVFE formulation applied in this research is derived from a finite-element point of view with a focus on the explicit expression for local flux. This discretization method has the advantage of easy handling unstructured geometry and higher order of accuracy (Young, 1978; Fung et al., 1991; Sukirman and Lewis, 1994; Yang, 2003; Fu, 2005; Matthai et al., 2005). In three-dimensional spaces, the fractures are modeled as two-dimensional surfaces, and the matrix is modeled as three-dimensional solid tetrahedrons in three-dimensional space. The basic concept of upstream flux-weighted CVFE is to use the fluid potential values for flux calculation; the flux so obtained is then used for mass balance calculations. To illustrate the concept of CVFE formulation, taking the example of the two-phase (water, oil), two-dimensional flow problem, governing equations derived above could be rewrote in Equation (2-16) as 0 = -V-p,U,+ - (2-16) where / represents phases. Figure 2.6 shows that the residual function for a control volume bi E fB with boundary T could be wrote by integrating Equation (2-16) to get Equation (2-17): 24 term in Equation (2-16) are omitted. As shown in Figure 2.7, a control volume is usually distributed across several triangular elements. During the computation, this distribution makes it difficult to evaluate the residual equation (Equation (2-17)). To solve this problem, an element-by-element add up method is applied to calculate the contributions from subvolumes (as Figure 2.8 shows). where n represent the unit outward normal on T'. For clarity, the phase term and source Figure 2.6 An example control volume mesh. The solid lines represent the triangulation T and the dashed lines represent the control volume (B Mathematically, the residual equation for the control volume bj in Figure 2. could be calculated by p i _ y 6 p i 1 £-im= 1 1 rn where subscript m denotes the sub volume of control volume 6, in Figure 2.8, represents the subvolume residual function (shown in Equation (2-19)) which is contributed by bi m and could be derived from Equation (2-17) as = 4 „ p u ' f i d s + Figure 2.7 A control volume with its boundaries across several triangular elements (2-18) (2-19) Figure 2.8 Decomposition of a control volume into several subvolumes As shown in Figure 2.9, for a subvolume bt in a triangle tm, the subvolume residual function (Equation (2-19)) is rewritten for bi>m to Equation (2-20) = f _ p u ■ n ds -f L - p u - n d s + f °®pS} dx (2-20) ‘ vi j *• Ojyjjj or where n is the unit outward normal of the corresponding boundary as shown in Figure 2.9. The first two terms in RHS of Equation (2-20) represent the flux flowing out of bim through the boundary cl] and ckl. These two flux terms can also be defined in Equation (2-21) as fi,cTj = 4 j P « '» ds and f iM = J _ ^ u • n ds (2-21) The flux-based upstream-weighting scheme (Yang, 2003) was applied in the CVFE formulation in this research program. The concept of this flux-based upstream weighting in CVFE is explained in the in the following. k Figure 2.9 Unit outward normals of subvolume bi m in triangle tm 27 Figure 2.10 demonstrates a triangle with constant flux across the three control volume boundaries. The upstream direction is determined by the sign of the each flux. From the flux definition in Equation (2-21), while flux /? > 0, the flux is flowing out of the control volume i. Therefore, for a control volume, the flux-based upstream operator up(ij) is defined in Equation (2-22) as The above flux-based upstream-weighting scheme in CVFE formulation was validated and verified for mass balances both global and locally by Yang (2003). Discretization of PDE in the last section results in lots of nonlinear equations. To solve the flow problem, global residual functions need to be linearized for the system. A nonlinear solver is required to generate the Jacobian matrix and assembly the multiphase residual vector. The simulator needs a linear solver to solve the optimized sparse system. Some concepts regarding CVFE reservoir simulator implementation are discussed below. (2-22) 2.5.3 CVFE Simulator k J Figure 2.10 Illustration of upstream nodes and flux directions 2.5.3.1 Assembly of Global Residual Functions As discussed in Section 2.5.2, partial residual functions need to be added up by the method of element-by-element. This element-by-element process utilizes local numbers of element for the computation convenience. However, for a reservoir simulation study, every node in the domain is assigned a unique ID. Compared with the local element-by-element number, this unique ID is called a global number. Therefore, a mapping between the local and global numbering system becomes necessary during the simulation study. This mapping system is normally accomplished by the local-to-global (LTG) array (Yang, 2003). By this LTG array, the global residual vector and Jacobian matrix can be assembled practically. 2.5.3.2 Nolinear Solver To solve the set of highly nonlinear global residual functions (Equation (2-18)) by CVFE discretization, the quasi-Newton framework is adopted in the reservoir simulator which was applied in this research program. In numerical analysis, Newton's method is a well-known iterative type of method for solving nonlinear equations. If the initial guess value is close to the root, Newton's method can often converge remarkably quickly. However, far from the desired root, Newton's method can easily lead an unwary user astray with little warning. Thus, good implementations of the method embed it in a routine that also detects and perhaps overcomes possible convergence failures. For the quasi-Newton framework in the CVFE simulator, Newton's method is always applied first. If the full Newton step improves the approximation, it is accepted. If not, a global line search is applied. Detailed information regarding quasi-Newton framework refers to Yang's work (2003). After this nonlinear process, the derivative of residual functions will 28 be assembled in a Jacobian matrix. The analytical form of the derivative should be used whenever possible for the best performance. In the reservoir simulation, most of the time the functions of relative permeability kr, formation factor B and viscosity are available, and their analytical form can be easily obtained. Compared to numerical derivatives, analytical derivatives are remarkable quickly. However, the problem associated with the analytical derivatives is that its derivation process is prone to human error. Simulators used in this research program have the capability of calculating derivatives both numerically and analytically. 2.5.3.3 Linear Solver After the process of linearization, obtained set of linear equations need to be solved using an appropriate conjugate-gradient type of method. For the CVFE simulator used in this research program, the linear solver from Portable, Extensible Toolkit for Scientific Computation (PETSc) is adopted to solve the set of linearized equations obtained from the last section. PETSc was developed at Argonne National Laboratory in the Mathematics and Computer Science Division as a general purpose suite of tools for the scalable solution of PDEs and related problems. PETSc provides a rich set of Krylov subspace methods as the linear solvers, for example, generalized minimal residual (GMRES), conjugate gradient (CG), conjugate gradient squared (CGS), and so on. The GMRES method is used in this research. Furthermore, PETSc offers several preconditioners to help optimize the matrix. Those preconditioners include Additive Schwartz, Block Jacobi, Jacobi, ILU, ICC, and so on. The ILU preconditioner is adopted by the CVFE simulator in this research work. The robustness and efficiency of PETSc have been widely tested and recognized. There are several features from PETSc that make it very convenient for the application programmer. For example, users can create complete application programs for the parallel solution of nonlinear PDEs without writing much explicit message-passing code themselves. Furthermore, PETSc enables a great deal of runtime control for the user without any additional coding cost. The runtime options include control over the choice of solvers, preconditioners and problem parameters as well as the generation of performance logs. Some basic PETSc functions used in the CVFE simulator are briefly introduced below. Detailed information refers to PETSc Users Manual (2007). VecSetValues(...): insert residual function values into a PETSc vector; VecNorm(...): calculate the norm of a vector; I MatSetVaIues(...): insert gradients into the Jacobian matrix; SLESCreate(...): create a PETSc linear solver; KSPSetType(...) : set method to solve the linear system; " SLESSolve(...): solve the set of linear equations 2.5.3.4 Modular Reservoir Simulator Due to its complexity, modem reservoir simulator developing work usually requires team effort which makes the object orientated concepts widely adopted. Supported by the U.S. Department of Energy (DOE), the Utah Finite Element Simulator (UFES) is a framework that has been developed during the past decade. The CVFE simulators applied in this research were developed from this UFES modular framework. 30 Figure 2.11 shows that the main structure of the framework consists of three major modules: discretization method (DM), physical model (PM), and Utility. The selection of discrete scheme (for example, CVFE in this research) and physical model (black-oil model for this research) determines how DM and PM are implemented. The Utility module provides "in-house" libraries and some special routes for interfacing the external libraries and the framework. The examples of the implementation layer in DM, the Interface policy provides driving forces between control volumes; the LineSouce policy provides the drive forces between control-volume and wells; the ControlVolume policy provides the size of the control-volume; the DiscretizationMethod policy provides the list of control-volume, line Figure 2.11 Illustrations of modular reservoir simulator framework source, and the connectivity map; the Object policy provides the control-volume and line source lists. On the implementation layer of PM, the PhysicalModel policy requires an implementation to decide how to get the solution of each time step and time step control; the ControlVolume policy requires the implementations to compute residual functions of each conservation equation, constraint equation and their derivatives; the LineSource policy requires the implementations of how the flows and their derivatives being computed between control-volumes and wells. Numerous submodules are included in the PM such as rock, fluid, rock-fluid modules, etc. The rock module provides information such as density, porosity, and permeability. In summary, during the simulation, all physical related behaviors will be provided by PM, and all numerical issues are governed by DM. ■ Based on this modular framework, several reservoir simulators were developed at the Petroleum Research Center (PERC) at the University of Utah. From the DM aspect, these simulators include CVFE, mixed finite-element (MFE), and finite-difference (FD). From the PM aspect, simulators include the black-oil model, thermal model, and compositional models. As mentioned above, the simulators applied in this research work are black-oil modeled simulators implementing CVFE discretization method. 2.5.4 Simulator Input and Output Reservoir simulation can be described as using computer system to study fluid flow in a reservoir. For a simulation job, various disciplines contribute to the preparation of the input data set. These contributions need to be integrated during the reservoir modeling process. The contribution from different disciplines make to reservoir modeling 32 is illustrated by Figure 2.12. The simulator is the contact point among disciplines. The integrated data set (reservoir domain information, properties of rock and fluid, well trajectories and operating conditions, production data and observation, reservoir initial conditions and boundary conditions) is sent into the simulator as the input file and the results will be output with the information of phase saturations, pressures, rates and even some visualization files. For the CVFE reservoir simulators used in this research program, the input file is a single Extensible Markup Language {XML) file. The XML is a general-purpose specification for creating custom markup languages. It is classified as an extensible language because it allows users to define their own elements. Its primary purpose of XML is to facilitate the sharing of structured data across different information systems. As shown in Figure 2.12, reservoir simulation input is a typical integrating information file from different resources which makes sense to adopt this format. Other reasons to Figure 2.12 Illustrations of reservoir simulation input and output information 34 choose this format are its well-formed syntax and mature parsing capability supporting from other software packages possibly be used. To demonstrate the input file format, a simple two-phase basement example is shown below: <?xml version="1.0"?> <Simulation name="Example "> <Partition name="Example-part"/> <Fluid name- "PVT-10 "> <oil-density unit= "API"> 45 <</nsity> <water-density unit= "Ibm/ftAJ "> 63 <</density> <oil type="function "> 1130 1 l.e-6 0.90 0.0 <</<water type="function"> 14.7 1 3e-6 0.96 0.0 <</ <</ <Rock name="Rock-10"> <oil-water> 0.00 0.00 1.0000 14.0 0.15 0.01 1.0000 13.0 0.20 0.02 0.9600 7.0 0.30 0.03 0.4300 4.0 0.40 0.04 0.2000 3.0 0.50 0.05 0.1000 2.5 0.60 0.08 0.0600 2.0 0.70 0.125 0.022 1.5 0.80 0.25 0.0100 1.0 <</ter> <compressibility> 14.7 l.e-5 <</ssibility> <</ <PointSource name-"Weill "> <Point> 3874 9669 <</ <PI> 1216 <<//PointSource> <Event time-"0"> <Element relay = "true " action="set thickness "> 35 0 7.0 <</t> <Element relay="true" action="set conductivity" name= "permeability"> 0 F 1000 0 1000 <</t> <ControlVolume relay="true" action = "set ControlVolume"> 0 PoSw 0.3 PVT-10 Rock-10 4300 0.15 <</lVolume> <PointSource relay="false" action-"set PointSource"> Weill open BHP production 200 <</ource> <</ <Event time="0.001 "> <ControlVolume action="write ControlVolume"property="SP"/> <PointSource action="write PointSource"property="Rate Cumulative"/> <</ <Event time="l "> <ControlVolume action = "write ControlVolume"property="SP"/> <PointSource action="write PointSource"property="Rate Cumulative"/> <</ <Event time="10"> <ControlVolume action="write ControlVolume"property- "S P"/> <PointSource action="write PointSource" property^"Rate Cumulative"/> <</ <Event time="100 "> <ControlVolume action="write ControlVolume" property="S P"/> <PointSource action="write PointSource " property= "Rate Cumulative "/> <PointSource relay= "false " action="set PointSource "> <!- example: welll open BHP injection PVT-10 2600 well2 open BHP production 1000 welll open LiquidRate injection PVT-10 26 0.8 well2 open LiquidRate production 2.6 -> Welll open BHP injection PVT-10 3751 <</ource> <</ 36 <Event time="900"> <ControlVolume action = "write Control Volume" property="S P"/> <PointSource action ~ "write PointSource " property^ "Rate Cumulative "/> <</ <Vertex> <!- CVID PNxy z -> 0 2862371.612708016 3796320.6108746817-15305.514648167276 1 2862223.7030747165 3796258.2797346474 -15107.410070856731 2 2862016.6404414456 3 796419.2723770784 -15223.334613773755 3 2861976.595561424 3796214.3124375804 -14872.203719736106 9668 9668 2859619.1288819755 3792828.3183453325 - 11757.863028030839 <</> <Element> <!- ElementlD Elementtype CVIDO CVID1 CVID2 (CVID3) -> 0F0 1 2 1F2 1 3 2 F 4 5 6 3 F 5 7 8 19079 F 241 9668 9666 <</t> . <</tion> . The data structure is clearly shown from the above example basement reservoir simulation input file. Every input file element is strictly defined by well-formed XML syntax with start-tag and end-tag. The above basement feature domain is described by meshed (domain mesh is separately discussed in the next section) 9669 nodes (under the Vertex tag from 0 to 9668) and 19080 triangles (under the Element tag from 0 to 19079). The double-precision Cartesian coordinates of each node is defined after its global ID with the order of x, y, and z. The basement features are composed of meshed triangles. Each triangle element is defined by the node global ID connections. The well in this basement reservoir is defined under the Point Source tag: the first number is the well 37 global ID in the Vertex tag, and the second number is assigned following last number of node global IDs. The well productivity index (PI) is calculated explicitly by Peaceman's (1977) well model: In Equation (2-23), kh is the product of permeability and feature's aperture with the direction of well penetrating fractures, rw is the well radius, and re is the equivalent control volume radius. re could be calculated by Equation (2-24) where A is the control-volume area being occupied by the well point-source. The above example basement domain is a very simplified reservoir model with homogeneous thickness, porosity, permeability, initial saturation and pressures. As shown in the file, some fracture properties such as fracture thickness and transimissivity (kh value) are defined with triangular elements; other fracture properties such as rock type, porosity, initial pressure, initial saturation are defined with the control volumes (node centered). Properties such as density, viscosity relative permeability, capillary pressure, rock compressibility are defined at the beginning of the input file by their own start- and end- tags. These tag names are very straightforward to be found. The output time schedule and well control schedule are under the tag of "Event." For the real simulation cases with multiwells operations, "event" may need to be defined on a day-by-day schedule based on well operation/production data. If this situation happened, the input file might become big. As part of this research program, this situation was solved PI - 2nkh (2-23) (2-24) by some preprocessing interfaces. Those preprocessing interfaces make the input file easy to assemble together. Under the tag of "Fluid" in this example, oil and water's formation factor and viscosity are given by equations with the reference pressure and their compressibility. These equations are (2-25) and (2-26): Si = Brefilexp [~CBil(P - Pr e f)} (2.25) Mi = t^re/,1 exP {P - (2-26) where subscript I denotes phases (oil, water), subscript ref denotes reference status. Example of oil tag , • <oil type="function"> 1130 1 l.e-6 0.90 0.0 <</where reference pressure is 1130 psi, at this pressure, the oil's formation factor is 1.0, the oil formation factor compressibility is 10 6 in this case. Also, at this reference pressure (1130 psi), the oil's viscosity is 0.90 cp and it will not change with pressure (since its viscosity compressibility equals to zero). For the output file from CVFE simulator, based on the "Event" setup on the input file, phase saturations, phase pressures are output by control-volume (node); well injection phase rate and production phase rate are output by well occupied point-source (node or control-volume). Phase saturations and phase pressures can be output as the visualization files which are in the format of vtk. Well injection/production performances are output by phase (oil, water, and gas) rates (STB/day) and cumulative amount to date (STB). Some reservoir behaviors can be analyzed through these well data, such as the oil 38 recovery (ratio of cumulative oil production over original oil in-place), water cut (ratio of water production rate and total fluid production rate), and so on. 2.6 Meshing Domain meshing is an important concept for reservoir simulation. To calculate the fluid flow through the porous media, modeled reservoir domain must be meshed based on numerical discretization method adopted. Today, the trend in simulator development is to separate the meshing part from the flow calculation (simulator). The meshing part is normally done by the meshing software packages. The output from the meshing package is used as the input for the simulator (as shown in the example input file on last section). Through this approach, one simulator may couple with different meshing programs. In a DFN modeled fractured reservoir, the fractures are modeled as two­dimensional surfaces in a three-dimensional space, and the matrix is modeled as three­dimensional solid in three-dimensional space. To get the fluid flow through fractures and matrix, at the beginning, the fractures have to be meshed into two-dimensional triangles. Then the matrix needs to be meshed as tetrahedrons in conformance with the meshed two-dimensional fracture triangles. As mentioned in Section 2.3.3, the bottleneck of simulating DFN characterized fractured reservoir is unstructured domain meshing. Low quality mesh will result in a significant mathematical problem during the flow simulation. Several commercial three­dimensional finite-element mesh software packages were used in this research program. These packages include CUBIT from Sandia National Lab, HyperMesh from Altair Co., RoseMesh from Golder Associate, MeshMaster from Golder Associate. As part of this 39 research work, the above mesh software packages were integrated into the multiphase and multiscale fractured reservoir simulation workflow. Some interfaces were developed among DFN characterization, mesh software and simulators. Furthermore, some mesh postprocess schemes were developed by this research work to analyze mesh quality and even reconstruct/reorder meshes. These efforts normally can improve the reservoir simulation convergence. Taking the example of a basement reservoir, the meshed triangles aspect ratios and control volume ratios usually play an important role during the reservoir simulation. During the reservoir simulation, these mesh shape elements are used in flux calculation, which can directly change the condition number of the Jacobian matrix. Large flux contrast ratio results in a large Jacobian matrix condition number. A problem with a low condition number is said to be well-conditioned, while a problem with a high condition number is said to be ill-conditioned. Normally ill-conditioned matrix results in a convergence problem. Most commercial mesh software offers various mesh quality control options which can maximize assure desired mesh quality being obtained. 2.7 Chapter Summary In this chapter, some fundamentals regarding fractured reservoir simulations were introduced. These fundamentals include basic concepts of black-oil reservoir models, fractured reservoir classifications, flow models to study fractured reservoir, background of the CVFE discretization method and the CVFE simulator development, simulation input/output files, and briefly introduction of unstructured domain geometry meshing. 41 2.8 Nomenclature 'y area, ft formation volume factor control volume of T control volume of node i compressibility, p s i1 ' residual function of node i, partial residual function of bit m total flux flowing out of bt within a triangle, fracture/fault aperture, ft i discretized fracture elements o, w, g, phases . side length matrix block in dual porosity model, ft permeability tensor, md components of the permeability tensor k, md relative permeability number of discretized fracture elements unit outward normal of a boundary fluid pressure, psi capillary pressure, psi •3 •> volume injected or produced, ft /(ft day) mass injected or produced, lbm/(ft3day) gas oil ratio, MSCF/STB 42 r = radius, ft Si = saturation of phase I Siw = connate water saturation Sro = residual oil sasturation T = boundary of polygonal domain T = control volume boundary b, t = time, day tm = triangle m u = Darcy's velocity, ft/day up = upstream function V = volume, ft w = water phase = fluid potential 0 = porosity p = density, lbm/ft3 fi = viscosity, cp 2 a = shape factor, ft' 2.9 Bibliography Aziz, K., Settari, A., "Petroleum Reservoir Simulation," New York: Elsevier (1979) Balay, S., Bushchelman, K., Eijkhout, V., Gropp, W., Kaushik, D., Knepley, M., Mclnnes L.C., Smith, B., Zhang, H., "PETSc Users Mannual," ANL-95/11 - Revision 2.3.3 (2007) . .. 43 Dershowitz, W.S., "Rock Joint Systems," Ph.D Dissertation, Massachusetts Institute of Technology, Cambridge, MA (1984) FracMan Reservoir Edition (FRED) User's Manual, Golder Associates (2007) Fu, Y., Yang, Y.-K., Deo, M., "Three-Dimensional, Three-Phase Discrete-Fracture Reservoir Simulator Based on Control Volume Finite Element (CVFE) Formulation," paper SPE 93292, (2005) Fung, L.S., Hiebert, A.D., and Nghiem, L.X., "Reservoir Simulation With a Control Volume Finite Element Method", paper SPE 21224 (1991) Matthai, S.K., Mezentsev, A., Belayneh, M., "Control-Volume Finite-Element Two- Phase Flow Experiments with Fractured Rock Represented by Unstructured 3D Hybrid Meshes", paper SPE 93341 (2005) Peaceman, D.W., "Fundamentals of Numerical Reservoir Simulation," New York: Elsevier (1977) Ronald A. Nelson, "Geologic Analysis of Naturally Fractured Reservoirs," Gulf Professional Publishing, second edition, (2001) Sukirman, Y.B. and Lewis, R.W.: "Three-Dimensional Fully Coupled Flow: Consolidation Modeling Using Finite Element Method", paper SPE 28755, (1994) Yang, Yi-kun, "Finite-element Multiphase Flow Simulation," Ph.D Dissertation, University of Utah, Salt Lake City, UT (2003) Young, L.C.: "An Efficient Finite Element Method for Reservoir Simulation", paper SPE 7413,(1978) SENSITIVITY STUDY OF OIL PRODUCTION FROM FRACTURED/FAULTED BASEMENT RESERVOIRS (A paper in preparation) Huabing Wang, Craig Forster, Chung-Kan Huang, and Milind Deo Department of Chemical Engineering University of Utah Salt Lake City, UT84112 CHAPTER 3 3.1 Abstract A fully-implicit, three-dimensional (3D), three-phase, discrete fault/fracture, black oil CVFE simulator provides new insight and understanding of oil production from reservoirs in fractured, low-permeability basement rocks. In this chapter, two three­dimensional, naturally fractured reservoirs with different fracture orientations and intensities were characterized by a discrete fracture network (DFN) geological model. Performances of three-phase, black-oil reservoir simulators based on the single porosity finite difference and upstream flux weighted CVFE discretization method were then evaluated on an orthogonal basement fracture network. Benchmark work showed that the results were in close agreement, even under three-phase flow conditions. CVFE is then used to simulate a complex network of intersecting faults that mimic a more realistic basement reservoir. The CVFE simulator provides a method that allows for complex fracture/fault geometries and spatial variations in the internal properties of faults. CVFE simulations of the realistic network illustrate the possible consequences of uncertainty in knowing fracture/fault properties (e.g., porosity, permeability, thickness, dip orientation, connectivity and flow transmissibility). For example, one of the possible consequences is the introduction of spatial variability in permeability within the fault planes (using spatially randomized patterns of 10, 100 and 1000 md), while retaining a constant geometric mean permeability of 100 md. This enhanced oil production is due to the high-permeability pathways. A 50:50 mix of 10 and 1000 md elements yielded 36% OOIP recovery while a 33:33:33 mix of 10, 100 and 1000 md yielded 24% OOIP production. These results were 25% and 13% greater, respectively, than those obtained by the uniform 100 md case (1 l%OOIP). This inherent variability, combined with the 45 uncertainty of knowing the detailed connections between different regions of the fault network, impacts the pattern of sweep in ways that are important to understand when attempting to develop innovative approaches to reservoir management. 3.2 Introduction Quantifying fractured reservoirs is much more difficult than quantifying nonfractured reservoirs due to the extreme complexity of the fractured reservoir. The complexity originates from the vast number of uncertain variables that dictate the final fractured reservoir response. Examples of these variables and their relatives in a fractured basement reservoir include the reservoir's storage from fracture thickness and porosity; fracture permeability/transimissivity distribution characterizations; and fracture geometry such as dip orientation, connectivity, area sizes etc. The single or combined effects of those fracture variables mentioned here can make the basement network act as both permeable pathways to fluid migration and also significant storage. This brings a wide range of uncertainties for basement fractured reservoir management. At the beginning of this chapter, a regularized fracture network was constructed that was amenable to gridding using ECLIPSE, the conventional finite difference reservoir simulator. A three-dimensional, three-phase, black oil modeled control-volume finite element (CVFE) simulator was used in this study. In order to validate the CVFE simulator, the results of regularized DFN direct simulations were compared with the results from ECLIPSE. Once the accuracy of the DFN approach was ascertained, some sensitivity of "real" basement fracture network with the multiphase DFN modeled basement reservoir was directly simulated through the CVFE simulator. 46 Permeability distribution in the fracture/fault zone, dip orientation and thickness variations will be considered in this sensitivity study. A typical fracture network is highly heterogeneous. To understand the impact of permeability heterogeneity between fractures, within fractures and the presence of preferred high permeability pathways between injectors and producers were quantitatively studied in this conceptual basement fractured network. The extent to which the preferred pathway network leads to early breakthroughs was also quantified by this chapter. The impact of shutting off the high permeability fracture using gel treatments on improving oil recovery will also be explored by this chapter. 3.3 Governing Equations and DFN Model As shown in Equation (3-1) to (3-3), the equations describing compressible three-phase flow are obtained by combining Darcy's law and mass conservation for each phase: OIL: - V . u o = £ ( 0 ^ ) + Qo (3-1) WATER: - V ■ ulv = £ + qw (3-2) GAS: -V ■ (Rsu 0 + ug ) = - ^0 + 0 + Rsq a + qfg (3-3) where subscript o, w, and g represent oil, water and gas phases, q is the source term, Rs is the gas-oil ratio, S is the saturation, B is the formation volume factor. The LHS in the above equations represent the flux term, and u denotes the Darcy phase velocity (by combining Equation (2-5) and (2-6)). 47 In this three-phase black-oil model, capillary pressures coupling phase pressures are listed in Equation (3-4) and (3-5) as = P. - (3-4) P.=(5J = - po (3-5) and the volume conservation equation is defined in Equation (3-6) S0 + S w + S g = 1 (3-6) where the subscripts o, w and g refer to oil, water and gas phases, u, p and p are the fluid velocity, density and viscosity. The rock is characterized by the porosity 0 and absolute permeability k. The multiphase properties are the relative permeability kr and capillary p c. Rs represents the gas-oil ratio which is the solubility of gas in oil as a function of pressure. B is the formation volume factor. Source terms are designated by q. The equations for three-phase flow are derived from the continuity equation: the left-hand side (LHS) is the flux term and the right-hand side (RHS) represents the cumulative term and the source term. For the basement reservoir study, all of the fluid flow through fractures is assumed to be governed by Darcy's law. The discretization of these equations for modeling flow in basement reservoirs entails the use of the DFN approach and the CVFE scheme. In the DFN modeled three­dimensional basement domain, all of the single fractures are represented by two­dimensional surfaces and those two-dimensional surfaces will be meshed into triangle elements. These dimensional reductions greatly decrease the node numbers, resulting in reduced computational time. A control volume is usually distributed across several triangular elements and calculated by summing together its subvolumes, as shown in Figure 3.1. 48 49 Figure 3.1 A control volume with its boundaries (dashed line) across several triangle elements and decomposition of a control volume into several subvolumes Figure 3.1 shows an example of a control volume with its boundaries (dashed line) across several triangle elements and decomposition of a control volume into several subvolumes. Detailed background of CVFE discretization formulation and simulator development refer to Section 2.5, the fundamentals of fractured reservoir modeling. 3.4 Proposed Domain Characteristics 3.4.1 Regularized DFN Basement Domain At the beginning of this study, a regularized well connected fracture network domain was generated for the purpose of verifying the CVFE simulator with the most popular commercial black oil simulator ECLIPSE. Since all of the features were connected orthogonally, the finite-difference based ECLIPSE simulator gave the most accurate solution. Two constraints were applied for the DFN regularized domain and its equivalent ECLIPSE single-porosity domain: the same OOIP and the same transmissivity (kh value). The regularized domain was characterized as a three-phase black oil model with initial solution gas in the oil phase. All simulations have the exact same fluid and rock properties. Examples of these properties include fluid density, viscosity, compressibility, PVT table, formation factor, gas/oil ratio table, rock compressibilities, relative permeability curves, initial conditions, well operation conditions and schedules, etc. Common model properties in the basement model are shown in Figure 3.2: • Impermeable matrix with 0 = 0 (Type I, basement reservoir system) • Domain = 1,000 ft by 1,000 ft by 200 feet deep • Total feature length = 30,000 feet • Feature property: k = 1,000 md, 0 = 1 4 %, width = 0.5 feet • Original Oil In Place (OOIP) = 53,580 STB • Injection Pressure = 4,300 psi 50 Figure 3.2 Regularized DFN basement model (left) and equivalent regularized ECLIPSE single porosity model 3.4.2 Irregularized (Hypothetical "Real") DFN Basement Domain As shown in Figure 3.3, the properties in the irregularized DFN basement model were exactly the same as the regularized domain except for the feature orientation and connectivity. To distinguish the regular "ideal" basement model discussed in last section, this irregular domain was named the hypothetical "real" basement model. Two questions to be answered were: Can this hypothetical "real system" be easily represented by the equivalent "ideal system"? How could this domain be best represented by the traditional finite-difference based reservoir simulators such as ECLIPSE? 51 Figure 3.3 Irregularized DFN basement model with the same model properties as regularized DFN basement model but different orientation and connectivity 3.4.3 Geological Characteristics Comparisons Between Hypothetical "Real" System and Regular "Ideal" System Certain fracture network geological characterization analyses were done on both "ideal" and "real" systems. Figure 3.4 shows that the regularized "ideal" system has only two orientations as (90, 0) and (180, 0); for the hypothetical "real" system, the main fracture pole orientations are distributed from (130, 0) to (350, 0). Figure 3.5 shows that the regularized "ideal" system has the uniform feature equivalent length (radius) of 252 feet, but the fracture sizes of hypothetical "real" system offered power law fitting of xrnin = 252.1 feet and exponential term b = 5.2. Detailed fracture network analysis is shown in Table 3.1. Fracture cluster analyses were performed on both "ideal" and "real" models. No isolated fracture was found in either system. The average fracture center was a little bit different: 1) Regularized "ideal system": (500.00, 500.00, 100) 2) Hypothetical "real system": (513.06, 488.70, 100) 3.5 CVFE Simulator Verification The presence of fractures makes the reservoir domain geometrically complicated, particularly when the issue of connectivity dominates recovery. One of the key advantages of the DFN model, compared with conventionally fractured modeling, is to present the fracture orientation and connectivity explicitly. Obviously, the traditional finite-difference based simulator offered the option of the dual porosity model to simplify connectivity and heterogeneity issues that are modeled much more accurately by the DFN 52 53 o Figure 3.4 Fracture pole orientation (Wulff equal-angle projection, lower hemisphere) between hypothetical "real system" and "ideal system" Trace Data Power Law Plot Figure 3.5 Comparisons of fracture size distribution (Power Law Plot) between the hypothetical "real system" and "ideal system" 54 Table 3.1 Fracture set statistics Regularized "Ideal System" Hypothetical "Real System" Number of fractures 30 107 Total fracture area 6000000 6000000 Total fracture volume 3000000 3000000 P32 (fracArea/volume) 0.02 0.02 P33 (fracVolume/volume) 0.01 0.01 Mean orientation 135,0 197.2836, 0 Equivalent radius mean 252.13 130.59 Equivalent radius std dev 1.06 28.34 Equivalent radius min 252.13 83.40 Equivalent radius max 252.13 237.99 Area mean 200000 56075 Area std dev 0 26292 Area min 200000 21850 Area max 200000 177946 Thickness, mean 0.5 0.5 Thickness, std dev 0 0 Thickness, min 0.5 0.5 Thickness, max 0.5 0.5 Permeability, mean 1000 1000 Permeability, std dev 0 0 Permeability, min 1000 1000 Permeability, max 1000 1000 Compressibility, mean 5e-006 5e-006 Compressibility, std dev 0 0 Compressibility, min 5e-006 5e-006 Compressibility, max 5e-006 5e-006 approach. The characterization by the DFN model requires rotatable permeability tensor along fractures/faults, which the finite-element discretization formula achieves especially well. A three-dimensional, three-phase black oil reservoir simulator was developed using a new CVFE formulation (Yang, 2003; Fu et al., 2005; Fu, 2007). In this new CVFE approach, flux-based upstream weighting was applied to ensure flux continuity both locally and globally. This weighting function solved the local mass conservation problem from traditional finite-element or CVFE simulations. Like all problems in the field of computation science and engineering, a new product needs to be "validated" and "verified" before it is actually implemented. The term validation means solving the right equation for a given problem and verification means solving the equations correctly. In the area of reservoir simulation, the same governing equations for a black-oil model have been widely used for the past few decades. Therefore, validation of a black-oil model is not considered for most new black-oil simulators. At the beginning, the focus of this study was to benchmark the CVFE simulator with the well-established commercial finite-difference based reservoir simulator ECLIPSE. This work could also be called as simulator verification by indexing method. On this verification work, the regularized "ideal system" was tested with 15 by 15 feature/faults strictly orthogonal in the domain of 1000 by 1000 by 200 ft. The matrix was considered to contribute zero porosity and permeability. The detailed reservoir properties are described in Chapter 3.4. This is a three-phase (oil, water and gas) study. Initially the reservoir pressure was over the bubble point which ensures all of the gas solute in the oil phase. One injection well and one production well were positioned at the diagonal direction of the reservoir domain as Figure 3.2 showed. Both wells horizontally penetrated the top of features. Bottom-hole pressure (BHP) controls were applied for well operation. Since the DFN approach presents feature/faults as a two-dimensional plane in the three-dimensional space, the domain was meshed into triangles which shared the same side or node at the region of feature/faults interaction. For ECLIPSE simulator, the features were presented by the stack of grid blocks which have the actual width of the features (0.5 ft) or equivalent thickness by the constraint of constant kh value (500 md-ft). To speed up calculations, the matrix blocks other than features in the ECLIPSE were set as nonactive. Another constraint for both the "ideal" and "real" systems was to set the constant total original oil in place (OOIP) as 53580 STB. The discretized domain for the DFN approach and conventional finite-difference are shown in Figure 3.6 and Figure 3.7. Figure 3.6 and Figure 3.7 show that both the DFN triangular mesh and the ECLIPSE grid blocks have five element/grid layers in the z-direction. The simulations were run for 900 days with the first 300 days as the primary production period. Starting from day 301, simulation was considered to be secondary recovery with the water injector being introduced. 56 y Figure 3.6 Three-dimensional view (left) and two-dimensional view (right) of discretized regularized, DFN modeled "ideal" basement reservoir domain for the CVFE simulator 57 ---------------------------------------------------* x or y Figure 3.7 Three-dimensional view (left) and two-dimensional view (right) of discretized, regularized finite-difference modeled "ideal" basement reservoir domain for the ECLIPSE simulator Simulation results from CVFE and ECLIPSE simulators such as oil production rates and oil recovery are compared in Figures 3.8 and 3.9. Figure 3.8 and Figure 3.9 show that results of both oil production rate and oil recovery from CVFE simulation are almost the same as the results from ECLIPSE. Since the regularized basement model is strictly orthogonal (which can be accurately calculated by ECLIPSE), this close agreement between results made a very positive conclusion about the CVFE simulator. Note that this case study was performed under the three-phase condition, and with the solution gas involved, the simulation results showed that CVFE simulator was benchmarked with ECLIPSE black-oil simulator remarkably well. This conclusion can also be observed from oil distribution snapshots shown in Figure 3.10. 58 Time (days) Figure 3.8 Comparisons of oil production rate between the DFN modeled CVFE simulator and the single porosity modeled ECLIPSE simulator for the regularized "ideal system" Tiine(days) Figure 3.9 Comparisons of oil recovery between the DFN modeled CVFE simulator and the single porosity modeled ECLIPSE simulator for the regularized "ideal system" 59 CVFE Simulator ECLIPSE Simulator Dav I) Day 300 End Pr imary Day 900 End SimuhUio Figure 3.10 Snapshot of oil saturations for the regularized "ideal" basement domain from both CVFE simulator (left column) and ECLIPSE simulator (right column) 60 The verification of constraints (constant kh value and constant OOIP) are shown in Figure 3.11. At day 900 (end of simulation), simulations with 0.5 ft to 20 ft fracture thickness had reached the same the residual oil distributions from the ECLIPSE simulator. The indexing verification methodology adopted by this study was commonly applied to verify a new reservoir simulator in the petroleum industry. Other more strict mathematical verification work has been done (Yang, 2003; Fu, 2007) on the CVFE formulation by implementing the manufactured solution method. In the manufactured solution approach, a solution function that satisfied both initial and boundary conditions was first synthesized or manufactured. That solution function was then substituted into the governing equation to find the corresponding source function. Mathematically, the solution obtained from the CVFE simulator showed an excellent order of convergence rate with the manufactured solution. Figure 3.11 Comparisons between feature/fault grids block thickness (left 20ft; right 0.5ft) by applying two constraints: OOIP & kh value. Through this indexing verification with the well-known ECLIPSE simulator, now the accuracy of the CVFE simulator is ascertained. Some sensitivities of DFN characteristics are discussed below. 3.6 Applications Uncertainty exists around the fractured reservoir due to the fracture network's heterogeneous properties. These fracture network properties include connectivity, permeability distribution, dip orientation, and size distribution. To demonstrate the impact effects of uncertain elements on the multiphase oil recovery scheme, as shown in Figure 3.12, this study used a "real" fracture network generated by a network of intersecting faults that mimic a more realistic basement reservoir. 3.6.1 Permeability Variations Studies Permeability is an important property of the porous medium and is a measure of the capacity of the medium to transmit fluids. Hydrocarbon reservoirs can have primary and secondary permeability. 61 Figure 3.12 Irregularized DFN basement model with the same model properties as regularized DFN basement model with different orientation and connectivity 62 Primary permeability is referred to as matrix permeability by reservoir engineers. Secondary permeability can be by either fractures or solution vugs. In this chapter, only the Type I (basement type) reservoirs are considered. Therefore, all the permeability discussed in this chapter is fracture permeability. Generally fracture permeability abides by cubit law of fracture thickness. Since the fracture is generated due to the stress and strain acting on the rock, the thickness of fractures cannot be homogeneous. This results in a highly heterogeneous permeability distribution within each fracture/fault. In this section, the hydrocarbon recovery scheme from three homogenous permeability models (1000 md, 100 md, 10 md) and two heterogeneous permeability models (with random geometric mean of 100 md) have been compared in order to demonstrate the fracture network uncertainty. The homogenous permeability models are easily understood. The results of oil saturation distribution snapshots from various homogeneous permeability models are shown in Figure 3.13. Figure 3.13 At day 900, oil saturation distribution snapshots from different homogenous permeability models 63 Figure 3.14 demonstrates that with lowering the homogeneous permeability by one or two magnitudes in the whole basement reservoir, water breakthrough times were delayed significantly. The simulated water breakthrough time, after the water injections, is 185 (1000 md), 3200 (100 md), and 14700 (10 md). The calculation geometric mean of a permeability data set [ki, k2, ..., kn] could be described by the Equation (3-7): dir=ikf)1/" = Vk i -k j ......k„ (3-7) where k is the permeability, and n represents the numbers of discretized elements. Mathematically, geometric mean is an average indicating the central tendency or typical value of a set of number. It can be understood in terms of geometry. The 1 -i 0.9 - 0.8 - 0.7 - * 3*1 0.6 - ter ( 0 iji 1 £ 0.4 - I1 « - * 0.3 - /i 0.2 - f0 |1 0.1 ■ 1«r 0 - i / 1000 md Uniform 100 md Uniform 5000 + 10000 Time (days) 15000 20000 Figure 3.14 Comparisons of water breakthrough and water cut among different homogeneous permeability models geometric mean of two numbers, a and b, is simply the side length of the square whose area is equal to that of a rectangle with side lengths a and b. that is what is n such that n2=axb. Similarly, the geometric mean of three numbers, a, b, and c, is the side length of a cube whose volume is the same as that of a rectangular prism with side lengths equal to the three given numbers. This geometric interpretation of the mean is probably what gave it its name. It is important to note that the geometric mean applies only to positive values to prevent the result of imaginary numbers. For example, no negative fracture permeability exists physically. However, physically there might be some "zero" permeability existing inside the fault. In that case, other average methods might be applied instead of the geometric mean. Figure 3.15 illustrates the geometric mean based intrafault permeability distributions. If the intrafault were composed of two types of permeability values: 10 and 64 50:50 (10, 1000 md) Geom. Mean = 100 md Figure 3.15 Illustrations of geometric mean based intrafault permeability distributions 65 1000 md, the 100 md geometric mean will give the half-half percentage of the two permeability values. Then the permeability of discretized elements can be characterized based on the percentage distribution. If the intrafault were composed of three types of permeability (10, 100 and 1000 md), the geometric mean of 100 md will result in one-third quantity distributions with each permeability value. Figure 3.16 shows residual oil saturation of three basement fracture networks (one uniform 100 md case and two various lOOmd geometric mean domains). As previously shown in the simulation, all three cases passed primary production stage (first 300 days) and secondary production (from day 301 to 900). Obviously, the water flooding patterns are totally different by different fracture network characterizations: uniform lOOmd case study offers the medium performance between the two types of 100 md geometric mean permeability values, the one-third permeability distribution cases result in the best water flooding and the half-half case offers the worst water flooding scenario. Gcom. Mean 100 md Gcom. Mean 100 md 100md 50:50 Rondomk 33:33:33 Rondom k Figure 3.16 At day 900, oil saturation distribution snapshots from different homogenous permeability models Figure 3.17 to 3.19 compared simulation results as the oil production rate, oil recovery and water cut. Obviously, the randomized intrafault permeability distributions brought totally different hydrocarbon recovery mechanisms to our proposed basement domain. The figures of the oil production rate, water breakthrough time, and water cut scales show possible large-range variations. In the random permeability distribution case, the oil production rate could be higher than in the uniform permeability case, due to possible high permeability channels (with those 1000 md characterized elements); the oil production rate also could be lower than in the uniform permeability case due to the low permeability in the main water flooding channel (with 10 md characterized elements). However, if a certain number of simulations were performed based on different random seeds, some uncertainty range might be identified. The more simulations performed, the more accurately the risk range might be described. Based on this hypothesis, a set of 66 E-i w o £- Time (days) Figure 3.17 Oil production rate comparisons among uniform 100 md and geometric mean 100 md cases 67 0 2000 4000 6000 8000 Time (days) Figure 3.18 Cumulative oil production comparisons among uniform 100 md and geometric mean 100 md cases Time (days) Figure 3.19 Water cut comparisons among uniform 100 md and geometric mean 100 md cases 68 geometric mean (33:33:33) 100 md simulations have been performed. The only difference among these simulations is the detailed element permeability distributions on each single feature. The water cut results of these simulations are shown in Figure 3.20. ' In the water flood case studies, water cut and water breakthrough time are the principal ways of making reservoir management decisions. Figure 3.20 shows a large range of water breakthrough times. For the same domain with the exact same operating conditions, the water flooding patterns were totally different. This type o f modeling mimics big impacts on the pattern of sweep, due to the uncertainty of the fracture/fault permeability distributions. In a way, this is important for understanding the fracture network uncertainty when attempting to develop innovative approaches to reservoir management. 1 0.9 0.8 0.7 b 0.5 0.4 0.3 0.2 0.1 0 0 2000 4000 6000 8000 10000 Time (days) Figure 3.20 Water cut comparisons for 33_33_33 geometric mean with random seeds Some conclusions can be drawn from the studies in this section: 1) The rock permeability is an important parameter of fluid flows inside the fracture networks. The variation of rock permeability distributions results in a huge uncertainty, particularly for the secondary production with water flooding. 2) Introducing spatial variability in permeability within the fault planes (using spatially randomized patterns of 10, 100 and 1000 md), while retaining a constant geometric mean permeability of 100 md, yields enhanced oil production due to the high-permeability pathways. And these kind of pathways result in different oil recovery factors for the long-term or short­term view, which might suggest totally different recovery plans for reservoir management. 3.6.2 Fracture Network Connectivity Studies The connectivity of fracture/fault in subsurface rock formations is a key factor in understanding and predicting fluid flow in hydrocarbon reservoirs. This is especially true in Type I basement reservoirs, since in this type of reservoir all hydrocarbons exist inside the fractures, and the fractures provide essential permeability. This permeability then makes the connectivity problem fatal for the hydrocarbon production. Some impacts, due to the connectivity, were already shown in the simulations in the above two sections. An example of residual oil resulting from the fracture network's connectivity is shown in Figure 3.21. In this figure, some oil traps (circled) result from the fracture connectivity. Obviously these main traps are limited in connection with the major water flooding way between the injector and the producer. In this section, the critical fracture/fault connections between the injector and the producer are identified by the trace map in Figure 3.22. Figure 3.22 demonstrates that the fracture networks were divided into two separate zones: west zone (left-side) and east zone (right-side). These two zones are connected by four east-west fractures/faults. The injector and the producer are located in different zones. The whole reservoir is assigned with uniform permeability at 1000 md. 70 Figure 3.21 Residual oil saturation distribution snapshot due to the fracture network connectivity (Case 1000 md uniform k with day 600) 71 Figure 3.22 Illustration of two separate fracture zones (separated by the acclivitous line in the middle of trace map) connected by four fractures (circles) The hypothesis for studying the fracture connectivity and how it impacts oil production was produced by comparing the productions from the original fracture network a