Advanced techniques for reservoir engineering and simulation

The outstanding surge in hydrocarbon production from unconventional reservoirs is unprecedented. Profitable oil prices and new technologies have untapped massive oil and gas resources in recent years. However, the correct exploitation of these resources has been dampened by the lack of understanding of these systems. Research efforts to understand and properly assess unconventional resources have exploded in the literature. In this research work, a series of advancements in reservoir production analysis, simulation modeling, and simulation development are made. A semi-analytical method based on conventional material balance was developed to approximate reservoir pressure distributions and permeability. One of the strengths of this method is that it only requires limited information to be viable. Reservoirs with dry gas and/or high gas oil ratios are handled with an additional average pressure correction factor that takes gas compressibility into account. Hence, this method can be used for any type of fluid and fluid flow as long as the correct material balance formulation and surrogate curves are employed. Verification of the method is made through comparison with synthetic data and a field case study. Furthermore, a standardized simplification workflow for hydraulically stimulated reservoirs was introduced. The aim of this workflow is to guide the engineer when developing a simplified reservoir simulation model with multiple wells and fractures. Simplified models have been around for a long time in the literature, however, their applicability to field-scale projects is very limited. Models that result from the application of this workflow are shown to retain the low simulation run-times characteristic of popular single-fracture models. In addition, fluid rate results from the proposed workflow models are in good agreement with results from full-scale simulation models. This is not the case for the single-fracture model which loses accuracy as the complexity of the project grows. Lastly, a new discrete fracture model formulation is implemented in a control-volume finite element simulator. This new fracture model provides fractures with their own control volumes and gives them freedom to be placed anywhere in the matrix domain. Verification of this implementation is made through comparison with analytical expressions and other well-established simulators.

ADVANCED TECHNIQUES FOR RESERVOIR ENGINEERING AND SIMULATION by Raul Velasco Guachalla A dissertation submitted to the faculty o f The University o f Utah in partial fulfillment o f the requirements for the degree o f Doctor o f Philosophy Department o f Chemical Engineering The University o f Utah May 2016 Copyright © Raul Velasco Guachalla 2016 All Rights Reserved The Unive r si t y of Utah Graduat e School STATEMENT OF THESIS APPROVAL The thesis o f ____________________ Raul Velasco Guachalla___________________ has been approved by the following supervisory committee members: ____________________ Milind Deo_____________________ , Chair 10/26/2015 Date Approved ______________ John David McLennan_______________ , Member 10/26/2015 __________________ Philip J. Smith___________________ , Member 11/23/2015 _______________ Anil Vasudeo Virkar________________ , Member 10/26/2015 Date Approved ________________Rasoul B. Sorkhabi_________________ , Member 10/26/2015 Date Approved and by ________________________ Milind Deo_________________________ , Chair/Dean o f the Department/College/School o f _______________ Chemical Engineering______________ and by David B. Kieda, Dean o f The Graduate School. ABSTRACT The outstanding surge in hydrocarbon production from unconventional reservoirs is unprecedented. Profitable oil prices and new technologies have untapped massive oil and gas resources in recent years. However, the correct exploitation o f these resources has been dampened by the lack o f understanding o f these systems. Research efforts to understand and properly assess unconventional resources have exploded in the literature. In this research work, a series o f advancements in reservoir production analysis, simulation modeling, and simulation development are made. A semi-analytical method based on conventional material balance was developed to approximate reservoir pressure distributions and permeability. One o f the strengths o f this method is that it only requires limited information to be viable. Reservoirs with dry gas and/or high gas oil ratios are handled with an additional average pressure correction factor that takes gas compressibility into account. Hence, this method can be used for any type o f fluid and fluid flow as long as the correct material balance formulation and surrogate curves are employed. Verification o f the method is made through comparison with synthetic data and a field case study. Furthermore, a standardized simplification workflow for hydraulically stimulated reservoirs was introduced. The aim o f this workflow is to guide the engineer when developing a simplified reservoir simulation model with multiple wells and fractures. Simplified models have been around for a long time in the literature, however, their applicability to field-scale projects is very limited. Models that result from the application o f this workflow are shown to retain the low simulation run-times characteristic o f popular single-fracture models. In addition, fluid rate results from the proposed workflow models are in good agreement with results from full-scale simulation models. This is not the case for the single-fracture model which loses accuracy as the complexity o f the project grows. Lastly, a new discrete fracture model formulation is implemented in a control-volume finite element simulator. This new fracture model provides fractures with their own control volumes and gives them freedom to be placed anywhere in the matrix domain. Verification o f this implementation is made through comparison with analytical expressions and other well-established simulators. iv To my family. TABLE OF CONTENTS ABSTRACT............................................................................................................................................ iii LIST OF TABLES................................................................................................................................ ix LIST OF FIGURES................................................................................................................................ x NOMENCLATURE...........................................................................................................................xiv ACKNOWLEDGEMENTS................................................................................................................xx Chapters 1. INTRODUCTION.............................................................................................................................1 1.1 Unconventional Hydrocarbon Production.......................................................................1 1.2 Research Motivation............................................................................................................. 3 1.3 Objectives.................................................................................................................................5 2. MATERIAL BALANCE APPLIED TO TIGHT FORMATIONS..................................... 7 2.1 Background.............................................................................................................................8 2.2 General Procedure ............................................................................................................... 10 2.2.1 Corrected Average Pressure..................................................................................... 12 2.3 Case Examples ...................................................................................................................... 15 2.3.1 Example 1: Multiphase Linear Flow into a Vertical Fracture.......................15 2.3.2 Example 2: Multiphase Linear Flow into a Vertical Fracture and High GOR............................................................................................................................... 19 2.3.3 Example 2 (Revisited): Multiphase Linear Flow into a Vertical Fracture and High GOR..................................................................................................................... 20 2.3.4 Example 3: Gas Linear Flow into a Fracture......................................................23 2.3.5 Example 4: Permian Basin Field Case..................................................................26 2.4 Other Applications.............................................................................................................. 29 2.5 Key Findings........................................................................................................................32 2.6 Method Development......................................................................................................... 33 2.6.1 Material Balance Theory........................................................................................... 33 2.6.2 Corrected Average Pressure Theory ...................................................................... 36 3. RESERVOIR SIMULATION SIMPLIFICATION WORKFLOW FOR HYDRAULICALLY FRACTURED RESERVOIRS..........................................................40 3.1 Background.......................................................................................................................... 41 3.2 Workflow Components......................................................................................................44 3.2.1 Interference Effects ..................................................................................................... 44 3.2.2 Boundary e ffe c ts......................................................................................................... 45 3.2.3 Multiple Hydraulic Fracture Representation......................................................46 3.2.4 Other Simulation Phenomena..................................................................................47 3.2.5 Symmetry Application...............................................................................................48 3.3 Proposed Simplification Workflow................................................................................ 49 3.4 Case Examples......................................................................................................................51 3.4.1 Single Well Case......................................................................................................... 51 3.4.2 Multiwell and Multifracture Case .......................................................................... 55 3.5 Results....................................................................................................................................58 3.5.1 Single Well Case ......................................................................................................... 58 3.5.2 Multiwell and Multifracture Case .......................................................................... 64 3.6 Fracture and Well Spacing Application........................................................................68 3.7 Key Findings........................................................................................................................76 4. A NEW DISCRETE FRACTURE MODEL IMPLEMENTATION................................ 77 4.1 Background.......................................................................................................................... 78 4.1.1 Natural and Hydraulic Fractures.............................................................................78 4.1.2 Fracture Representation in Numerical Simulation............................................ 80 4.2 Governing Equations...........................................................................................................88 4.3 Numerical Methods ............................................................................................................. 90 4.3.1 Transmissibility............................................................................................................91 4.3.2 Mobility Term Upstream Weighting.................................................................... 98 4.3.3 Formulation o f Residual Functions........................................................................ 99 4.3.4 Current CVFEM-Based DFM Implementation................................................. 99 4.3.5 New CVFEM-Based DFM Implementation.....................................................102 4.3.6 Well M od e l.................................................................................................................111 4.4 Verification Studies...........................................................................................................114 4.4.1 Finite-Conductivity Fracture Behavior.............................................................. 115 4.4.2 Multistage and Multiwell Problem..................................................................... 124 4.5 Key Findings.......................................................................................................................124 5. SUMMARY.................................................................................................................................... 127 5.1 Original Contributions......................................................................................................127 5.2 Recommendations for Future Work.............................................................................128 Appendices A. RESERVOIR ROCK AND FLUID PROPERTIES...........................................................130 B. CONVENTIONAL MATERIAL BALANCE..................................................................... 134 vii C. GOODNESS OF F IT ................................................................................................................... 136 D. ECONOMIC ANALYSIS...........................................................................................................138 REFERENCES....................................................................................................................................139 viii LIST OF TABLES 2.1 Reservoir and operational parameters for simulations........................................................ 16 3.1 Summary o f reservoir model and operational parameters................................................. 52 3.2 Symmetric and multiplier factors for elements 1 and 2 for the "Proposed model 1" . 59 3.3 Comparison o f simplified models to a full-scale model through statistical analysis for the single well c a s e ............................................................................................................................... 63 3.4 Run time comparison for different models for the single well c a s e ............................... 64 3.5 Symmetric and multiplier factors for a multiple horizontal well case........................... 65 3.6 Comparison o f simplified models to a full-scale model through statistical analysis for the multiwell c a s e ................................................................................................................................ 67 3.7 Run time comparison for different models for the multiwell c a s e ..................................68 3.8 Spacing study reservoir simulation data..................................................................................71 3.9 Economic analysis capital and operating co sts......................................................................72 4.1 Simulation key parameters........................................................................................................ 121 LIST OF FIGURES 1.1 Map o f U.S. most prominent unconventional plays (Source: EIA )...................................2 1.2 U.S. tight oil (left) and dry shale gas (right) production (Source: E IA ).......................... 3 2.1 Original Oil In Place calculations for two systems showcasing the "Distance" variable..................................................................................................................................................... 12 2.2 Difference between equilibrium average pressure (material balance pressure) and "real time" average pressure........................................................................................................................ 14 2.3 Linear flow into a vertical fracture........................................................................................... 15 2.4 Average pressure profile matched by a surrogate curve.....................................................17 2.5 Calculated material balance pressure profiles versus simulation....................................18 2.6 Calculated material balance pressure profiles versus simulation....................................20 2.7 Pressure profiles determined using corrected average pressures..................................... 23 2.8 Pressure profiles determined using the material balance method....................................25 2.9 Pressure profiles determined using corrected average pressures..................................... 26 2.10 Transient material balance method application to a field c a s e ...................................... 28 2.11 Field case calculated average pressure decline................................................................... 29 2.12 Calculated material balance permeability for an infinite reservoir............................... 30 2.13 Pressure decline over time at a point 50 ft away from a draining fracture.................31 2.14 Calculated oil saturation profile...............................................................................................32 2.15 Typical error function pressure profile..................................................................................34 2.16 Average pressure profile plot calculated by varying L .....................................................35 2.17 Material balance average pressure compared to an average pressure surrogate curve.......................................................................................................................................................... 36 2.18 System consisting o f tanks 1 and 2 .........................................................................................37 2.19 Equivalence between Tanks 1 and 2 at producing state and Tanks 1 and 2 at equilibrium state.................................................................................................................................... 38 3.1 Aerial view o f typical multiwell and multistage reservoir simulation model and how it may be simplified..................................................................................................................................43 3.2 Types o f internal and external fractures and their corresponding volume assignments (top v iew )................................................................................................................................................ 46 3.3 Example o f a solution process diagram that shows various force contributions to flow in a simple hydraulically fractured reservoir................................................................................. 50 3.4 Horizontal well with 50 equally-spaced fractures............................................................... 51 3.5 Aerial view o f a full-scale horizontal well and its simplification approaches as presented by the single fracture and proposed models. Dashed lines represent no-flow boundaries for each model...................................................................................................................53 3.6 Solution process diagram simplification for a single well model (a) and multi-well model (b)................................................................................................................................................. 54 3.7 Full-scale multiwell simulation model schematic.................................................................56 3.8 Aerial view o f a full-scale multi-well model and its simplification approaches as presented by the single fracture and proposed models. Dashed lines represent no-flow boundaries for each model...................................................................................................................57 3.9 Cumulative oil comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b )..........60 3.10 Oil rate comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b )......................60 3.11 Internal and external fracture cumulative GOR comparison for a 50 nD case (a) and a 1000 nD case (b )...................................................................................................................................61 3.12 Cumulative GOR comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b) ... 62 3.13 Cumulative oil production comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b ).............................................................................................................................................65 3.14 Oil rate comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b )......................66 3.15 Cumulative GOR comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b) ... 66 3.16 Net present value for a 50 nD matrix permeability spacing study................................ 73 3.17 Net present value for a 500 nD matrix permeability spacing study..............................74 xi 3.18 Net present value for a 1000 nD matrix permeability spacing study........................... 75 4.1 Fault and joint visual representation.........................................................................................79 4.2 Schematics o f dip slip and strike slip faults............................................................................79 4.3 Top view o f a single porosity fracture model with linear grid refinement................... 82 4.4 Dual porosity model representation..........................................................................................83 4.5 Matrix-Fracture connectivity schematic for the (a) dual porosity model, (b) subdomain model and (c) dual porosity / dual permeability model..............................................................86 4.6 Reservoir triangular element discretization around a fracture (red lin e ).......................87 4.7 Transmissibility illustration for the classic finite difference method.............................93 4.8 Example configurations for a point-distributed (left) and block-centered grid (right) systems ..................................................................................................................................................... 93 4.9 Transmissibility illustration for the corner-point method.................................................. 95 4.10 Simple triangular element and its associated subcontrol volumes................................ 96 4.11 Assembly o f control volume associated to vertex " 1" ......................................................96 4.12 Example discrete fracture representations......................................................................... 101 4.13 Discrete fracture elements for the (a) current and (b) the suggested implementation...................................................................................................................................103 4.14 Discrete fracture (in red) fragmentation for (a) a simple and (b) a complex network................................................................................................................................................. 104 4.15 Exaggerated fracture block representations for (a) regular and (b) irregular geometries.............................................................................................................................................105 4.16 Fracture representation in linearized systems................................................................... 106 4.17 Fracture grid blocks surrounded by matrix sp a c e s..........................................................109 4.18 Representation o f matrix nodes and fracture blocks interactions............................... 109 4.19 Flux from a matrix node onto a fracture grid block........................................................ 111 4.20 Well model representation in a simulation m od e l...........................................................112 4.21 Well representation for a Cartesian grid (left) and triangular-element grid (right) 113 xii 4.22 Finite-conductivity fracture sy stem ..................................................................................... 116 4.23 Fracture fluid flow m od e l....................................................................................................... 117 4.24 Analytical and simulation pressure drop results for a 10-segment fracture.............122 4.25 Analytical and simulation pressure drop results for a 20-segment fracture.............122 4.26 Areal simulation pressure distribution for the sy stem ................................................... 123 4.27 Analytical and simulation dimensionless flux results for a 20-segment fracture ... 123 4.28 Simplified multi-stage multi-well simulation model......................................................125 4.29 Simulation results comparison...............................................................................................126 xiii NOMENCLATURE Symbol a Diffusivity foot2/sec Pc Transmissibility Conversion Factor - O Phase Potential psi Porosity - AP Pressure Difference psi At Time Difference day Ax Length Difference foot ^g Gas Viscosity cp ^o Oil Viscosity cp ^w Water Viscosity cp Total Viscosity cp Q Domain - A Flow Area foot2 Bg Gas Formation Volume Factor RB/SCF BHP Bottom Hole Pressure psi Bo Oil Formation Volume Factor RB/STB Bw Water Formation Volume Factor RB/STB Cf Formation Compressibility 1/psi cg Gas Compressibility Co Oil Compressibility Cr RoCk Compressibility Ct Total Compressibility Cw Water Compressibility FBHP Flowing Bottom Hole Pressure Cfdf Dimensionless FraCture ConduCtivity GOR Gas Oil Ratio GOIP Gas Originally In-PlaCe Gp Cumulative Gas Produced h Formation ThiCkness k Absolute Permeability k Multiplying FaCtor kf FraCture Permeability km Reservoir Absolute Permeability Krg Gas Relative Permeability Kro Oil Relative Permeability krog Oil Relative Permeability in Gas-Oil System krow Oil Relative Permeability in Water-Oil System krw Water Relative Permeability kx Reservoir Absolute Permeability in the X-direCtion kxy Matrix Horizontal Permeability ky Matrix Permeability in the Y-direction kz Matrix Permeability in the Z-direCtion xv 1/psi 1/psi 1/psi 1/psi psi SCF/STB MMSCF SCF foot mD mD mD 1/psi mD mD mD mD L Reservoir Extent ft lx, ly, lz Spacings o f Fractures Planes - N Number o f Fracture Normal Sets - Np Cumulative Oil Produced STB nf Number o f fractures - nw Number o f wells - OOIP Oil Originally In-Place MSTB P Pressure psi p Phase - Pb Bubble Point Pressure psi Pbp Bubble Point Pressure psi Pf Pressure at Fracture psi Pi Initial Reservoir Pressure psi P Average Reservoir Pressure psi P wf Flowing Bottom Hole pressure psi qg Gas Rate SCF/day qo Oil Rate STB/day qw Water Rate STB/day re Distance o f Reservoir External Boundary foot Rp Produced Gas Oil Ratio SCF/STB cpRp Cumulative Gas Oil ratio SCF/STB Rs Solution Gas Oil Ratio SCF/STB Rsb Gas Oil Ratio at Bubble Point Pressure SCF/STB Rsi Initial Solution Gas Ratio SCF/STB rw Wellbore Radius foot S Symmetric Factor - S* Effective Saturation - Sg Gas Saturation - Sgc Critical Gas Saturation - groS Residual Oil Saturation in Gas-Oil System - Sorw Residual Oil Saturation in Water-Oil System - Sw Water Saturation - Swc Irreducible Water Saturation - Swc Connate Water Saturation - t Time day T Reservoir Temperature °F T Transmissibility STB/day-psi V Volume foot3 w Width foot We Cumulative Water Injected STB Wp Cumulative Water Produced STB X Reservoir Dimension in X-direction foot Xf Fracture Half-length foot Y Reservoir Dimension in Y-direction foot Ymodel,i Modeled Value Unit o f Output Yobs,i Observed Data Unit o f Output Yobs,max The Maximum Value o f Observed data Unit o f Output Yobs,min The Minimum Value o f Observed data xvii Unit o f Output Subscripts b Bubble point bp Bubble point CD Conductivity cv Control Volume e Number o f elements f Formation, fracture frac Fracture g Gas i Initial j Number count l Fluid phase m Matrix max Maximum min Minimum o Oil obs Observed r Rock R Reservoir t Total w Water, well x X direction y Y direction xviii z Z direction Abbreviations BHP Bottom Hole Pressure FBHP Flowing Bottom Hole Pressure GOR Gas Oil Ratio GOIP Gas Originally In-Place OOIP Oil Originally In-Place LGR Local Grid Refinement MB Material Balance MBPP Material Balance Pressure Profiles NRMSE Normalized Root Mean Square Error NTG Net to Gross Ratio PVT Pressure Volume Temperature IMPES Implicit Pressure Explicit Saturation FEM Finite Element Method DFM Discrete Fracture Model RMSE Root Mean Square Error xix ACKNOWLEDGEMENTS I would like to thank my advisor, Dr. Milind Deo, for his inCredible support, guidanCe, and encouragement during my graduate program at the Chemical Engineering Department. I am very grateful to have had the chance to grow intellectually and personally under his wing. I am also thankful for the advice and support from the members o f my supervisory committee: Drs. John McLennan, Anil Virkar, Philip Smith, and Rasoul Sorkhabi. I thank all the faculty and staff at the Chemical Engineering Department for sharing their knowledge during classes and being resourceful whenever I needed help. I am grateful to the Energy and Geosceince Institute for honoring me with a fellowship and many other awards and programs that supported me financially. I thank all my colleagues and classmates including Palash, Mustafa, Manas, Jing, Hyukmin, and Justin for making this a welcoming and great workplace. I am fortunate to have great friends outside o f my work circle that made Salt Lake City truly a special place, in particular, Yulia for her support during the writing o f this work. Finally, I thank my parents, Vania and Raul for their encouragement, unconditional support, and motivation, and my siblings Ximena and Ricardo for always being there for me. CHAPTER 1 INTRODUCTION 1.1 Unconventional Hydrocarbon Production North American unconventional oil and gas resources have proven to be an important energy asset. So much so that the United States is projected to become the single largest oil producer in the world due to unprecedented development from unconventional reservoirs.1 The U.S. Energy Information Administration (EIA) has noted that the definition o f "unconventional" oil and gas resources is nebulous and depends on varying technologies and economies. In other words, conventional resources are the hydrocarbon fluids that are easy and cheap to extract and exploit, whereas unconventional resources consist o f a large and wider variety o f sources that includes oil sands, extra heavy oil, gas to liquids, and other liquids that require extra technology to produce. By this definition, as oil prices rise and technological advances are made, unconventional resources can migrate into the conventional category. In the following sections, the term "unconventional" will refer to low permeability formations in the range o f micro- to nanodarcies, also known as tight formations. Some o f the most important unconventional reservoirs are the Bakken and the Eagle Ford shown in Figure 1.1. Since hydrocarbon flow can be very low in tight and ultra-tight formations, advanced 2 Figure 1.1 Map o f U.S. most prominent unconventional plays (Source: EIA) technologies are essential for economic exploitation, such as horizontal drilling and hydraulic fracturing. Even though these technologies add an extra price tag to oil companies, increasing oil prices made it profitable for plays like the Monterey/Santos in southern California, the Bakken and the Eagle Ford to be recently developed as reported by the EIA in July 2011.2 Oil and gas production from unconventional reservoirs in the United States is indisputably the main driver for American hydrocarbon production surge. Oil production from tight formations more than tripled in three years, increasing from about 250 Million Stock Tank Barrels (MSTB) per day at the beginning o f 2009 to nearly 900 MSTB per day by 2011. About 84 percent o f all tight oil production in November 2011 came from the Bakken formation in North Dakota and Montana, and the Eagle Ford shale in South Texas. 3 All seven regions: Eagle Ford, Bakken, Niobrara, Permian, Haynesville, Utica, and the Marcellus accounted for 95% o f domestic oil production growth and all o f domestic natural gas production growth from 2011 to 2013.3 Figure 1.2 shows the increase o f oil and gas production from the most important shale plays in the United States. Despite current crude oil price decline resulting in changes to oil production, rates are is still expected to rise for the next 30 years.3 Hence, unconventional reservoirs are proving to be a very important energy asset for the future. 1.2 Research Motivation Even though hydrocarbon production from shale plays has been critical to overall North American hydrocarbon production, ultimate recoveries remain very low. In addition, these resources cannot be properly assessed due to limited information about the reservoir and specifically the lack o f knowledge o f key parameters such as permeability. Sometimes, this Figure 1.2 U.S. tight oil (left) and dry shale gas (right) production (Source: EIA) information can be very expensive or impractical to obtain in these tight formations. This is why developing cheap and fast analytical tools to quickly assess a reservoir's hydrocarbon potential is important. The widespread use o f computer simulation to characterize and predict hydrocarbon recovery has become standard in the industry. However, computer simulation is still handicapped by limitations o f computer power that may result in weeks and even months o f simulation run time. This is why model simplifications o f highly complex systems are sometimes mandatory at the cost o f accuracy. Unfortunately, a well-structured and thorough simulation simplification workflow does not exist. Another important factor in reservoir simulation is the representation o f hydraulic and natural fractures. Historically, there have been many attempts to represent fluid flow through fractures accurately without sacrificing computational time. The Discrete Fracture Network (DFN) technique is one o f the most widely used methods used for representing fractures in simulation. However, like all other methods, DFN has some disadvantages that need to be accounted for. This research work develops semi-analytical techniques to quickly determine reservoirs' potential without the necessity o f expensive well pressure testing or well logging techniques. Another aim o f this research work is to aid reservoir engineers history-match and predict hydrocarbon production by using simplified yet accurate models. Lastly, a modification o f the DFN is suggested to better represent fluid flow through hydraulic and natural fractures 4 1.3 Objectives In view o f the aforementioned research topics and motivation, the following research objectives are presented: • Develop an analytical technique to help quickly and easily identify essential reservoir and fluid properties. o This technique should be able to handle distinct flow patterns. o Different fluids such as oil, water, and gas should be incorporated into this method. o Provide potential reservoir characterization analytics. • Present a comprehensive reservoir simulation workflow that simplifies a simulation model without sacrificing accuracy. o Show that some o f the current widely used simplified models do not accurately represent reality. o Are there any physical phenomena that could be ignored during reservoir simulation? o Develop rules and methods to determine whether reservoir model simplification is achievable. • Modify current Discrete Fracture Network (DFN) implementation in Control Volume Finite Element Method (CVFM) simulators in order to better represent fluid flow in fractures. o Add fracture control volumes to a current DFN implementation. o Keep or improve current DFN meshing techniques to avoid troublesome explicit fracture meshing procedures. 5 6 The stated research objectives are addressed in the following chapters. Background information is first presented and proposed methods and procedures are explained in detail. Application o f such methods are shown in the form o f examples. These examples use reservoir and production information from simulation, field case studies or referenced work. Finally, validation and/or derivation o f the proposed methods are presented in corresponding chapters. Discussion and conclusions highlight the suggested methods' strengths and describe weaknesses. Additionally, important observations and possible improvements are made in respect to the objectives o f this research work. CHAPTER 2 MATERIAL BALANCE APPLIED TO TIGHT FORMATIONS The black oil material balance is a simple method used to estimate original hydrocarbon in place in conventional reservoirs. It has served as a basic fast tool used to assess the production potential o f a reservoir and as a way to either support or question numerical simulation results. Because typical application o f material balance principle is in conventional boundary-dominated reservoirs, low and ultra-low permeability reservoirs (which have seen significant activity in the last few years) have not benefited from the cheap analytical power offered by this method. Hence, a new material balance formulation is required for analyzing these systems. In this study, the conventional material balance is used in a semi-analytical method to construct transient-flow pressure profiles and estimate important reservoir parameters such as permeability. Since the conventional material balance method only requires pressure-volume- temperature data (PVT) and production information, it serves as an essential tool for systems where information is limited. This method is used with a high degree o f accuracy for single-phase flow systems and gives reasonable results for two and three-phase systems. Validation is made through comparison with output from numerical simulation. Examples o f oil, dry gas, and a field data case are discussed after establishing the method steps. 2.1 Background Production o f fluids from tight formations (shales) has completely changed the energy equation for the United States. Liquid production from the Bakken (North Dakota), Eagle Ford (Texas), and plays in the Permian Basin (Texas) accounted for over a third o f the total U.S. oil output in 2014. The development has been very rapid, and the recoveries - particularly o f liquids is low. It is recognized that the matrix permeability is extremely important in initial rates and in ultimate oil recovery. Obtaining good diagnostic information about how the recovery process is unfolding is critical to the success o f the production method. In this section, a semi-analytical conventional material balance method is developed for application in tight formations often referred to as shales. Material balance is one o f the oldest methods developed to assess the hydrocarbon potential o f a reservoir. This formulation was developed by Schilthuis (1936)4 and has since been one o f the main analysis tools for engineers, especially with the introduction o f the straight line equation.5 Although this is arguably the era o f the numerical simulation, material balance remains an accurate and effective way o f determining a reservoir's potential early into production.6 Important work to extend the conventional material balance formulation and apply it to unconventional reservoirs has been done in recent years. In 1994, Walsh et al.7 introduced a general material balance equation that accounts for condensates and it is only recently that commercial numerical simulators are incorporating this concept. A material balance model was presented by Penuela et al. (2001)8 to address naturally fractured reservoirs using a dual-system approach for initially under-saturated reservoirs. Under the same concept, Niz et al. (2004)9 extended this approach for reservoirs with an initial gas 8 cap. Finally, Sandoval et al. (2009)10 included condensates into their framework by incorporating a volatilized oil component in gas. There are other important advancements in the material balance theory that include the dynamic material balance and ‘flowing' oil and gas material balance among others.11-15 Even though significant work has been put into expanding the material balance formulation to unconventional reservoirs, conventional material balance cannot be used for transient systems. This is because material balance is a zero-dimensional approach which assumes that the reservoir behaves like a single tank where all properties are exactly the same throughout the system. This is generally true for homogeneous formations with high permeabilities (typically conventional reservoirs). However, pressure propagation in tight reservoirs is very slow; taking months and even years before the pressure front propagating from the well reaches reservoir boundaries. This means that tight reservoirs remain in the transient flow stage for most o f their productive liv es15 making the material balance technique inapplicable. A novel technique is proposed to help determine a reservoir's potential by using conventional material balance and surrogate pressure profiles. This semi-analytical method fundamentally relies on the material balance formulation (assuming the reservoir behaves like a tank) and combines it with surrogate pressure profiles to obtain pressure behavior at any time during a w ell's transient-flow life. This is accomplished in three main steps: 1. Construct an average pressure profile for the system by using conventional material balance. 2. Determine a pressure profile based on the average pressure curve. 3. Compare pressure profiles to equations from the literature and back-calculate 9 10 important reservoir properties. Pressure profiles for radial single-phase transient flow systems are represented by the exponential integral16 while linear single-phase flow systems are described analytically by the error function17. These expressions are solutions to differential equations coupled with Darcy's law and are broadly used in well testing techniques18. Two and three-phase transient flow systems are typically solved by numerical simulation. In this work, equations such as the error function are used to calculate reservoir parameters when compared to material balance pressure profiles. It is important to highlight that the application o f this method is not limited by the type o f fluid flow (radial or linear) as long as an appropriate pressure equation is used. Also, dry gas material balance can be applied after a modification is introduced into the method to account for compressibility. Since this method is based on the material balance formulation, it shares the same assumptions. There are also other assumptions that depend on the reservoir geometry and rock/fracture properties as described in the upcoming sections. Numerical simulations as well as a real field case were used to test and validate this semi-analytical method. 2.2 General Procedure Before the method is introduced, an important variable which makes this concept different from the conventional material balance is addressed. All material balance expressions contain a parameter known as "Oil Originally In Place" (OOIP) or "Gas Originally In Place" (GOIP). This parameter is typically the unknown during the application o f conventional material balance. OOIP may also be calculated volumetrically if the reservoir's geometry is known. In the case o f transient-flow, however, the extent o f the reservoir cannot be determined due to the infinite-reservoir behavior which is characteristic o f this type o f flow. Hence, rather than calculating OOIP volumetrically, the extent o f the reservoir is simply treated as a variable. This variable is referred to as "Distance" in this work and is defined as a vector normal to the well or hydraulic fracture as shown in Figure 2.1. The detailed method theory is provided in Appendix A. The general application o f the method steps are as follows: 1. Set up the material balance equation such that "distance" is the dependent variable. In other words: d = f (A r, P, Sw, Cw, Cf, Np, Gp, Wp, . . . ) 2.1 All the parameters are described in the Nomenclature Section. 2. Plot P versus distance using equation 2.1 and a production data point at time "t". 3. Extrapolate the average pressure curve (P vs distance) to the known bottom hole pressure at P (0)=Pwf The use o f a surrogate curve that matches the average pressure curve and interpolates to Pw f is highly recommended and used in future examples. 4. Determine a pressure profile expression from the average pressure curve by using a volumetric average pressure relation. 5. Back-calculate permeability using surrogate curve equation. 6. Determine another pressure profile and permeability by following the steps above using a different production data point. This effectively means that the time variable will be equal to effective production time. This method can be used for any kind o f fluid flow and fluid type. The level o f accuracy 11 12 Figure 2.1 Original Oil In Place calculations for two systems showcasing the "Distance" variable o f pressure profiles calculated using this method will depend directly on the material balance equation being used and other factors which will be discussed in the upcoming sections. 2.2.1 Corrected Average Pressure In some specific cases, the material balance may not be as accurate as desired and a modification to the method must be done. There are two main reasons why accuracy in this method may be affected: • The near-fracture effect. Since interpolation to the known bottom-hole pressure takes place near the fracture or well, any phenomenon in this region is ignored (mostly gas coming out o f solution). It was found that this has a 13 negligible effect depending on the Gas Oil Ratio (GOR). • Material balance average pressure. Perhaps the most impactful on this method's accuracy is the difference between the calculated material balance average pressure (equilibrium average pressure) and average pressure during production. Figure 2.2 shows average pressure decline in a reservoir where a well has been shut after one year o f production. Average reservoir pressure continues to decline even after shutdown until it reaches an equilibrium point as determined by numerical simulation. Material Balance average pressure is the equilibrium average pressure at t = rc>, rather than at t = 1 year. The difference, AP, is a direct result o f fluid compressibility, partial bubble point transition, and so forth. Average pressure difference is only substantial for cases o f high GOR or dry gas reservoirs and may affect pressure profiles calculated by the material balance method. In cases where the material balance may not be sufficiently accurate, the corrected average pressure method presented earlier is applied to obtain better results as shown in the upcoming examples. To correct for the average pressure difference seen in Figure 2.2, the following relation is applied to gas reservoirs in conjunction with the method: The development o f this equation can be found in the upcoming sections in more rigorous forms. Equation 2.3 shows the oil equivalent o f equation 2.2: 2.2 2.3 14 5350 5300 - 5250 - a.o t 5200 - <2> 5150 - 5100 - 5050 Shut-in Average ~ \ j Pressure \ Equilibrium \ Average Pressure \ / - >-AP - 50 psi 1 1 1 1 1 1 1 1 1 4 5 6 Time, years 10 Figure 2.2 Difference between equilibrium average pressure (material balance pressure) and "real time" average pressure Equation 2.3 can be rearranged in the following manner as shown in the corrected average pressure theory section: wh [ L $ ( P ( x ) ) S 0 (P( x ) ) 2.4 5 .6 1 4 5 8 J0 5 0 (P (x ) ) dx = OOIP- Np The application o f these equations to the material balance method will increase accuracy. However, some knowledge o f relative permeability curves is needed when dealing with multiphase systems. Example 1 shows the application o f the material balance method to a multiphase case where no modification is needed. Example 2 shows a typical case where material balance alone may not meet accuracy requirements and the corrected average pressure modification is introduced. Example 3 shows the application o f the material balance with average pressure correction for a single-phase compressible fluid. 2.3 Case Examples 2.3.1 Example 1: Multiphase Linear Flow into a Vertical Fracture Consider a single hydraulic vertical fracture draining from a low permeability matrix media with homogenous fluid saturations as depicted in Figure 2.3. This represents a portion o f a multistage hydraulically fractured reservoir. The following assumptions are made with respect to the physical system. • Infinite conductivity vertical fractures. • Fractures can be adequately represented by a vertical plane. • All fractures are identical and have similar production. 15 Figure 2.3 Linear flow into a vertical fracture 16 Using the production data and fluid information provided in Table 2.1, the method steps in the previous section are followed: 1. Setting up the conventional black oil material balance equation in the form o f equation 2.1: d = 5.61458- B0i(NP(B0 + (RP - Rs)Bg) - (We - Wp)Bw) 2.5 2$wh ((1 - SWc) ((B0 - Boi) + (Rsi - Rs)Bg) + {CwSwc + Cf )APB0l^ Note that since there is no gas cap in the system then m = 0. 2. Using production data at an arbitrary time o f 90 days, an average pressure plot was built as shown in bold in Figure 2.4. 3. A surrogate curve based on the single phase linear flow solution was used to match the average pressure curve and interpolate it to the known bottom-hole pressure o f 500 psi. This curve is shown in Figure 2.4 and was obtained by iterating through a in equation 2.6. ' ' ‘ ' ' 2.6 P(x) = Pwf + {Pi - Pwf) (6.285311x\ 0.0897632 _ / 39.5051X2 erfl-----------)+----------------Vat\ e at - 1 \ Vat ) x \ 4. Once a satisfying value for a was found, a pressure profile expression was determined by substitution and differentiation o f equation 2.6 as shown by equation 2.7. Table 2.1 Reservoir and operational parameters for simulations Oil Production STB Gas Production MSCF GOR at Pb scf/stb BHP psi 90 days 1 year 90 days 1 year Example 1 2147 4218 3050 5850 700 500 Example 2 2286 4540 8120 15300 1200 1000 Example 3 - - 71782 231000 - 500 17 Figure 2.4 Average pressure profile matched by a surrogate curve P = f i P d V / > ( * ) dx 2.7 J > L The resulting expression is Miller's single phase linear flow solution as shown in equation 2.8. P(x) = PWf + (Pi - PWf ) erf ( '6 .2 8 5 3 1 1 x \ V a t ) 2.8 A set o f pressure profiles at arbitrary production times o f 90 days and 1 year were produced. Validation o f these results was done through numerical simulation as shown in Figure 2.5. 5. Equation 2.8 represents Miller's single phase linear flow pressure behavior, where a is diffusivity as shown in equation 2.9. 18 Figure 2.5 Calculated material balance pressure profiles versus simulation _ k 2.9 a where Ct = SoCo + SgCg + SwCw + Cf, and can be estimated using weighted average values. By simple substitution, permeability was calculated to be around 180 nD, a very good estimation when compared to the actual absolute permeability o f 200 nD. The agreement between the calculated MBPP's (Material Balance Pressure Profiles) and numerical simulation pressure profiles is remarkable as shown in Figure 2.5. The slight variation near the well is due to gas coming out of solution. In certain cases where GOR is not very high, gas contribution to this process can be safely ignored by simplifying equation 2.5 into equation 2.10. , c „ „ co B0i(NPB0 - (We - Wp)Bw) 2.10 d = 5 .6 1 4 5 8 -------------------------------------------------- -------------------------- r- 2$w h ( (1 - SWc) (B0 - Boi) + ( CwSwc + Cf)APBoij Equation 2.10 comes in handy for quick calculations where GOR is not considerably high. As GOR increases, the accuracy o f MBPP's and permeability calculations may be reduced. However, even in cases where GOR is quite high, MBPP's are still a decent approximation to a full numerical simulation model. This is shown in Example 2. It should be noted that a ‘manual' interpolation from the average pressure curve to the known bottom-hole pressure is also acceptable. The advantage o f having surrogate curves is their easy incorporation into the methodology. Also, surrogate curves are in many cases solutions to known fluid flow phenomena that contain important physical properties that can be back-calculated. In other words, parameters such as permeability or diffusivity are readily obtainable with this method. However, one should be cautious as most o f these surrogate curves are based on one-phase flow. Other surrogate curves similar to Miller's expression can be found in the oil and gas literature and transport books such as Bird's Transport Phenomena (1960).19 19 2.3.2 Example 2: Multiphase Linear Flow into a Vertical Fracture and High GOR In this case, the same steps as in Example 1 were followed using production data from Table 2.1 to produce pressure profiles shown in Figure 2.6. As mentioned earlier, the match between calculated pressure profiles and output from numerical simulation is not as good as in Example 1. However, MBPP's are a very good approximations given that no knowledge about matrix permeability, relative permeabilities nor production rates are required for this method. 20 7000 Simulation at 1 year Simulation at 90 days Material Balance 6000 1000 0 0 100 200 300 400 500 Distance, ft Figure 2.6 Calculated material balance pressure profiles versus simulation 2.3.3 Example 2 (Revisited): Multiphase Linear Flow into a Vertical Fracture and High GOR The implementation o f equation 2.4 to the material balance method for increased accuracy is shown in this section. In this procedure, a new term "L", defined as an arbitrarily large number, is introduced. 1. Using surrogate pressure equation 2.8 and an arbitrary value for a, P(x) for x E[0,L] was plotted. 2. The product in the left hand side o f equation 2.4 was calculated numerically using a pressure profile as determined in step 1. 2.8 21 wh If L<rVp( PK( xJ)J) S0OK( P(K x )J)J ^x = 0 Ql p - N 2.4 5 .6 1 4 5 8 iJ0 B0 (P( x ) ) ' p Porosity is a function o f pressure, but it may be assumed constant and be taken out o f the integral for most cases. Also, before gas comes out o f solution oil within the reservoir, saturation profile is best represented by an error function. After gas comes out o f solution, oil saturation profile can be approximated by a linear function. Hence, if instantaneous gas oil ratio remains constant, oil saturation may be estimated by: /<' 6.285311x\ 2.11 S0(x) = S0*(BHP) + (Soi - PS0(BHP)) erf ( 6 .2 8 5 3 1 1 x \ 2.12 / Va t where So*(BHP) is defined as: Sr0 (^Br>PrA) - Sr0 ie r^f (/-6--.-2--8=1--1--x-)\ c * (n u p \ __________________ V (n n r ) = /6 .2 8 5 3 1 1 x \ 1 - e r f (- V i H Oil saturation at the bubble point is: s 0iBP) = 1 - Sw‘( 1 - c; (P* - p ‘) ) 2 1 3 ' 1 + Cr (P„p - Pt) If instantaneous gas oil ratio shows signs o f gas coming out o f solution within the reservoir, a linear relationship describes oil saturation as shown in equation 2.14. (1 - Sw(BHP) - S JBHP ) - S0(B P ) ] ( x bp - x) 2 1 4 S0(x) = S0(BP) + ± ------------------------------------------------------- ---------------- Xfop In this case, gas and water saturations at the bottom-hole pressure can be calculated either analytically or iteratively by solving equations 2.15, 2.16, and 2.17 at the same time. Kro (^rroJ)Swc 22 Krnw \ I Kroa \ s.. .. \1 2.15 ( v ) + ^ rw) i ( K + ^rg\ (Krw + Krg) \ ^V^roJSwc ' '-(Aro) Swc ' -1 ^rw - WOR 2.16 Kro fto ^0 Kra V-qBq 2.17 - ^ -(GOR - Rs) ^ - f ^ro 3. OOIP was volumetrically determined by assuming that "L" from step 1 defines the extent o f the reservoir. This was then used to calculate the right hand side o f equation 2.4. 4. A value o f a was iteratively found in equation 2.8 so that the equality in equation 2.4 was satisfied. 5. Steps 1 through 4 were repeated for arbitrary production times o f 90 days and 1 year. Resulting pressure profiles are shown in Figure 2.7. It is clear from Figure 2.7 that pressure profiles determined using corrected average pressures are more accurate than regular material balance pressure profiles. The trade off, however, is a more complex and involved method that requires some knowledge about relative permeability data. The choice o f method is then left to judgement according to acceptable accuracy and available resources. Equation 2.4 may only be satisfied with the correct pressure profile using oil production. The gas equivalent is shown in equation 2.18. Both equations can be used separately or in combination in a multiphase system. wh fL f<p(P(x))S0(P(x))Rs(P(x)) _ ^(P(x))Sg(P(x))\ n ^ 2.18 ■ +---------- r-■-"T----- ) OX - uUIr - Gy, 5.61458J J0 { B0(P(x)) Bg (P(x)) Jv 23 Figure 2.7 Pressure profiles determined using corrected average pressures 2.3.4 Example 3: Gas Linear Flow into a Fracture This is a good example where application o f the material balance method requires the corrected average pressure for substantially better accuracy. Fortunately, the procedure in this case is much simpler than Example 2 due to the fact that this is a single-phase system. The following steps summarize the application of the material balance method for dry gas systems. 1. The dry gas material balance equation was set in the form o f equation 2.1. d = 5.6 1 4 5 8 - Bgi (GpBg + WpBw) 2.19 2 0w h ( ( 1 ^Wc) (Bg Bgi) + (CWSWC + Cf )APBgi) 2. Using an arbitrary production time o f 90 days from Table 2.1, a plot P o f vs d was built. 3. A surrogate curve based on the gas linear flow solution was used to match the average pressure curve and interpolate it to the known bottom-hole pressure o f 500 psi. This curve was obtained by iterating through a in equation 2.20. 2.20 24 .C I V + ^ - W f f 2 8 5 3 1 1 * '1) dx _ wf + ( r ' rwr ) e , ' ( Va t ) P(x ) = - L 4. Once a satisfying value for a was found, a pressure profile expression was determined by substitution and differentiation o f equation 2.20 as shown by equation 2.7 in Example 1. I 7 /6 .2 8 5 3 1 1 x \ 2.21 P ( x ) = )Pw/2 + (P.2 - P w / 2) e r f ( v _ Equation 2.21 describes an estimate to compressible linear flow solution and was plotted for arbitrary production times o f 90 days and 1 year for comparison with numerical simulation in Figure 2.8. It is clear from Figure 2.8 that MBPP's do not match numerical simulation and the usage o f corrected average pressures may be required. For this purpose, steps similar to Example 2 were followed: 1. Using surrogate pressure equation 2.21 and an arbitrary value for a, P(x) for x E [0,L] was plotted. Note that "L" is defined as an arbitrarily large number. I " ~ ~ /6 .2 8 5 3 1 1 x \ 2 2 1 P (x ) = j p w /2 + ( P,2 - P„/ 2 ) e r f ( ^ ^ - ) 2. The product in the left hand side o f equation 2.22 was calculated numerically using a pressure profile as determined in step 1. Average Pressure, psi 25 Distance, ft Figure 2.8 Pressure profiles determined using the material balance method wht i C 1 2.22 1 dx = GOIP - G r f 5 .6 1 4 5 8i JJn0 BBJgP(~P(( x ) ) ~ " p Note that equations 2.4 and 2.22 are similar in nature, except that equation 2.22 is easier to calculate as it applies to a single-phase system. 3. GOIP was volumetrically determined by assuming that "L" from step 1 defines the extent of the reservoir. This was then used to calculate the right hand side of equation 2.22. 4. A value o f a was iteratively found in equation 2.21 so that the equality in equation 2.22 was satisfied. 5. Steps 1 through 4 were repeated for arbitrary production times o f 90 days and 1 year. Resulting pressure profiles are shown in Figure 2.9. 26 6000 5000 '\A Q^) 4000 i-3iyf>l © t 3000 6«0 >V < 2000 1000 0 Figure 2.9 Pressure profiles determined using corrected average pressures Inspection of Figure 2.8 and Figure 2.9 clearly show substantial improvement by using corrected average pressures in the material balance method. For single-phase cases, this implementation is simple and is highly recommended for compressible fluids. 2.3.5 Example 4: Permian Basin Field Case This field case study consists of a multistage hydraulically fractured horizontal well with over 90 hydraulic fractures. The reservoir is composed o f several layers with different petro-physical properties. However, the ranges for these properties are narrow and a weighted average was calculated for the whole system. Artificial lift is introduced to help increase drawdown after 70 days o f production. The Distance, ft bottom-hole pressure takes about 100 days to stabilize. After this time, the bottom-hole pressure reaches a more stabilized value o f about 500 psi as determined by empirical multi­phase flow correlations. Even though MBPP's may be determined with varying bottom-hole pressures, a constant value o f 500 psi was used to keep the method simple. Using a procedure identical to Example 1, permeability was calculated as a function of time. Figure 2.10 shows normalized fluid production rates and calculated permeability for a production time o f almost 2 years. "Material balance permeability" trends see an initial increase until it levels o ff at around 600 nanodarcies. This is expected given that bottom-hole pressure stabilizes during this time. The estimated permeability o f 600 nD is consistent with the range o f permeabilities measured for samples from this well. Although permeability calculations by using material balance are possible, the analytical tools made available by this method may be of more interest. One important aspect of this plot happens at 600 days, where the permeability skyrockets. Looking at oil and gas rates there is no reason for the permeability to escalate, but a big surge in water production occurs at this time. Because the conventional material balance used in this particular case accounts for multiphase behavior, the spike in water production translates in higher observed permeability. Average reservoir pressure decline was also calculated by assuming a stimulated reservoir volume based on fracture dimensions as shown in Figure 2.11. This field example demonstrates the applicability of the method developed in this research work. It should be noted that the production rates did not stabilize for some time and that the well was shut down from time to time. The bottom-hole pressure also varied over the production time studied. Despite these complications in data collection, material 27 28 Figure 2.10 Transient material balance method application to a field case balance permeability was calculated and analyzed. Not only does the method provide an estimate of permeability that is consistent with production and pressures, but it also helps analyze major disturbances (like the water surge as seen in Figure 2.11. Material balance permeability also seems to be sensitive to long shutdowns as seen at around 300 and 500 days of production. Other helpful production analytics can be implemented into this workflow thus making this a potential topic for future research work. Material Balance Permeability, nD 29 Figure 2.11 Field case calculated average pressure decline 2.4 Other Applications As discussed earlier, surrogate curves used in this method are derived from one-phase equations developed by Miller (1962)17 and presented by Katz et al. (1959).16 Therefore, there are some physically important variables that can be extracted from these curves as is the case o f permeability and diffusivity. Even though calculated permeabilities using the material balance method were very good estimates, the actual power o f the method lies in its potential use for analytics. This was shown in Example 4, where permeability was calculated using production data points for about 2 years. A material balance permeability versus time plot was used to quantify the impact o f artificial lift on production and assess the impact o f a production water surge. As a rule o f thumb, Figure 2.12 shows the material balance permeability behavior o f a true infinite reservoir. The flat portion represents the estimated permeability and any deviations could be a result o f significant fracture interference, fracture closure or boundary dominated flow. Hence, this plot can be further studied and important reservoir behavior can be derived as a consequence. Another important piece o f information that can be derived from this method is pressure decline. Figure 2.11 shows that reservoir average pressure decline computations are possible with this method. Pressure decline at particular points in the reservoir are also possible as shown in Figure 2.13. Comparison with pressure sensor information would provide decent insight into the system's behavior. 30 Figure 2.12 Calculated material balance permeability for an infinite reservoir 31 Figure 2.13 Pressure decline over time at a point 50 ft away from a draining fracture As described in Example 2 (revisited), fluid saturation calculations are also possible after pressure profiles are established. Even though important assumptions regarding the system were made, this serves as an example o f potential opportunities for analysis. Figure 2.14 shows the oil saturation profile developed after an arbitrary production time. These are a few examples that serve as supplementary analytical tools for reservoir assessment. Analytical options using the transient material balance workflow are significant and easily available with minimum information about the reservoir. 32 l o 0.3 - 0.2 - 0.1 - 0 -I----------------------------1---------------------------------------------1---------------------------------------------1---------------------------------------------1--------------------------------------------- 0 10 20 30 40 50 Distance, ft Figure 2.14 Calculated oil saturation profile 2.5 Key Findings In this work, it was shown that pressure profiles for a transient-flow system can be determined using a method based on conventional material balance. The method's basic steps and theory were presented and its results were compared to numerical simulation and applied to an unconventional field case. Material Balance pressure profiles seem to be a decent approximation to single and multiphase systems. Important reservoir information such as permeability and average reservoir pressure decline can be extracted from the use of this method. When compressible fluids are flowing to the wellbore, a modified step was presented to improve accuracy. The applicability of the method was demonstrated by comparing results o f several cases with results from a reservoir simulator. A field case from the Permian Basin further validated the method. 33 The material balance approach combined with pressure profile calculations uses cumulative fluid production data and fluid properties. In low-permeability shale reservoirs, early transient flow is expected to last for a long period o f time, and this method allows for the extraction o f important reservoir information. While weaknesses on this method have to do with the assumptions o f the material balance formulation, it was shown that this workflow can be used as a reservoir assessment tool requiring only limited data. 2.6.1 Material Balance Theory In this section, a brief summary o f the material balance theory is introduced. Consider a case where a single-phase incompressible fluid flows linearly into a fracture as shown in Figure 2.3. If the flow is transient, a pressure profile is developed during production and may be represented by equation 2.23 as described by Miller. Assume that half o f this system is represented by Figure 2.1 and that the "distance" variable in this system is an arbitrary large number "L". Consider a second system identical to the first one, but its "distance" variable is "L*", where L* > L. Both systems are identical, except that the first system is smaller than the second. If both systems are flowing during the transient stage, then pressure fronts have not yet reached the boundaries o f the system at L * or L . A sample transient pressure profile developed in both systems is shown in Figure 2.15. 2.6 Method Development 34 Figure 2.15 Typical error function pressure profile It is clear from Figure 2.15 that P(L) = P(L*) = Pi, hence the pressure for any distance greater than L will be the initial pressure. It follows from this idea that any reservoir size larger than L will span identical pressure profiles. The same concept can be applied to average pressure profiles. Through application o f volumetric average pressure relation 2.7, equation 2.24 is obtained and an average pressure profile can be plotted as shown in Figure 2.16. p = I o P d V _ t i P (x ) dx 2 7 P(X) - Pwf + (Pj - Pwf) erf ^6.285311x^ V kt + L 0.0897632 X -v<pc(\ 39.5051X2 e kt 2.24 1 35 Pi L L* ------------------------------------ > Ln Distance away from the fracture Figure 2.16 Average pressure profile plot calculated by varying L If an average pressure curve such as the one shown in Figure 2.16 is obtained by an alternative method, a pressure profile may be "back-calculated" without knowledge o f permeability, viscosity, and so forth. Also, permeability and diffusivity values can be back-calculated by simple substitution. The novel idea behind the transient flow material balance is the fact that a conventional material balance formulation is used to independently calculate average reservoir pressures at different values o f "distance" or in this case L, L* and Ln as shown in Figure 2.17. The material balance method provides with a good average pressure profile until it breaks down in the near-well region. There are several reasons for this behavior, one o f which is partly due to the concept that reservoirs of varying extents (quantified by the use 36 Material Balance Pressure Average Pressure / / I / / / L L* ------------------------------------ > Ln Distance away from fracture/well Figure 2.17 Material balance average pressure compared to an average pressure surrogate curve of the "distance" variable) span identical pressure profiles. This only holds as long as the extent of the reservoir is "large" enough for the flow to be considered transient. Since material balance makes use of production data, different pressure profiles can be calculated at different production data points. Diffusivity and permeability calculations can also be made for a range of production data points for analysis. 2.6.2 Corrected Average Pressure Theory It was found that there is a difference between material balance (equilibrium) average pressure and average pressure while producing. The aim here is to find the average pressure difference, so that a corrected average pressure can be determined. To do this, a simple 37 equation relating equilibrium and producing average pressure is needed, and is presented as follows. Consider a simple case of a reservoir consisting of two tanks. Both tanks have 100% porosity, 100% oil saturation, and have the same volume. Tank 1 has a pressure of 2000 psi and tank 2 has a pressure of 5000 psi as shown in Figure 2.18. Since both tanks are identical, except for their pressures, the instantaneous average pressure of this system is 3500 psi, which is the arithmetic average of 2000 and 5000. However, simulation shows that the equilibrium average pressure of this system at equilibrium is 3540, a 40 psi difference from the instantaneous average. To equate both, instantaneous and equilibrium systems, their corresponding oil volumes are brought to atmospheric pressures (produce this fluid). The atmospheric fluid volume of systems must be exactly the same. This is equity is visualized in Figure 2.19. Mathematically, Figure 2.19 is expressed in equation 2.25. Bo i B o2 B0 (P) V-l + V-l = Vi + V2 2.25 Or, Bo 1 B o2 2.26 Tank 1 Tank 2 Porosity = 1 Oil Saturation = 1 Volume = V Pressure = 2000 psi Porosity = 1 Oil Saturation = 1 Volume = V Pressure = 5000 psi Figure 2.18 System consisting of Tanks 1 and 2 38 Figure 2.19 Equivalence between Tanks 1 and 2 at producing state and Tanks 1 and 2 at equilibrium state By taking into account that both tanks have the same volume, and expanding this to an arbitrary number of tanks, the following relation can be established. Z " da,ti LL 227 i = i B^i = B0jJ(FP)) where di and L are individual tank length and sum of all tank lengths correspondingly. Equation 2.27 can also be expressed as an integral. 2.28 I 1 1 L d x = Jo B0 (P) B0 (P) A more rigorous expression that takes into account porosity and oil saturation can be similarly derived. 2.29 Io L$ ( P ) S 0 (P) L $ (P ) S 0 (P) OX = Jo B0 (P) B0 (P) Similarly, this relation can be obtained for the gas case as shown by equation 2.30. J0 39 L0 (P )S g ( P ) ^ = L0(P)Sg (P) 2.30 /o W X W Equation 2.31 results from equation 2.30 for the case of a rectangular reservoir as shown in Figure 2.3. wh f L0 (P )S o (P) ■ox = 0 0 /P - W*, 2.31 5.61458J Jo0 B0(P) --------------- For a rectangular multiphase reservoir, a gas analogous expression is similarly derived. wh (■V0(P)So(P)fis (P) , </>(P)Ss ( P ^ , „ 2.32 +----- „ . - ) ox = GO/P - Gr 5.61458 JIo tV B0(P) ' Bfl(P) 7 .................. P Equations 2.31 and 2.32 may be used separately or in combination to obtain more accurate results than the basic transient material balance approach. The solutions to these equations can be reached numerically or analytically depending on the PVT data format. CHAPTER 3 RESERVOIR SIMULATION SIMPLIFICATION WORKFLOW FOR HYDRAULICALLY FRACTURED RESERVOIRS Production from unconventional formations, such as shales, has significantly increased in recent years by stimulating large portions of a reservoir through the application of horizontal drilling and hydraulic fracturing. As a result, reservoir numerical simulation of tight and ultra-tight reservoirs has become the standard tool to assess and predict production performance from these unconventional resources. Because many of these unconventional fields are immense, consisting of multistage and multiwell projects, simplification of simulation models is common both in the industry and academia. It is important to represent the results from a full-scale model by performing computations on smaller models that capture the physics while keeping the computational time manageable. First, the most widely used simplified models are shown not to be representative of a full-scale reservoir model. These models do not account for many important factors such as fracture interference, well interference, and so forth, resulting in inaccurate rate predictions, especially when interference or boundary dominated flow take place. A novel, rigorous workflow is proposed, and compared to other literature models and to a full-scale model. Procedures, results, observations, and limitations of this workflow are then discussed. The models that result from the application of the proposed simplification workflow can predict the fluid rates, cumulative production, and gas oil ratio with higher accuracy than popular simplified models while retaining low run times. 3.1 Background For over five decades, numerical reservoir simulation has been one of the most important reservoir engineering tools. Simulation continues to be the best way to quantitatively describe multiple phase fluid flow behavior in highly complex, heterogeneous systems.20 Over time, reservoir simulation has become the final hydrocarbon potential assessment tool, if not the only one, due to its incredible flexibility and success during exploitation of underground natural resources. This is why continued research efforts are made to increase the accuracy and flexibility of these simulators and to decrease the computational time required to run single or multiple studies. Generally speaking, the process of building a reservoir model can vary depending on the objectives of a project. Reservoir models may encompass full, comprehensive fields that are built based on geological, seismic, and petrophysical interpretations. These static models are then upscaled to meet computational limits and handed over to engineers for fluid recovery optimization through well management strategies, tertiary oil recovery, and production schedule based on the current economic climate. Simulation models may also be very simple to either study particular aspects of smaller reservoirs or for research purposes where analytical expressions are validated or supported. The latter situation is very common, particularly for projects where a full reservoir simulation study is either too costly or large enough for simulation run times to be impractical. In these cases, simplified models are used both in academia and the industry. 41 Unfortunately, the success accrued by numerical simulation has been somewhat taken for granted, leaving space for models that do not accurately represent physical reality. This is the case for the simplified "single fracture" model which is used widely in the industry and academia.21-27 Recent drilling technologies in low permeability reservoirs include horizontal drilling and multistage hydraulic fracturing. Such systems are typically modeled as shown in Figure 3.1. Such systems are widely known to be impractical due to immense simulation run times and also the effort put forth to build such models is considerable. For this and other reasons, simplification techniques have been proposed to avoid full scale studies; such is the case of the single fracture model. The single fracture model is a result of a concept known as Stimulated Reservoir Volume (SRV). This is the volume stimulated by hydraulic fracture half lengths, fracture height, and number of propped fractures.28 The single fracture approach consists of a simulation model containing a single fracture draining from an allocated volume of reservoir. Typical dimensions of the model correspond to the length of the fracture and interfracture spacing as shown in Figure 3.1. The advantage of this approach is the substantially low computational requirements by assuming all fractures in the field behave similarly, hence they can be represented by a single fracture. In the single fracture model, total fluid production of the entire horizontal well is trivially calculated by multiplying production from the single fracture model by the number of fractures in the well as described by equation 3.1. X - W-f %Single fracture 31 where, X is well fluid (oil, gas or water) flowrate or cumulative production; nf is the number 42 43 Figure 3.1 Aerial view of typical multiwell and multistage reservoir simulation model and how it may be simplified of fractures, and xsingle fracture is the single fracture fluid flowrate or cumulative production. Similarly, total production for a multiwell project is calculated from the single fracture model as shown in equation 3.2. X - ^w^-f^Single well model 3 2 where nw corresponds to the number of wells. The problem with the single fracture model is that while it properly represents early transient flow, also called infinite-acting reservoir behavior29, it fails as soon as boundary dominated flow takes over. Therefore, one of the objectives of this work is to establish the validity of the single fracture model and to propose a new, more accurate simulation modeling technique. This new technique retains the low simulation run times which is key to the single fracture model, while accurately accounting for fracture and well interference effects. 3.2 Workflow Components In this workflow, an alternative to the popular single fracture model is presented. Simulation models resulting from the proper application of this new workflow should retain the single fracture's short run times while achieving a higher degree of accuracy. The primary objective of this workflow is to establish a standard in the literature in terms of full-scale simulation model simplification. This simplification workflow standard is meant to guide and provide the simulation engineer with educated criteria in the sometimes blind mission to simplify full-scale models. The structure of the simplification workflow revolves around the consideration of fluid flow phenomena taking place in numerical simulation. The workflow is formally introduced after important fluid flow phenomena is addressed. 3.2.1 Interference Effects As presented earlier, one of the fluid flow phenomena ignored by the single fracture model is the interference effect that result when two transient pressure fronts meet either from neighboring wells, fractures or both. As shown in Figure 3.1, the single fracture model is limited to dimensions corresponding to the horizontal fracture half-length and interfracture spacing. Assuming two neighboring fractures share the same operating conditions, such as flowing bottom-hole pressures and fluid rates, these two fractures are expected to reach mutual interference at some point between them. At this point, the single fracture model imposes a no-flow boundary condition, which is expected to mimic interference effects. Unfortunately, in the upcoming sections, it is shown that this representation may not be accurate enough when a single fracture model is compared to a 44 full-scale model. Interference may substantially affect production from fractures and even wells depending on the operating conditions, fluid type, spacing, and diffusivity. Elliot (1951)30 has shown the great effects of well interference on production to the point where wider well spacing was highly recommended in conventional reservoirs. Given the evidence supporting interference effects on fluid production, it is hard to neglect inter-fracture and interwell interference on simulation models, especially if operating conditions differ from fracture to fracture or well to well. 3.2.2 Boundary effects Conditions imposed to simulation boundaries are several and serve a number of purposes. Such is the case of the constant pressure boundary condition which is applied to represent aquifers or the no-flow boundary representing nonpermeable faults or discontinuities in a reservoir. Model boundaries are always present in some form or another in any simulation model; the most common version of boundary conditions is the no-flow boundary which works by simply setting outside fluxes to zero, therefore nullifying any fluid contribution from the boundary. Under this condition, boundary dominated flow begins as a new flow regime which can be identified through production analytics. Reservoir boundaries and their conditions not only greatly affect hydrocarbon flow, but also delineate a reservoir's original hydrocarbon in place. This is yet another disadvantage of the single fracture model that accounts for the volume stimulated only by internal fractures. External fractures, which are the first and last fractures in a horizontal well, have greater hydrocarbon potential depending on reservoir boundaries as shown in Figure 3.2. 45 46 Figure 3.2 Types of internal and external fractures and their corresponding volume assignments (top view) 3.2.3 Multiple Hydraulic Fracture Representation Hydraulic fracture creation and propagation are challenging topics of research in geomechanics which force stimulation and simulation engineers to work with what is available. Microseismic fracture mapping is the closest the oil and gas industry has come to imaging hydraulic fracture shape; however, this method does not account for fracture closure after treatment and may overestimate the contact area between the reservoir and the well. For most practical applications, geophysicists and engineers agree with an adequate fracture half-length and height for a rectangular shaped fracture to be used in numerical simulation. Under the assumptions made about hydraulic fractures' dimensions, fracture generalization in horizontal wells is widely practiced in the literature. If the injection operating conditions are identical for all fractures in a given horizontal well, it is believed that all fractures will possess similar dimensions and rock-flow properties. This has generated the assumptions that all fractures are identical and will behave similarly as stated in most analytical and numerical research works by both the industry and academia. Multiple hydraulic fracture representation through a single fracture is one of the most important assumptions in computer modeling of simplified hydraulically fractured reservoirs. In this workflow, this concept is also used, but instead of a single fracture generalization, two or even more generalizations can be made. In summary, these generalizations come in the form of internal fractures and external fractures. Types of internal and external fractures can be easily identified as shown in Figure 3.2. In this work, the fluid flow behavior and contribution to production from internal and external fractures are taken into account to achieve more accurate results. 3.2.4 Other Simulation Phenomena There are other important components that should be considered when making decisions about a computer simulation model. Most simulation phenomena such as reservoir boundary conditions and interference effects that influence fluid flow may be common while others may not show up as often, such as vertical gridding (gravity segregation), grid geometry, natural and hydraulic fractures, property heterogeneity, and special well completions, among others. Because the objective of this paper is to establish a generalized workflow rather than presenting every possible fluid flow mechanics applicable to particular models, only a word of caution is presented here. Engineering judgement, dimensionless analysis, and even a trial and error analysis may suggest whether 47 48 a particular simulation phenomenon should be neglected or not. In the following sections, application of the aforementioned simulation phenomena is discussed. The core of this simplification workflow revolves around fluid flow phenomena and the use of symmetry to come up with a simplified simulation model truly representative of a full-scale model. Other simulation phenomena not specifically discussed here should be included in the particular workflow used to build these particular models. 3.2.5 Symmetry Application As mentioned earlier, the single fracture model relies on the assumptions that all fractures propped in a horizontal well and under the same operating conditions behave identically. Therefore, as shown in equations 3.1 and 3.2, the behavior of a single well is represented by multiplying the results from a single fracture model by the number of fractures in the well. This application of symmetry, however, does not capture important fluid flow phenomena that was discussed in the previous sections. Models can significantly benefit from any symmetrical behavior as long as all contributing fluid flow phenomena is captured. A simplified simulation model is a small, representative piece of a full-scale simulation model that captures all significant physical and chemical phenomena, thus decreasing required run times while remaining reasonably accurate. Under this definition, the modeler must focus on all contributing forces to fluid flow if a simplified system is to be built. In the following section, the proposed simplification workflow is introduced where important fluid flow components are revisited and combined with a symmetry application that results in a reliable simplified model. 49 3.3 Proposed Simplification Workflow The solution to most numerical simulation models is dependent on the contribution of various physical forces and fluid flow phenomena. Some of these factors may be deemed negligible for practical reasons and excluded from the simulation altogether. A simple example of how these factors combine to reach a solution for a specific problem is shown in Figure 3.3. After all significant factors have been identified for a particular project as depicted in Figure 3.3, the simplification process may begin. As mentioned in the previous section, one of the many assumptions made in the literature is the multiple hydraulic fracture representation which states that all hydraulic fractures propped under similar operating conditions behave similarly. Hence, a single (or multiple) representative fracture is enough for numerical and analytical solutions which means that the "Fracture variation" component of the solution process in Figure 3.3 may be crossed out. Similarly, other contributing factors may be crossed out based upon assumptions or previous analysis. Once a solution process is built and simplified, the concept of symmetry is applied by capturing all remaining factors. This is explained more thoroughly in the upcoming examples. Lastly, results from the simplified model must be modified to represent full-scale simulation. Equations that make these conversions are intuitive as in the case of the single fracture model. However, as models grow in complexity, these equations may become less obvious. Equation 3.3 is designed to help keep track of important elements in simple and complex simulation models. 3.3 50 Solution \__________ V 00 Stage/Fracture Well * Side Interference ■ Interference T boundary L J -------- --------- J * . Top/Bottom + Well + Heterogeneity boundary trajectory V J v J * Completion (layers) V ' ^ * Fracture variation L wsJ + Internal Fractures v J * External Internal External Fractures Wells iP Wells L-- i------- j L J \ ---------- ---------J Figure 3.3 Example of a solution process diagram that shows various force contributions to flow in a simple hydraulically fractured reservoir. where e is the number of elements; Si is the symmetric factor for the ith element; k is the multiplier factor for full model representation and Xi is the cumulative production or production rate from the ith element. The application of this equation is shown in the next section. The simplification workflow steps can be summarized as follows: 1. Build a solution process diagram by accounting for all present phenomena in the system. 2. Simplify the solution process diagram by crossing out negligible factors based on assumptions or previous analysis. 3. Prepare a unit model representative of a full-scale model by applying symmetry and accounting for all remaining factors in step 2. 4. Transform simplified results to full-scale results by applying equation 3.3. 51 3.4 Case Examples 3.4.1 Single Well Case In this instance, the single fracture model and a model prepared after application of the simplification workflow are compared to a full-scale, black-oil, single horizontal well model. To give some context of the problem at hand, the full-scale model is shown in Figure 3.4 with basic model properties as described in Table 3.1. The single fracture version of the full-scale model shown in Figure 3.4 consists of a model with no-flow boundaries at the fractures edges and half interfracture spacing as shown (bounded by red dashed lines) in Figure 3.5. Figure 3.4 Horizontal well with 50 equally spaced fractures 52 Table 3.1 Summary of reservoir model and operational parameters Reservoir Top (ft): 9600 Matrix Permeability, kx, ky, (nD): 50, 100, 500, 1000 Matrix Permeability, kz (nD): xk * 0. Fracture Orientation YZ plane Number of Fractures 50 Initial Reservoir Pressure (psia): 4500 Rock Compressibility (1/psia) @5000 (psia) 4X10-6 Initial HC Saturation (%): 84 ( Single phase) Reservoir Porosity (%): 8 Flowing Bottom hole Pressure (psi): 500 Bubble Point Pressure (psia) 1965 Oil Gravity (API) 52 Reservoir Temperature (OF) 245 53 Proposed model 1 Single fracture model Figure 3.5 Aerial view of a full-scale horizontal well and its simplification approaches as presented by the single fracture and proposed models. Dashed lines represent no-flow boundaries for each model. From this model, it is straightforward to capture fracture interference through reservoir boundaries located exactly halfway between two neighboring fractures. By application of equations 3.1 and 3.2, the results of this model are modified in an attempt to represent a full-fledged horizontal well model. Note that the single fracture and the proposed models simulate only a portion of the full model. On the other hand, application of the proposed simplification workflow by following the steps described in the previous section results in a different model. Firstly, important phenomena, albeit not all (for the purpose of simplicity), are shown in Figure 3.3. Secondly, factors in the solution process diagram are crossed out on the basis of the model driving forces. Because this model is homogenous, with identical fractures spaced equally from one another and consists of a single well, components such as interwell interference, heterogeneity, and so forth can be crossed out. If the model is thin enough and/or gravity segregation is deemed negligible based on previous analysis, layering and layer-based completions can also be neglected under this basis. The final version of the solution process diagram is shown in Figure 3.6 (a). Lastly, the model is built based upon the solution process and symmetry application. According to Figure 3.6 (a), fracture interference needs to be included in the reservoir simulation as well as contribution from external and internal fractures. Clearly, under these requirements, the single fracture model is insufficient and a new model is necessary. Based on the model requirements and taking advantage of symmetry, "Proposed model 1" is built as shown in Figure 3.5. 54 Top/Bottom boundary Fortunes (a) (b) Figure 3.6 Solution process diagram simplification for a single well model (a) and multi­well model (b) As presented in "Proposed model 1", there is a fracture located at the far left of the model and another one located at the right at a regular fracture spacing distance. The left fracture is representative of interior fracture behavior while the right fracture represents exterior fracture behavior as required by the solution process diagram in Figure 3.6 (a). With this model, all the requirements in the solution process diagram are satisfied such as fracture interference and interior and exterior fracture behavior. Contrary to the single fracture model, this proposed model, result of the application of a new standardized simplification workflow, consists of two fractures. "Proposed model 2" is the result of the same solution process diagram shown in Figure 3.6 (a) if the "External fractures" component were also crossed out. Because several horizontal wells consist of a substantial number of hydraulic fractures, it may be feasible to ignore the contribution of external fractures to flow without sacrificing much accuracy. "Proposed model 2" also consists of two fractures at each side of the model, both serving the purpose of capturing the interference effect of internal fractures. Results of these three models (full-scale model, single fracture model, and proposed model 1 are discussed in the Results section. It is important to note that while results of the full-scale model are readily available, results from the simplification approaches need to be calculated either using equations 3.1, 3.2, or 3.3. The application of these equations is discussed as well for clarity. 3.4.2 Multiwell and Multifracture Case Similarly to the single well case, other more complex cases can be explored with this new simplification workflow. As a more thorough example, a three horizontal well model 55 is considered for simplification this time. Fluid as well as rock properties are found in Table 2.1 and the model schematics can be found in Figure 3.7. For the purpose of this example, the solution process diagram for the single well case shown in Figure 3.3 will be also used for this instance. After all important driving forces are taken into consideration and the solution diagram is simplified, the resulting diagram is shown in Figure 3.6 (b). Similarly to the single well case, a model is built by taking into consideration fracture interference effects as well as interwell interference. The resulting model is shown in Figure 3.8. As shown in the simplification schematics, the single fracture model does not change substantially while the proposed model adapts to capture well interference. Fractures elements are labeled as 1, 2, 3, and 4 which represent external fractures in external wells, 56 Figure 3.7 Full-scale multiwell simulation model schematic 57 Single fracture mode! Figure 3.8 Aerial view of a full-scale multi-well model and its simplification approaches as presented by the single fracture and proposed models. Dashed lines represent no-flow boundaries for each model. internal fractures in external wells, external fractures in internal wells, and internal fractures in internal wells, respectively. In order to clarify the objective of these fracture elements, one must consider the fluid flow contributions made by internal and external fractures in the single well case. By applying this concept to wells, the concept of internal and external wells is born. Hence, by application of symmetry, the proposed model is shown in Figure 3.8 and satisfies all components of its corresponding simplified solution process diagram. Note, that the proposed model only takes into account half a fracture for the internal well, and a full fracture of the external well. This was done with the purpose of capturing boundary effects as soon as external wells perceive them and to capture interference between wells and fractures while applying symmetry. 58 3.5 Results 3.5.1 Single Well Case After all three reservoir models (simplified single fracture, proposed simplified model 1, and a single well full-scale model) were run, the final step of the simplification process consists of scaling results. As mentioned before, the single fracture results are modified trivially as shown in equations 3.1 and 3.2. On the other hand, results modification from the "Proposed model 1" is not as obvious. Knowing that the number of fracture elements in "Proposed model 1" as shown in Figure 3.5 is 2, equation 3.3 becomes: Taking into consideration the fracture element labels 1 and 2 in Figure 3.5, Table 3.2 lists value assignments for variables in equation 3.4. Upon inspection of Figure 3.5, the symmetric factor for element 1 is 1 because the entire fracture element is modeled by the simplified model. Similarly, the symmetric factor for fracture element 2 is 2 because only half of the fracture behavior is considered in the simplified model. The multiplier factor represents the number of elements present in the full-scale model. Since element 1 represents external fractures, its multiplier factor is 2. Fracture element 2 represents internal fractures, hence its corresponding multiplier factor is the total number of fractures minus 2. Upon value substitution of equation 3.4 based on Table 3.2, equations 3.5 and 3.6 are used to determine modified oil rates and oil cumulative production, respectively. X = S1k 1x 1 + S2k 2x 2 3.4 Roil = 2 Rl ,oil + 2 R2,oil(P-f 2) 3.5 Np = 2Ni,p + 2N2,p(nf - 2) 3.6 59 Table 3.2 Symmetric and multiplier factors for elements 1 and 2 for the "Proposed model 1" 1 2 S 1 2 K 2 nf-2 Equation 3.7 shows typical cumulative gas oil ratio calculation according to data from Table 3.2. QQg _ ^0as _ ^1,5as + ^2,,gas(n/ - 2) 3 7 ^oii ^1,oii + ^2,oii(n/ - 2) Four cases for each model (single well full-scale model, single fracture model, and "Proposed model 1") were run with horizontal permeabilities varying from 50 nD to 1000 nD. Based on the above equations and the results obtained by numerical simulation, oil rates, cumulative production, and cumulative GOR were calculated and compared as shown in Figure 3.9 and Figure 3.10. As observed in Figure 3.9, cumulative oil production for all cases is almost identical for an initial period of time because during transient flow, fracture interference is not experienced in any model. However, as soon as boundary dominated flow and/or interference effects take place, there is a substantial difference between the full-scale model and the single fracture model. Since the proposed model does account for exterior fracture and interior fracture interference effects, it matches almost perfectly with the full-scale model. A similar situation is observed in oil rates as seen in Figure 3.10 where the quality of results from the single fracture model is even less relatable to the full-scale model. 60 (a) (b) Figure 3.9 Cumulative oil comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b) (a) (b) Figure 3.10 Oil rate comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b) One interesting feature observed for models considering internal and external fracture behavior is shown in Figure 3.11. Cumulative GOR for both fracture cases show upper and lower bounds which encapsulate the cumulative GOR behavior of the entire model. In other words, when the contribution of internal fracture is more significant than that of the external fractures, the entire model will show a GOR behavior close to the external fracture (lower bound). Similarly, when the contribution of external fractures is more significant, the entire model shows a GOR behavior closer to the internal fracture (upper bound) behavior. Hence, it is concluded that the full-scale GOR behavior cannot be higher than the internal fracture GOR or lower than external GOR behavior regardless of the number of fractures. This is proven as shown in Figure 3.12. 61 Figure 3.11 Internal and external fracture cumulative GOR comparison for a 50 nD case (a) and a 1000 nD case (b) 62 0 2000 4000 6000 0 Time (day) (a) (b) Figure 3.12 Cumulative GOR comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b) The qualitative difference between the full-scale model or the proposed model and the single fracture model are obvious from these results. A more quantitative approach was also applied to help determine the effectiveness of the proposed model which is product of the simplification workflow presented. Table 3.3 shows calculated Average Absolute Relative Error (AARE) and Normalized Root Mean Square Error (NRMSE) of both simplified models results when compared to full-scale results. The highest error calculated for the proposed model is about 2% for the GOR in the 50 nanodarcy case, whereas the highest error for the single fracture model is almost 100% for the single fracture case. It is again shown that error accrued by the single fracture model can be quite substantial. Table 3.4 shows run times for all three models proving that the new models retain valuable short times without sacrificing accuracy. Time (day) 63 Table 3.3 Comparison of simplified models to a full-scale model through statistical analysis for the single well case Reservoir Proposed Model Single Fracture Model Permeability (nD) NRMSE (%) AARE (%) NRMSE (%) AARE (%) 50 0.28 3.72 11.95 98.25 etea P=H Oil 100 0.00 0.06 12.49 98.01 500 0.52 2.79 17.43 97.96 1000 0.01 0.58 1.00 10.46 Oil e vi tialu l mu C 50 0.08 1.33 3.46 5.02 100 0.03 0.03 4.34 4.73 500 0.17 1.44 3.96 3.84 1000 0.20 0.08 2.73 2.98 R OGe vi tialu l mu C 50 2.25 1.03 15.43 7.32 100 0.77 0.90 6.25 5.83 500 1.14 0.83 7.63 4.31 1000 1.67 2.14 4.32 2.68 64 Table 3.4 Run time comparison for different models for the single well case Reservoir Permeability (nD) Run Time (Minutes) Single Fracture Model Proposed Model Full-Scale Model 50 1.6 1.6 17.9 100 1.4 1.5 17.8 500 1.4 1.5 17.9 1000 1.4 1.4 19.7 3.5.2 Multiwell and Multifracture Case Another advantage of this new technique is that it can be applied to well spacing. Well and fracture spacing studies are very important when it comes to oil production optimization. An accurate and simple simulation model is crucial to conduct timely studies. Even though full-scale reservoir models are desirable for the most accurate description of fluid flow behavior, multiwell fields are usually handled separately without considering well and fracture interference. The single fracture model's results from the last section can be used to calculate results for three wells. This is not the case for the proposed model which accounts for more fracture elements representing different components that contribute to fluid flow. Similarly to the previous section, symmetric and multiplier factors are determined using fracture elements as labeled in Figure 3.8 and shown in Table 3.5. Applying equation 3.3 to this case, the oil rate expression becomes: tfoil = 4<fa,oil + 4(nf - 2')^2,olI + 4(P-w - 2')^3,olI + 4(nf - tyfaw - 2')^4,olI 3.8 65 Table 3.5 Symmetric and multiplier factors for a multiple horizontal well case 1 2 3 4 S 1 2 2 4 K 4 2(nf-2) 2(nw-2) (nf-2 )(nw -2 ) Results comparison were made in the same fashion as the single well case as shown below. Cumulative oil production, oil rates, and cumulative GOR are shown in Figure 3.13, Figure 3.14, and Figure 3.15, respectively. Even though the error in the proposed model is slightly higher than the single well case, it is still significantly lower than the single fracture approach as shown qualitatively by plots and quantitatively in Table 3.6. Computational run times were also determined for each case as shown in Table 3.7. a) b) Figure 3.13 Cumulative oil production comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b) 66 (a) (b) Figure 3.14 Oil rate comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b) a) b) Figure 3.15 Cumulative GOR comparison for 50nD and 100 nD (a) 500nD and 1000 nD (b) 67 Table 3.6 Comparison of simplified models to a full-scale model through statistical analysis for the multiwell case Reservoir Proposed Model Single Fracture Model Permeability (nD) NRMSE (%) AARE (%) NRMSE (%) AARE (%) 50 2.48 2.95 2.18 21.43 etea P=H Oil 100 2.60 2.94 0.78 18.97 500 3.11 4.69 0.95 14.40 1000 3.37 6.96 0.98 16.89 Oil e vi tialu l mu C 50 0.92 3.85 3.74 4.86 100 0.99 4.27 5.39 5.27 500 1.32 5.55 1.71 2.26 1000 1.95 6.49 3.09 3.17 R OGe vi tialu l mu C 50 3.57 2.08 15.35 6.76 100 2.80 1.94 13.01 6.01 500 2.15 2.46 4.66 2.75 1000 1.67 2.14 4.32 2.68 68 Table 3.7 Run time comparison for different models for the multiwell case Reservoir Permeability (nD) Run Time ( Minutes) Single Fracture Model Proposed Model Full-Scale Model 50 1.6 2.7 81.3 100 1.4 2.7 77.7 500 1.4 2.4 74.7 1000 1.4 2.5 73.1 Again, it is shown that models resulting from the application of the proposed simplification workflow are reliably accurate and save valuable project time. Similarly, more complex systems can also be simplified by application of the present simplification workflow by accounting representative fracture elements and applying equation 3.3. This workflow is particularly helpful for well spacing optimization studies where modeling time and computational run times are impractical for full-scale models. 3.6 Fracture and Well Spacing Application One clear application of the simplification workflow is the study of well and fracture spacing. Several studies have addressed hydraulic fracture spacing with the aim to optimize profitability of hydrocarbon production. The fact that increased oil production is not directly proportional to the added number of wells due the drainage interference of wells is well recognized from the early research on vertical well spacing in conventional reservoirs.3311-3344 Zuber et al. (1995)35 and Baker et al. (2012)36 presented well and fracture spacing studies where simulation and economic analyses were conducted for natural gas coalbeds. Meehan (1995)37 conducted simulation studies to identify optimal fracture treatment designs and well spacing configurations for heterogeneous reservoirs. Meyer et al. (2010)38 presented approximate analytical production solutions for multiple patterned transverse hydraulic fractures where they looked at Net Present Values (NPV) as a function of number of fractures and propped fracture lengths. Hards et al. (20 1 3)39 used a fully compositional simulation to optimize fracture design, fracture spacing, and well spacing for the Cardium formation. Jin et al. (2013)40 estimated ultimate recovery based on correlations developed using fracture parameters to optimize fracture spacing in oil reservoirs. Eburi et al. (2014)41 looked at well interference effect on estimated ultimate recovery to seek optimum well spacing in liquid rich shale plays. Several studies were conducted to optimize well spacing configurations based on NPV in various fields such as the Eagle Ford42,43 and the Bakken.44 Integrity of well spacing with fracture half-length and number of fractures was studied in tight gas reservoirs.45,46,47 Mechanical properties in reservoir such as in-situ stress were considered in fracture spacing optimization48,49 by maximizing the fracture network between hydraulic fracture and natural fractures.50 Optimum fracture spacing of 200 ft. was found in gas-condensate reservoirs considering the Knudson flow through micropores.51 Fewer wells were drilled in liquid lean reservoir compared to liquid rich reservoir as shown in the simulation study on the effect of fluid compositions on well spacing for fixed hydraulic fracture geometry.52 Besides the deterministic approach using simulations, stochastic methods are also popular using empirical relationships such as decline curve analysis to study well spacing 69 optimization in oil reservoirs.53,54 In this section, a brief application of the simplification workflow is presented. After the workflow is applied, a resulting simplified model that represents a typical multiwell and multistage configuration is used for simulation and a brief spacing economic analysis is performed. The strengths of the workflow are shown in the form of accurate results and low simulation times for several simulation case studies where an economic assessment can be made quickly and accurately. For this case study, a 640 acre section is considered with fluid properties and operating conditions shown in Table 3.8. Fracture propped lengths are constant as they are not considered for this study. The two main variables for the simulation models are fracture spacing and well spacing. Due to geometrical considerations, fracture spacings of 40, 60, 120, 176, and 240 feet are considered. The number of wells in this scope are set to 1, 2, 4, 6, and 8. As shown in example 3.4.2, a multistage and multiwell system can be represented by 4 elements. Spacing between these 4 elements determines fracture and well spacing as can be seen in Figure 3.8. One set of simulations were run with all spacing combinations and a matrix permeability of 50nD. Two more sets were performed matrix permeabilities of 500 and 1000nD totaling 75 unique simulation cases. After simulations were run, results processing was done by application of equation 3.8 with symmetric and multiplier factors values from Table 3.5. It was clear that recovery factors are proportional to the number of fractures and number of wells. Hence, an economic analysis was carried out to determine the optimal configuration for this specific case. 70 71 Table 3.8 Spacing study reservoir simulation data Reservoir Top (ft): 12100 Matrix Permeability, kx, ky, (nD): 50, 500, 1000 Fracture half-length (ft): 305 Matrix Permeability, kz (nD): 0.1 * kx Initial Reservoir Pressure (psia): 4500 Rock Compressibility (1/psia) @5000 (psia) 4X10-6 Initial HC Saturation (%): 84 ( Single phase) Reservoir Porosity (%): 8 Flowing Bottom hole Pressure (psi): 500 Bubble Point Pressure (psia) 1965 Oil Gravity (API) 52 After 30 years of production, Net Present Values (NPV) for all 75 cases were calculated using capital and operating cost information found in Table 3.9. Surface functions of well and fracture spacing are plotted. Ideally, the optimum NPV is found at the maxima of the surface by numerical calculation of: dNPV(nw,nf ) dNPV(nw,nf ) 3.9 dnw ' dnf As shown by Figure 3.16, Figure 3.17, and Figure 3.18, net present values for this study were found to be indirectly proportional to fracture and well spacing for most cases. 72 Table 3.9 Economic analysis capital and operating costs Land acquisition ($/acre): 2,500 Permitting ($): 2,700 Site Construction ($): 400,000 Horizontal drilling cost ($/ft): 450 Vertical drilling cost ($/ft); 200 Cost per fracture ($/nf): 350,000 Completion ($/well): 200,000 Production to gathering station ($/well): 450,000 Royalty (%): 15 Discount Rate (%): 10 Depth, (ft): 12,100 Lateral length, (ft): 5,280 Oil Price ($/STB): 50 Gas Price ($/MSCF): 2.8 Monthly operating cost ($/well): 5,000 NPV ($MM) 73 Figure 3.16 Net present value for a 50 nD matrix permeability spacing study NPV($MM) 74 Figure 3.17 Net present value for a 500 nD matrix permeability spacing study NPV ($MM) 75 1400 1200 1000 800 600 400 200 0 ■ 0 0 FRACTURE SPACWG, FT NUMBER OF WELLS Figure 3.18 Net present value for a 1000 nD matrix permeability spacing study The most important finding is the definite impact of matrix permeability to fracture and well spacing. The relationship between NPV and fracture spacing change dramatically as fracture spacing becomes less and less important as permeability increases. This is due to interference being reached faster in high permeability reservoirs. Without any further calculation, it is clear from the surface plots that optimum economic configurations tend to be in the short spacing range. The total run time associated with the 75 cases studied in this section is less than a day, thus proving the efficiency of simplified models that result from the application of the present standardized workflow. 3.7 Key Findings A new standardized simplification workflow was presented and proven to greatly reduce simulation run times while achieving accurate results for production from low-permeability formations with hydraulically fractured wells. The simplification steps were explained in detail by accounting for phenomena that contribute to fluid flow, building a solution process diagram and crossing out phenomena deemed not relevant. A simulation model, which is a representative unit of a full-scale model and based on the simplified solution process, is built and its results are then modified to represent the full-scale behavior. It was shown that models resulting from this process require very low simulation run times while producing results that match full-scale model results. It was also shown that other simplified models may not be accurate representations of a full-scale system, thus proving the necessity for a new standardized simplification workflow. 76 CHAPTER 4 A NEW DISCRETE FRACTURE MODEL IMPLEMENTATION Natural fractures have been identified as important fluid flow drivers for most unconventional reservoirs. In most cases, history matching through reservoir simulation is not possible without due consideration of fracture contributions to the flow. In light of this, several attempts have emerged to include and properly represent natural and hydraulic fracture physics into numerical simulation. Amongst the most popular fracture representations found in the literature, three models seem to be standard for most engineering practices: single continuum, dual porosity, and the discrete fracture model. Even though studies have shown the advantages and disadvantages of some of these models,55 the relative low simulation run times and explicit representation of fracture networks make the discrete fracture model the preferred method. Implementations of the discrete fracture methodology have been made to finite element method (FEM) simulators in the last decades with promising results. A brief introduction to reservoir fractures from a geological point of view is followed by a discussion on the finite element numerical methods used in this research work. With this background information, the DFM implementation to FEM simulation is introduced and a new approach to the discrete fracture concept is proposed. Verification studies are performed on the newly developed approach by comparison to 78 analytical solutions and single continuum simulation results. The strengths, weaknesses and possible improvements to the suggested implementation are then discussed. 4.1 Background 4.1.1 Natural and Hydraulic Fractures Reservoir fractures can be defined as macroscopic planar discontinuities where a loss of rock cohesion has taken place through geological processes such as overburden or tectonic forces.56 Hydraulic fractures, on the other hand, are artificially created by injection of water into the rock, eventually causing the rock medium to crack. These fissures are then held open by proppant agents. Hydraulic fractures may be engineered to acquire certain penetration, half lengths, and widths within the reservoir and are generally depicted as shown in Figure 2.3 and Figure 3.1. Even though natural and hydraulic fractures are complex phenomena still under research, fractures in general can be described as ruptures in reservoir rocks. Most fractures can be characterized as faults and joints as shown in Figure 4.1. Faults are fractures along which one side has moved relative to rock on the other side.57 When no movement has occurred, the fractures are then known as joints. Faults can be categorized in one of two big groups: Dip slip and strike slip faults. Dip slip faults separate two rock blocks known as the hanging wall and footwall where the motion of the hanging wall relative to the footwall block occurs in a direction parallel to the dip of the fault plane. Depending on the movement direction of the hanging wall, inclined slip faults can be classified as normal faults or reverse faults as shown in Figure 4.2. 79 Figure 4.1 Fault and joint visual representation Figure 4.2 Schematics of dip slip and strike slip faults. Strike slip faults, on the other hand, do not have hanging walls or foot walls, this is because the motion of the pair of rock blocks occurs in a direction parallel to the strike line of the fault plane. This is also illustrated in Figure 4.2. The theory on fractures from a geological standpoint is too robust to be included in this section. Ample information in regard to the formation, classification, and evaluation of fractures can be found elsewhere. 57, 58 4.1.2 Fracture Representation in Numerical Simulation Fractures need to be characterized before consideration into numerical simulation is considered. Fracture information such as height, half-length, and conductivity is obtained by geophysicists and stimulation engineers. Seismic surveying at the time of hydraulically fracturing the target reservoir is one of the most trusted tools to gather information about the morphology and growth of fractures.59 Information of the injection schedule, proppant, and fracking fluid is also used to determine fracture size through simulation when microseismic mapping is not available. Because fracture permeability and width are hard to come by independently, fracture conductivity has become an important measure for fracture flow effectiveness. Fracture conductivity is defined as follows: Cf = wkj- 4.1 Equation 4.1 is usually expressed in field units of [md-ft]. To help understand fracture conductivity's relationship with the matrix, the dimensionless conductivity is introduced: wkj- 4.2 nkXf Cinco-Ley et al. (1978)60 showed that a dimensionless conductivity of 10 or more 80 reduces fracture pressure drop considerably. Hence, values for dimensionless conductivity at this range are considered essentially infinite. Non-Darcy effects are accounted for by correcting the dimensionless fracture conductivity as shown by Gidley (1991).61 Despite advances being made in the analysis of fractures, numerical simulation remains as the only way to handle complex fracture networks. However, fracture representations in the simulation space are varied and have been implemented in various type of simulators. For the most part, there are three common methods used to model fractures: 1. Single porosity model 2. Dual porosity model 3. Discrete Fracture Model (DFM) The following sections are dedicated to reviewing these models. 4.1.2.1 The Single Porosity Model This fracture model has been vastly used in the literature and consists of simply representing fractures explicitly. The fracture is meshed together with the matrix and the properties are given explicitly to the grid blocks hosting the fracture. Grid refinement around the fracture is necessary for convergence and can be distributed logarithmically or linearly as studied by Panja et al. (2014).62 A simple single porosity model is shown in Figure 4.3 The main disadvantage of the single porosity model is the large (sometimes enormous) number of grid blocks associated with the model that result in very long run times and large model sizes. Even though this problem can be somewhat alleviated with the introduction 81 82 Figure 4.3 Top view of a single porosity fracture model with linear grid refinement of Local Grid Refinement (LGR) techniques, other modeling techniques have emerged to address proper physical representation of the fracture while keeping run times relatively low. There is a second approach to the single porosity model where the permeability tensor for each grid block is modified to include the influence of fractures to fluid flow.63 The modified permeability is obtained using upscale methods as described by Oda (1985).64 However, this approach is mostly used to model short fractures as the characteristic length of the fracture is smaller than the characteristic length of the hosting grid-block. 4.1.2.2 The Dual Porosity Model The main concept of the dual porosity model is that matrix blocks and fractures are represented by two different continua. Fluid flow is carried through connecting fractures while the reservoir volume is represented by matrix blocks. A shape factor describes the connectivity of flow between matrix blocks and fractures. Mathematically, this model can 83 be expressed as: ft - + Qf 4.3 A schematic of the dual porosity concept is shown in Figure 4.4. This model was first introduced by Warren et al. (1963)65 to model natural fractures in single phase flow systems. Kazemi et al. (1976)66 later introduced a dual porosity application to a two phase immiscible system. The governing equations for fractures and matrix blocks are shown as follows: d f(pSv \ ( k r v \ = 4.4 d f y S p \ 4.5 ^ \ H / m = - q p 'mf matrix blocks fractures Figure 4.4 Dual porosity model representation 84 where qp,m/ in equations 4.4 and 4.5 is the matrix-fracture transfer function and can be calculated with equation 4.6. = ( k r n\ „ 4.6 where a is the shape factor which can be calculated as shown by equation 4.7 4 N (N + 2) 4 7 * = - f2------ N is the number of fracture normal sets. Calculation of I for N = 1, 2, and 3 is shown as follows: lx N = 1 4.8 2lJ y N = 2 I = lX + ly 3 lXly lZ JV - 3 v Ix^y + ly^z + where lx, ly and lz are spacings of fractures planes for each direction. The shape factor can also be calculated based on the Gilman-Kazemi67 formulation as shown in equation 4.9. ( 1 1 1 \ 4.9 ° = 4 \\ TLx2 + lLy2 + ll2z /l The concept of the shape factor, however, is controversial given the fact that a rigorous theoretical base for this concept is nonexistent. This remains one of the main drawbacks of the classic dual porosity model even though some work was made to address this particular issue.68,69,70 The original dual porosity model has other disadvantages that include the lack of gravity drainage, capillary continuity, reinfiltration, and viscous displacement. Fortunately, several researchers have proposed various improvements to account for gravity 85 segregation.71-73 Other research efforts went into the development of a subdomain model67,72-76 and a pseudo function method.73-81 A new implementation to the classic dual porosity model extends the original concept of matrix-to-fracture flow to incorporate matrix-to-matrix and fracture-to-fracture flow. This implementation is known as the dual porosity / dual permeability model and can be compared to the classic dual porosity model as shown in Figure 4.5. This formulation is governed by equations 4.10 and 4.11. The additional term in the matrix equation expresses the matrix-to-matrix flow portion of the dual porosity / dual permeability formulation. In fact, this additional term results in a model that requires greater computational effort than the classic dual porosity model. The limitations of the dual porosity model are well recognized for several applications. Because fractures are not modeled explicitly, hydraulic fracture modeling can become challenging. This problem is only accentuated by considering the lack of a rigorous basis for shape factors as previously discussed. For this reason, fracture interpretation through the discrete fracture model seems to be a solution. 4.10 4.11 86 Figure 4.5 Matrix-Fracture connectivity schematic for the (a) dual porosity model, (b) subdomain model and (c) dual porosity / dual permeability model. 4.1.2.3 The Discrete Fracture Model In the discrete fracture model formulation, the matrix is an n-dimensional domain that contains fractures represented by (n-1)-dimensional elements. For instance, a two­dimensional discretized reservoir contains one-dimensional fractures shown as lines as visualized in Figure 4.6. One of the earliest papers that used the DFM formulation to study fluid flow in a porous medium was published by Wilson et al. (1974).82 In this paper, they studied steady-state seepage in a fracture system beneath a dam. The first model consisted of an unstructured single porosity model where fractures are represented by triangular finite elements. In the second model, one-dimensional finite elements represent fractures in an impermeable medium. A two-dimensional model was implemented in a fractured medium using upstream weighted finite element method by Noorishad et al. (1982).83 Baca et al. (1984)84 followed 87 X Figure 4.6 Reservoir triangular element discretization around a fracture (red line) a similar approach to study two-dimensional single phase flow with heat and solute transport. Single phase approach was then extended by Bourbiaux et al. (1999)68 to multiphase flow, where a joint-element technique was used to represent fracture networks. Kim et al.55,85,86 and Yang (2003)87 used a similar approach as Noorishad et al. and Baca et al. to develop a two-phase black oil model with a parallel computing option. Karimi (2001)88 applied the same concept and developed an Implicit Pressure Explicit Saturation (IMPES) two-phase black oil model simulator. Fu (2007)89 extended this application to a three-dimensional fully implicit multiphase flow simulator where fractures are rep