Significant parameter identification and characterization of complex in situ reservoir simulations

The energy industry is facing major challenges as world energy demand continues to increase. Unconventional resources, including oil shale, have the capacity to meet that demand as production technology develops. There are many challenges associated with production of fuels from these resources including economics, environmental considerations, sustainability, and government policies. Improving predictive capabilities for energy production from these resources could play a major role in addressing these challenges. Simulation of thermal/reactive reservoir systems is complex. Heat transfers through the rock and fluids by conduction and convection, chemical reactions occur, phases change, multiphase fluids ow, and rock mechanical properties change. Each of these physical processes occurs simultaneously, affecting each other. Important physical processes also occur at widely varying length and time scales. Even where physical process models may be appropriate, the input data are highly uncertain. Useful modeling tools must find a balance between solution accuracy and computational effeciency. Simplifying assumptions are typically made according to data and experience. However, these simplifying assumptions may not be justified where data are sparse and where system characteristics change as a process unfolds. This research explores methods to expose the most important physical parameters and models for making efficient and useful predictions. Reservoir simulation methods for approaching thermal and reactive problems have been explored and analyzed showing heating by conduction alone is slow, and removal of products in a rubble system is not trivial. Experimental design methods have been implemented to expose some important parameters for predicting liquid fuel production, and surrogate models have been created with those parameters that approximate full simulation within 15% accuracy. Expected variations in kinetics and relative permeability modeling parameters predict a normal probability with a mean of 295 bbls (about 39 wt% initial kerogen in place) ultimate oil recovery with a standard deviation of 40 bbls (about 5 wt% initial kerogen in place). A case study involving a hybrid process involving in situ pyrolysis, in situ combustion, and CO2 enhanced oil recovery has shown that the energy input requirement for oil shale heating can be reduced by more than 335 million BTU, or 55%, and produce 160 bbls or 36% more oil. Finally, novel methods for performing these complex simulations are discussed. Progress and challenges with these novel methods are also discussed.

SIGNIFICANT PARAMETER IDENTIFICATION AND CHARACTERIZATION OF COMPLEX IN SITU RESERVOIR SIMULATIONS by Jacob H. Bauman A dissertation submitted to the faculty of The University of Utah in partial ful llment of the requirements for the degree of Doctor of Philosophy Department of Chemical Engineering The University of Utah August 2012 Copyright c Jacob H. Bauman 2012 All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of has been approved by the following supervisory committee members: , Chair Date Approved , Member Date Approved , Member Date Approved , Member Date Approved , Member Date Approved and by , Chair of the Department of and by Charles A. Wight, Dean of The Graduate School. Jacob H. Bauman Milind D. Deo May 8, 2012 Philip J. Smith May 9, 2012 John D. McLennan May 10, 2012 Raj K. Rajamani May 10, 2012 Eric G. Eddings May 10, 2012 JoAnn S. Lighty Chemical Engineering ABSTRACT The energy industry is facing major challenges as world energy demand continues to increase. Unconventional resources, including oil shale, have the capacity to meet that demand as production technology develops. There are many challenges associated with production of fuels from these resources including economics, environmental considerations, sustainability, and government policies. Improving predictive capabilities for energy produc- tion from these resources could play a major role in addressing these challenges. Simulation of thermal/reactive reservoir systems is complex. Heat transfers through the rock and uids by conduction and convection, chemical reactions occur, phases change, multiphase uids ow, and rock mechanical properties change. Each of these physical processes occurs simultaneously, a ecting each other. Important physical processes also occur at widely varying length and time scales. Even where physical process models may be appropriate, the input data are highly uncertain. Useful modeling tools must nd a balance between solution accuracy and computational e ciency. Simplifying assumptions are typically made according to data and experience. However, these simplifying assumptions may not be justi ed where data are sparse and where system characteristics change as a process unfolds. This research explores methods to expose the most important physical parameters and models for making e cient and useful predictions. Reservoir simulation methods for approaching thermal and reactive problems have been explored and analyzed showing heating by conduction alone is slow, and removal of products in a rubble system is not trivial. Experimental design methods have been implemented to expose some important parameters for predicting liquid fuel production, and surrogate models have been created with those parameters that approximate full simulation within 15% accuracy. Expected variations in kinetics and relative permeability modeling parameters predict a normal probability with a mean of 295 bbls (about 39 wt% initial kerogen in place) ultimate oil recovery with a standard deviation of 40 bbls (about 5 wt% initial kerogen in place). A case study involving a hybrid process involving in situ pyrolysis, in situ combustion, and CO2 enhanced oil recovery has shown that the energy input requirement for oil shale heating can be reduced by more than 335 million BTU, or 55%, and produce 160 bbls or 36% more oil. Finally, novel methods for performing these complex simulations are discussed. Progress and challenges with these novel methods are also discussed. iv For my wife, Carly. CONTENTS ABSTRACT : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : iii LIST OF FIGURES: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : viii LIST OF TABLES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : xi ACKNOWLEDGMENTS: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : xii CHAPTERS 1. INTRODUCTION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Oil Shale Processing Strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Ex Situ Technologies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 In Situ Technologies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Thermal Reservoir Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Research Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2. OIL SHALE RESERVOIR MODELING : : : : : : : : : : : : : : : : : : : : : : : : : : : 9 2.1 CMG STARS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.1 STARS Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Conductive Heating Process Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 Rubble Bed Process Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4 Conduction/Convection Heating Process Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3. PARAMETER SPACE REDUCTION AND PROXY MODELING : : : 41 3.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.1.1 Parameter Space Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.2 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3 Key Findings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4. COMBINED PYROLYSIS, IN SITU COMBUSTION, AND CO2 STORAGE PROCESS : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 56 4.1 Combined Process Simulation Background . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.2 Results and Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.3 Key Findings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5. THERMAL RESERVOIR MODELING DEVELOPMENT : : : : : : : : : : 74 5.1 Parameter Space Reduction Using Experimental Designs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.2 General Parameter Space Reduction With Hyperspace Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.3 Monte Carlo Sampling With Central Limit Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.4 Advanced Reactive Transport - Kinetics and Geomechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6. CONCLUSIONS : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 92 REFERENCES: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 95 vii LIST OF FIGURES 1.1 Green River Formation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.1 An aerial view of the heating and production well locations is shown on the left. The simulated wedge is shown on the right. . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 The weight percent of kerogen from the U059 well Fisher Assay data. . . . . . . 13 2.3 Cumulative oil and gas production. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4 Cumulative oil and gas production rates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.5 Energy supplied to reservoir and energy lost to under/over burden. . . . . . . . . 16 2.6 Temperature history comparison for three distances from heaters. . . . . . . . . . 18 2.7 Kerogen concentration comparison for three distances from heaters. . . . . . . . . 18 2.8 Oil saturation comparison for three distances from heaters. . . . . . . . . . . . . . . . 19 2.9 Coke concentration comparison for three distances from heaters. . . . . . . . . . . 19 2.10 Finite di erence grid representation of the Ecoshale process. . . . . . . . . . . . . . . 21 2.11 Heat distribution after 200 days from a horizontal heater. . . . . . . . . . . . . . . . . 22 2.12 Oil production result with 0.3 md rubble permeability. . . . . . . . . . . . . . . . . . . 22 2.13 Oil production result with an additional producer at the bottom of the capsule. 23 2.14 Oil production result with 10 md rubble permeability, including the additional producer at the bottom of the capsule. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.15 Relative permeability curves used in these simulations. The curves on the left endpoints have irreducible uid saturations while the curves on the right have endpoints extending to 0 and 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.16 Oil production result with relative permeability endpoints extended to 0 and 1. 26 2.17 Oil production result with reduced production pressure at 300 days. . . . . . . . . 26 2.18 Oil production result where inert gas is injected through heaters to drive oil to producers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.19 Best case scenario mass balance at the end of one year. . . . . . . . . . . . . . . . . . . 27 2.20 Temperature pro le of the rock after 7300 days (20 years) of heating. . . . . . . . 32 2.21 Best case scenario mass balance at the end of one year. . . . . . . . . . . . . . . . . . . 33 2.22 Oil cumulative production and rates results with varied kmul and with and without use of the STARS geomechanics module. . . . . . . . . . . . . . . . . . . . . . . 34 2.23 Permeability after 300 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.24 Permeability after 900 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.25 Permeability after 1800 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.26 Permeability after 6000 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.27 Permeability after 7300 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.28 Deformation after 300 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.29 Deformation after 900 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.30 Deformation after 1800 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.31 Deformation after 6000 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.32 Deformation after 7300 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.1 Normal distribution of activation energies for reaction 1. . . . . . . . . . . . . . . . . . 45 3.2 Oil/water and liquid/gas relative permeability curves. Sw is water saturation, and Sl is liquid saturation. Gas (krg), water (krw), oil/water (krow), and oil/gas (krog) relative permeability parameters are used in Stone's II relative permeability model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3 Eight run fractional factorial experimental design. . . . . . . . . . . . . . . . . . . . . . . 48 3.4 Pareto chart for parameter e ects on simulation time. . . . . . . . . . . . . . . . . . . . 49 3.5 Aerial view of simulated wedge. The distance between heating wells was reduced to 26.5 ft. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.6 Normal probability plot of the e ects on ultimate recovery of oil from 27􀀀4 fractional factorial design. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.7 Experimental design for 6 to 8 factors without confounding of individual parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.8 Pareto chart from 16 run fractional factorial design for 8 factors. . . . . . . . . . . 51 3.9 Histogram of Monte Carlo calculations of response surface. . . . . . . . . . . . . . . . 55 4.1 Horizontal Permeability after 400 days of in situ pyrolysis heating. . . . . . . . . . 61 4.2 Kerogen concentration after 400 days of in situ pyrolysis heating. . . . . . . . . . . 61 4.3 Coke concentration after 400 days of in situ pyrolysis heating. . . . . . . . . . . . . 62 4.4 Cumulative energy input for combined and pyrolysis only processes. . . . . . . . . 64 4.5 Oil production comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.6 Gas production comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.7 Oil pools in hot zones of the reservoir if heating and permeability generation is insu cient. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.8 Cumulative energy input showing results with additional increased pyrolysis heating rate case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.9 Cumulative energy input where pyrolysis is switched to in situ combustion at di erent times. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.10 Cumulative oil production where pyrolysis is switched to in situ combustion at di erent times. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.11 Layers based on Fischer Assay U059 core data. . . . . . . . . . . . . . . . . . . . . . . . . 70 4.12 Average richness uniformly dispersed throughout section. . . . . . . . . . . . . . . . . 70 ix 4.13 Disconnected kerogen rich layers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.14 Cumulative oil production results for three kerogen layer richness representations. 73 5.1 Models at the top of the pyramid are large scale highly coupled problems. Models at the bottom of the pyramid include small scale uncoupled problems. 79 5.2 Visualization of two parameter ranges subdivided into response surface subre- gions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.3 Arbitrary nonlinear three-dimensional surface. . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.4 Arbitrary nonlinear three-dimensional surface with one response surface. . . . . 84 5.5 Arbitrary nonlinear three-dimensional surface with four response surfaces. . . . 85 5.6 One-dimensional heat equation solution demonstration. . . . . . . . . . . . . . . . . . 87 5.7 One-dimensional heat equation solution combining the Central Limit Theorem with random Monte Carlo sampling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.8 One-dimensional heat equation solution combining the Central Limit Theorem with random Monte Carlo sampling in space while moving forward in time. . . 91 5.9 Histogram of sample mean set for a single calculated value at an arbitrary point and time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 x LIST OF TABLES 3.1 High and low molecular weights and associated stoichiometry for reactions 1 and 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.2 Range of activation energies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3 Uncoded parameters for screening design. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.4 Summary of calculated e ects. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.1 Representative component molecular weight . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.2 Combustion kinetic parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 ACKNOWLEDGMENTS I would like to acknowledge some important people who have helped me produce this dissertation. I appreciate help from my academic advisor, Dr. Milind Deo, who has given useful research guidance without sti ing creativity. I appreciate the kind natures my thesis committee members have shown as they have provided appropriate and constructive feedback. I appreciate the help and friendship of members that have come and gone from our research group. I thank my parents, especially my mother, for making my educational opportunities possible with great personal sacri ce. I thank my wife, Carly, for her patience and support throughout my time at the university. I will also thank my wife and my little daughter, Victoria, for making my life joyful in spite of life's changes and challenges. CHAPTER 1 INTRODUCTION 1.1 Motivation Sustainably meeting world energy demand while protecting and improving the environ- ment is truly the challenge of this age. Accessible energy has proven to be one of the most crucial drivers for improving the standard of living in developed, and in developing countries. However, this increasing energy demand has increased concerns about the world's ability to supply that demand. Nonrenewable fuel resources will not be sustainable inde nitely. Also, over the past decade, a great deal of attention has been given to the amount of greenhouse gases emitted into the atmosphere by the utilization of fossil fuels, and the uncertain potential e ects. Several organizations and governments throughout the world have mostly unsuccessfully attempted to nd alternatives to accepting the consequences of utilizing fossil fuels without harming economies or negatively a ecting the quality of life. Liquid fuels are especially valuable because of their high energy density. Liquid fuels can be relatively easily and safely transported and stored. Control technologies exist to reasonably deal with harmful criteria pollutants, though there are opportunities for signi cant improvements. The United States alone uses between 15 and 20 million barrels of oil each day, most of it being imported from other countries. Growing economies in some countries, like China, have caused the global demand for oil to rise, which has caused the price of oil to increase. In the summer of 2008, the price of oil rose to a maximum of about $140 per barrel. During this time interest in unconventional and alternative energy sources became urgent. Months later, however, market forces caused the price of oil to decline due to a worldwide recession. The need for sustainable energy sources was illustrated through this time period. Heavy oil, oil sands or tar sands, oil shale, and coal are examples of unconventional sources for liquid fuels. Each of these unconventional sources require some form of pro- cessing to produce liquid fuels. The amount of global recoverable liquid fuels is enormous when compared to conventional oil reserves. The Green River Formation alone, located in 2 Colorado, Utah, and Wyoming, contains an estimated 1.8 trillion barrels of high quality oil shale [22, 8]. Figure 1.1 is a map showing the Green River Formation. Similar scales of reserves of oil sands are located in Alberta, Canada, and production from these reserves have seen tremendous growth over the past 20 years. Coal resources are also massive throughout the world. Much of this coal could be suitable for liquid fuel production. This study focuses mostly on simulation of oil shale processes. However, principles learned during these studies should be applied to other thermal and reactive reservoir systems. Some other thermal and reactive reservoir systems of interest include geothermal reservoirs, underground nuclear waste storage, underground waste remediation, and carbon dioxide sequestration. Each of these applications is very unique and di erent, and this study primarily is focused on in situ oil shale reservoir simulation. However, modeling each of these diverse systems has some similarities. One similarity is they all include the solution of underground heat and mass transport equations. They all are highly complex and nonlinear systems. They all su er from the challenges of acquiring expensive and utilizing sparse data, and intricate heterogeneity. Solution of these problems is challenging, even with advances in high performance computing. One major part of these challenges is veri cation, validation, and uncertainty quanti cation. Su cient physical and modeling data are required to e ectively perform these tasks. Both physical and modeling data for these systems can be expensive and sparse. Large uncertainty often makes operations expensive and prohibitive. Reducing the uncertainty through modeling exercises, or through Figure 1.1. Green River Formation. 3 physical data acquisition is an important challenge. These e orts should be e cient for e ective applications and advances in technology and industry. There is an important need to improve physical understanding and modeling capabilities in thermal and reactive underground systems. 1.2 Oil Shale Processing Strategies The main hydrocarbon constituent in oil shale is called kerogen. Kerogen is an immature hydrocarbon that has not been converted to oil or gas through natural heating and pressure over a long period of time. Kerogen is a solid at room temperature, and is typically classi ed as not soluble in organic solvents. Production of liquid fuels from kerogen consists of heating the kerogen until the heavy hydrocarbon kerogen molecules crack into lighter hydrocarbon molecules, which are liquids or vapors. There are two principle options for oil shale processing into liquid fuel. These are ex situ and in situ. Ex situ refers to processes where oil shale is mined to the surface followed by heating. In situ processing refers to some form of heating while the oil shale is underground with subsequent production of the generated uids to the surface. Both ex situ and in situ processes have their bene ts and drawbacks. Ex situ processes require large mining footprints, a great deal of material handling, and large waste material volumes. The advantages of most ex situ strategies are greater control, improved e ciency, and higher energy yield. In situ processes have a reduced surface footprint, less material handling and less solid waste. In situ processes are more di cult to control, heating e ciency is challenging, and recovery of the generated products is not trivial. 1.2.1 Ex Situ Technologies Ex situ retort technologies heat mined oil shale in a variety of ways. Indirect heating methods are less common than direct gas, or solid heat carrier technologies. Indirect heating methods heat the oil shale through the walls of a pipe, kiln, or other retort vessel. One indirect heating technology discussed in greater detail in this research is the Redleaf EcoShale process. In this process, oil shale is mined and the rubble pieces of shale are buried in a containment capsule for groundwater isolation. Hot gases circulate through pipes installed in this containment capsule to heat the oil shale [5]. Direct heating methods heat the oil shale directly with hot gases or with solid heat carriers. Examples of hot gas heating technologies are Paraho or Shale Tech International Gas Combustion Retorts [4] and the Petrobras Petrosix retort [2]. Examples of solid heat carrier technologies include the 4 Alberta Taciuk Process [3] and the Ene t process [7]. Each of these processing technologies continues to evolve. 1.2.2 In Situ Technologies In situ technologies are appealing because land disturbance can be reduced when com- pared to ex situ mining operations. Spent shale wastes remain in place underground. Resources inaccessible to mining may be accessible by in situ heating and production methods. However, there are unique challenges associated with in situ processing. Some of these include groundwater isolation, external energy supply (for heating) constraints, energy e ciency, and recovery e ciency. Oil shale is typically impermeable and a poor heat conductor, making heating e ciency a crucial challenge. Shell oil has been active in developing in situ oil shale conversion technologies for more than 50 years. In 1971-1972 they conducted a pilot scale test in the Piceance Basin, a rich subsection of the Green River Formation, where hot water was injected into vertical wells to solution mine nahcolite, a mineral produced as baking soda or sodium bicarbonate. Dissolution of the nahcolite created permeability and porosity allowing steam injection for heating, and subsequent production pathways for generated uid products. The results of the pilot test were reported as generally disappointing in terms of production. There were two main problems noted by the pilot results. First, heating cause the shale to disaggregate into void spaces causing ow problems as the disaggregated material plugged permeable pathways. This disaggregation of oil shales is sometimes referred to as rubblization. Once the permeable pathways were plugged, steam could not be injected at design rates and temperatures. There were also issues with corrosion in the well piping [32]. In subsequent testing Shell moved to in situ electric heaters using knowledge gained through Swedish experience in the 1940s and 1950s with the Ljungstrom in situ electrother- mic method. The Swedish experience with electric heating reported 60 vol% Ficher Assay liquid recovery e ciency, and energy e ciency of 3.1 units of energy produced per 1 unit of energy supplied to the process [35]. The result of more decades of development has resulted in the current version of the Shell In Situ Conversion Process (ICP). The current version of the Shell ICP consists of a hexagonal pattern of in situ heaters surrounding production wells. The electric in situ heaters are spaced approximately 8 feet apart. No steam is injected. Heating depends on conduction through the reservoir rock, and is therefore slow heating of the kerogen. It takes months to heat the kerogen to pyrolysis temperatures. When the kerogen is converted to liquid and vapors it is believed that permeable pathways 5 are generated, allowing the pyrolysis products to be produced [49, 34]. ExxonMobil has developed the Electrofrac process. In this process hydraulic fractures from horizontal wells are injected with electrically resistive uid. When electricity is sup- plied, these planar fractures act as conductive heaters to the reservoir. Extending the heat source beyond the well into fractures increases the heating surface area. This makes heat transfer into the reservoir more e cient, and fewer wells must be drilled. Still, conductive heating through rock is quite slow, and after uids are generated from converted kerogen, production of those uids by multiphase ow through porous media is not trivial [43, 40]. American Shale Oil, or AMSO, owns another lease in Colorado for oil shale development. AMSO's in situ oil shale processing technology is called Conduction Convection and Re ux (CCR). Pairs of horizontal wells are drilled in an arrangement similar to steam assisted gravity drainage (SAGD) designs where a horizontal heating well is located some distance above a horizontal production well. The horizontal heating well heats the impermeable oil shale initially by conduction. It is believed that as solid kerogen converts and the pore space becomes unoccupied by it, permeable pathways will be generated. Also, heating will cause thermal compressive stresses in the rock leading to microfractures or mechanical failures, providing additional permeability. When this permeability is created, boiling oil near the heater rises as a hot vapor and heats the reservoir by convection. When the hot vapors are cooled far from the heater they condense and re ux until they are produced [16]. Chevron also owns a development lease in Colorado. Chevron has not been as aggressive advertising their methods and technology as the other lease owning companies developing in situ oil shale. Chevron's current planned strategies resemble eld experiments in the 1970s more closely than the conductive heating strategies described. Chevron plans to generate the required permeability for kerogen conversion and subsequent production of uids by placing explosives in the well bore. After the explosives have rubblized a pocket of oil shale rock, a proprietary chemical is injected into the reservoir to chemically convert the kerogen to oil and gas without a need for heating. Another alternative for generating permeability that Chevron is considering is cold uid injection causing the oil shale rock to fracture by thermal tensile failure [27]. Several other in situ technologies are being developed and advertised. Independent Energy Partners has developed an in situ geothermic fuel cell to generate electricity and fuels [6, 21]. Several companies including Schlumberger, Raytheon-CF, Phoenix-Wyoming, and Pyrophase have been developing radio frequency heating applications for oil shale. Radio 6 frequency heating technology is tuneable, which allows the process to focus heating energy on a target material. Radio frequency heating is also faster than heating by conduction [37, 21]. Mountain West Energy company is developing an in situ vapor extraction (IVE) process. This in situ process injects hot methane to heat the oil shale. The hot methane sweeps out the vaporized fuels, which are then condensed at the surface [23, 21]. 1.3 Thermal Reservoir Simulation Thermal and reactive reservoir simulation, in situ oil shale reservoir simulation for example, is very complex. Reservoir simulation predicts production of oil and gas at relevant time and length scales. Production in oil and gas wells occurs over the course of several years, and in over several feet to miles length scales. Physical phenomena at the scales of fracture widths, or pore spaces in the rock, need to be approximated with models. These models can be physical models or empirical models. In thermal and reactive reservoir simulation, the following physical phenomena are important, and need to be modeled: thermal conduction and convection, thermal geomechanics, chemical reaction kinetics, multiphase uid ow, and multicomponent phase behavior [42]. Each of these physical phenomena is capable of interacting with others. For example, in situ chemical reactions can occur over the course of years at reservoir temperatures, but at higher temperatures caused by reservoir heating by conduction or convection with injected or reservoir heat transfer uids, the same chemical reactions can occur in small fractions of a second. Another example is where heat transfer causes thermal stresses in the rock leading to failure and fracturing. Fluid ow through fractures is orders of magnitude di erent from uid ow in porous media, especially low permeability porous media like shales. Phase behavior depends on composition, temperature, volume, and pressure. Flow calculations depend on phase behavior, and uid phase properties. Phase changes can also have signi cant heat e ects. Modeling these complex interacting systems requires an e ective and e cient physics coupling strategy. Even after development of a fully coupled physical model that accurately solves each of these physical phenomena, quality data are required inputs to the simulators. Unfortunately input data in these types of reservoirs are expensive and sparse, if not totally unknown. Also, these data can vary signi cantly across the scales of interest. The rock and uid physical properties required as inputs for predictive calculations su er from a great deal of uncertainty and variability. Models incorporated into reservoir simulators are forced to 7 approximate physical properties depending on the desired knowledge to be gained from simulations, and expertise and experience with each physical case study. Statistical ap- proximations are often use to approximate unknown properties. Simulators are occasionally simpli ed to only incorporate physics that are critical for general insight into a system, or as screening tools [42]. In other words, many of the complexities in a process are ignored because the success of a process only depends on a few critical aspects. For example, in some oil shale resources experience may show that resources can be screened according to organic richness or presence of fractures without appreciating certain changes in mineralogy. In other resources, mineralogy may be a crucial indicator for engineering decisions and optimization. In development of a general purpose simulator, sensitivity of desired results to all physical inputs should be evaluated before the simulator is reduced to a screening tool. However, computational software and hardware are still limiting factors in the development of a truly general purpose simulator for thermal and reactive systems. Another development technique is to tune simple general models, or to tune case speci c models by history matching. Another major challenge with thermal and reactive reservoir simulation is computational e ciency. Representing all relevant physics at the time scales determined by important process dynamics, while only outputting data at time scales useful for decision making and optimization. These useful output time scales are not known because they depend on all possible variations of process dynamics. More importantly, ne resolution in spatial and temporal time scales is prohibitive when information over years and several miles is desired. Finally, though a simulator may be accurate, general, and e cient, the output results may will not be trivial to interpret. With hundreds or thousands of possibly interacting parameters throughout the course of a process, the many physical contributions to a result at any time may not be obvious. For example, in an oil shale reservoir, the production of liquid fuels should depend greatly on the richness of the resource, but if insu cient heat is supplied, and if permeable pathways do not develop or are plugged, the richness of the resource may have little impact on the amount of oil produced. Corollaries between physical contributions and calculated results can very easily be misinterpreted as cause and e ect relationships. The more intricate the relationships between physical contributions, the more di cult it is to sort through process dynamics and make conclusions about results. Small physical contributions at some particular moment during a process could have major implications years later. If the physics of that signi cant moment are considered unimportant because the same physics do not generally contribute signi cantly over the course of a process, then 8 poor conclusions are drawn by misinterpreting results. Consequently, current thermal reservoir simulators can be used to highlight physical trends, sensitivities, and general reservoir behaviors. Alternatively, thermal reservoir simu- lators can be developed for very specialized applications. Some current thermal reactive reservoir simulators being developed include STARS from Computer Modeling Group (CMG), Eclipse from Schlumberger, VIP-Therm from Landmark Halliburton, Finite Element Heat and Mass (FEHM) from Los Alamos National Labs, and MOOSE from Idaho National Labs. The Petroleum Research Center (PERC) at the University of Utah has also developed Ufes - CKT and Ufes - ARTS versions of thermal and reactive reservoir transport models. 1.4 Research Objective Ultimate goals for thermal and reactive reservoir simulation to address the challenges discussed should have awareness of the following problem characteristics. Several physical processes may have signi cant impacts on results at di erent times during the process. All important physical processes must be represented. Practicality and e ciency require justi ed assumptions to be made. Important physical processes occur at a wide range of time and length scales. Accurate physical insight into a complex process is important; not only accurate calculated results. This research gives a general description of CMG STARS, and the models and numerical methods used to solve these problems. Several processes have been modeled and analyzed with CMG STARS including a conductive heating process, a rubble bed process, and conductive/convective process. A methodology using experimental designs to create surro- gate models to more e ciently approximate full simulations, and to estimate uncertainty is demonstrated. A case study combining thermal reactive processes to address major challenges for in situ oil shale development is discussed. Finally, novel approaches for solving complex thermal and reactive reservoir problems are discussed, and steps for future development of these approaches are listed. CHAPTER 2 OIL SHALE RESERVOIR MODELING 2.1 CMG STARS CMG STARS is a reservoir simulator developed by CMG capable of modeling thermal and reactive processes. STARS has been used by modelers to study oil shale processes [36, 9]. The CMG STARS user manual lists the equations the simulator solves with a brief discussion on solution methods [1]. The conservation equations are summarized. Conductive heating, rubblized bed, and conductive/convective in situ technologies have been explored with STARS. Observations about simulating these processes are discussed. 2.1.1 STARS Equations Reservoir porosities and volumes are de ned as follows in STARS. The di erences between the de nitions of void volume or void porosity and uid volume or uid porosity are subtle but important. V = Vr + Vs + Vw + Vo + Vg (2.1) Vf = Vw + Vo + Vg (2.2) Vv = V 􀀀 Vr (2.3) v = Vv=V (2.4) f = Vf=V (2.5) Consequently, uid porosity is calculated with a relationship with solid concentrations and solid densities. These solids occupy a portion of the void porosity in the rock, but do not ow. f = v(1 􀀀 Xcsi si ) (2.6) Phase saturations are de ned as the volume of the uid phases divided by the uid volume. 10 Sw = Vw Vf = Vw fV (2.7) So = Vo Vf = Vo fV (2.8) Sg = Vg Vf = Vg fV (2.9) Sw + So + Sg = 1 (2.10) The mass and energy conservation equations are the principle equations that are solved in these thermal reservoir simulations. Accumulation of species is equal to the net in ow and net source terms in any reservoir volume. V @ @t [ f ( wSwwi + oSoxi + gSgyi) + vAdi] = ( wvwwi + ovoxi + gvgyi + wDwi wi + oDoi xi + wDgi yi) +( wqwkwi + oqokxi + gqgkyi) +(V (s0 ki 􀀀 ski)rk) + wqaqwk (2.11) V @ @t [ f ( wSwUw + oSoUo + gSgUg) + vcsUs + (1 􀀀 v)Ur] = ( wvwHw + ovoHo + gvgHg + K T) +( wqwkHw + oqokHo + gqgkHg) +(V Hrkrk) + (HLk + HLv + HLc) + (HACV + HACD)k (2.12) The volumetric ow rates for uid phases are calculated with Darcy's law. vj = T( krj jrj ) j (2.13) Flow and thermal transmissibilities, T and K respectively, and dispersion coe cients Dji , are calculated with geometric considerations accounted for and e ective parameters used between two calculated regions. E ective parameters are geometric or harmonic averages of permeability, thermal conductivity, and dispersion coe cients. For example, ow transmissibility is calculated as follows, where keff is the e ective permeability or area weighted harmonic average of the absolute permeability of two adjacent regions. 11 ( A l )effkeff (2.14) Solid components do not ow. The mass balance accumulation term for solid components is V @ @t [ vci]. Wells are treated as source and sink terms in the energy and mass conservation equations. Well phase rates are calculated with well index correlations, and are functions of the well pressure, and the pressure in the region containing the well. Finally, approximation models are used to calculate heat loss and aquifer source/sink terms. Phase mole fractions are related by K values. yi = Kgo i xi (2.15) xi = Kow i wi (2.16) wi = Kwg i yi (2.17) Phase mole fractions are also subject to the following constraints. yi = 1 (2.18) xi = 1 (2.19) wi = 1 (2.20) Sw + So + Sg = 1 (2.21) STARS is generally uses a nite di erence method for solving these partial di erential equations. Nonlinear equations are solved with Newton's method, and the linear system of equations with a general sparse linear system solver. An adaptive implicit solution method (AIMSOL) is used for e ciency. More details on numerical methods can be found in the STARS user guide [1]. 2.2 Conductive Heating Process Modeling Shell has developed an oil shale conversion process using electric heaters discussed in the introduction of this dissertation. STARS, a commercial reservoir simulator was used to simulate this process, and to gain understanding of the process characteristics, and simulation approaches. The well geometry based on the Shell ICP pilot where six heating wells surround one production well as shown in Figure 2.1. In initial simulations the heaters were spaced 53 feet apart. The thickness of the simulated reservoir was 50 feet. Due to 12 Figure 2.1. An aerial view of the heating and production well locations is shown on the left. The simulated wedge is shown on the right. symmetry, only a triangular wedge was simulated as shown in Figure 2.1. The simulated section can be repeated where reservoir properties are appropriate to estimate results for larger areas. Initial conditions need to be speci ed for these simulations. Initial conditions required in STARS include uid saturations and solid concentrations, reservoir temperature, and reservoir pressure. Gamma-ray log data from the Utah Geological Survey for the U059 well in the Uinta Basin [48] were used to estimate the weight percent of hydrocarbons (kerogen) initially in place. A kerogen rich section in the Mahogany zone of the well is from 665 feet to 715 feet deep, and the kerogen wt% varies from 12.5 wt% to 25 wt%. Figure 2.2 shows the weight percent of kerogen at di erent depths in the well which were used to calculate the initial kerogen volume, or solid concentration, at each depth. The remaining volume was assumed to be inorganic rock. The rock porosity in each layer was calculated with the assumption that kerogen nearly lled all the pore space in the rock, and that the inorganic rock density was constant throughout the simulated region. The initial pressure and temperature assigned to the region were a constant 1000 psi and 80 degrees F. In this simulation the reservoir was directly heated with two vertical injection wells to simulate resistive heating wells. The heaters heated uniformly from the top to the bottom of the well. The heaters each supplied 50,000 BTU/day to the reservoir for a four year time period. The production was pressure controlled by the producer well. This base case scenario used a bottom hole pressure (BHP) of 100 psi. 13 Figure 2.2. The weight percent of kerogen from the U059 well Fisher Assay data. 14 Representing the chemical kinetics of oil shale conversion to products has been studied by several researchers. Various conversion mechanisms have been proposed. Upon heating, kerogen is pyrolysed to produce oil, gas, and residue. Kerogen is a complex material, and the resulting pyrolysis products are also quite complex. Kerogen pyrolysis was simulated using properties of lumped representative components. The following reaction mechanism, adapted from a previous study [13] was used in these simulations. Kerogen ! HeavyOil + LightOil + Gas + CH4 + Char HeavyOil ! LightOil + Gas + CH4 + Char LightOil ! Gas + CH4 + Char Gas ! CH4 + Char Char ! CH4 + Gas + Coke All reactions are assumed to be rst order, and kinetic parameters in this simulation were taken from a previous study [13]. The heats of reaction were assumed to be 46.5 kJ/gmole for each reaction based on similar reactions from the template input les provided for STARS. Overall heats of reaction of kerogen pyrolysis to oil have been reported [18]. Stoichiometry was approximated based on the molecular weights and hydrogen to carbon ratios chosen for each component to force a mass balance. It should be noted that stoichiometric coe cients used in this reaction scheme are not unique, and are entirely dependent on chemical lumping assumptions. They are simply estimated to force mass and elemental balances based on approximated molecular weights and hydrogen to carbon ratios of each representative component. The results of this general simulation are shown in Figures 2.3 through 2.9. The cumulative oil and gas production over a four year period are plotted in 2.3. The production rates shown in Figure 2.4 show a maximum oil production rate of approximately 1.2 bbl/day. This maximum oil production rate occurs approximately two years after the heating is initiated. No signi cant oil is produced until after 400 days of heating. This time delay represents the time required to convert solid kerogen to producible oil with the given heating rate, well geometry, reservoir characteristics, and process parameters. These rates are multiplied by approximately 30 to convert the rate to bbl oil/day/surface acre. With the pyrolysis kinetic parameters and mechanism used in this simulation, much of the kerogen can convert to residue and gas at the local heating histories. Figure 2.5 shows the heat supplied and heat lost to overburden and underburden during the course of the process. After four years 50% of the heat supplied to the reservoir is lost to overburden and underburden. 15 Figure 2.3. Cumulative oil and gas production. Figure 2.4. Cumulative oil and gas production rates. 16 Figure 2.5. Energy supplied to reservoir and energy lost to under/over burden. 17 Processes could be engineered to minimize the heat lost to overburden and underburden by changing heating patterns, histories, and strategies. Figure 2.6, Figure 2.7, Figure 2.8, and Figure 2.9 show a comparison of three simulated grid blocks: one near the heater (block 18, 10, 11), one far from the heater (block 10, 1 11), and one in the middle of the simulated wedge (block 14, 5, 11). The reservoir was heated rapidly near the heater, and it can be seen that conduction through the reservoir is slow. Kerogen conversion is rapid at the high temperatures near the heater. The temperatures near the heater are unrealistically high, but high temperatures seem to be required for enough heat to conduct through a reservoir this size within a reasonable time period. The actual temperatures are dependent on the heat input strategy used in the simulation. It may take up to 700 days to supply su cient heat for pyrolysis far from the heater under these conditions. Coking can be signi cant with the given process parameters. 2.3 Rubble Bed Process Modeling Some of the major challenges with in situ oil shale processing are: e cient reservoir heating, generating permeable pathways for product uids to ow through, and isolating spent shale reservoirs from groundwater. RedLeaf Resources in Utah, for example, has developed a strategy to try to nd solutions to these challenges with their EcoShale process. A heating pit is constructed with heating pipes and an impermeable containment barrier, and mined oil shale rubble is dumped into the prepared heating capsule. The capsule is then sealed and buried. In a full scale process, the residual heat in a spent shale capsule is recovered and used in another capsule. This process may be a good option for addressing these challenges, however, mining and its associated challenges are not avoided, and the surface footprint could be extraordinary with a large scale process. Though the capsules are buried, subsequent land use would be limited. This process has some in situ characteristics, and STARS has been used to study these characteristics. In these simulations, the same kinetic mechanism is used where heavy hydrocarbons are thermally cracked into lighter hydrocarbons. This system can be di cult to model with nite di erence grid blocks. Finite di erence methods will not appreciate geometric com- plexities of the rubble pile. However, parameters possibly can be chosen that consistently approximate the overall capsule behavior. Also there are possibilities to gain understanding of truly in situ systems that are rubblized by explosives, solution mining of carbonates prior to heating, or other techniques where the oil shale is rubblized prior to heating. In this 18 Figure 2.6. Temperature history comparison for three distances from heaters. Figure 2.7. Kerogen concentration comparison for three distances from heaters. 19 Figure 2.8. Oil saturation comparison for three distances from heaters. Figure 2.9. Coke concentration comparison for three distances from heaters. 20 system heat transfers by conduction and convection to heat the outer surface of rubble pieces of rock. Then heat transfers by conduction through the rubble pieces, converting the kerogen to liquid and vapor that transports to the surface of the rock pieces, and then must ow through the rubble bed to be produced. Figure 2.10 shows how the geometry of a rubble bed can be discretized with a nite di erence reservoir simulator like STARS. In this gure the initial permeability of each grid block is shown. The blue blocks, which represent the rubble pieces of rock, have initial permeability of 10 md in this case, and the red blocks, representing the gaps between the rubble pieces, have an initial permeability of 1000 md. The blocks representing the rubble pieces are assigned properties that simulate oil shale rock characteristics, and the blocks representing gaps are assigned properties of some kind of vapor in a packed bed. Horizontal heaters placed in various locations add heat by conduction, as is done in the EcoShale and other processes, or by hot vapor injection. An example of the temperature distribution in the capsule due to the horizontal heating pipe is shown in Figure 2.11. In these simulations the target temperature for the capsule was 700 oF, and the initial temperature of the rubble bed was 220 oF. This initial temperature is high, but was used to reduced the simulation time required to gain some intuition from the results. The following set of gures, Figure 2.12 through Figure 2.18, illustrate how in complex thermal simulations, changing the required input parameters can signi cantly a ect the results, even when the parameters chosen are in a range that could be expected in a process. Figure 2.12 shows the liquid oil production results after one year of heating where four pairs of heaters and producer pipes are located several feet from each wall. The permeability with each rubble piece block was 0.3 md, and the Fisher Assay, or richness of the oil shale was 23.7 producible gallons per ton of rock. The production results should be multiplied by about 120 to predict the production for a pilot scale capsule. Still, the production in one year is insigni cant though production rate seems to be increasing. Upon inspection of the results, it seems that pressure management within the capsule is very important and challenging. Figure 2.13 shows how adding an additional producer at the bottom of the capsule to try to collect the liquids that may migrate to the bottom produces additional oil. Still insigni cant recovery of the potential oil is seen. The results in Figure 2.14 show that increasing the rubble piece permeability to 10 md has a major impact on the amount of oil produced in a year. Even though the oil shale is 21 Figure 2.10. Finite di erence grid representation of the Ecoshale process. 22 Figure 2.11. Heat distribution after 200 days from a horizontal heater. Figure 2.12. Oil production result with 0.3 md rubble permeability. 23 Figure 2.13. Oil production result with an additional producer at the bottom of the capsule. Figure 2.14. Oil production result with 10 md rubble permeability, including the additional producer at the bottom of the capsule. 24 rubblized into pieces, the generated liquids and vapors within the pieces must migrate to the surface of the pieces. Another source of uncertainty relates to the treatment of multiphase ow through porous media. Depending on the wetting characteristics of the porous media, and the interactions between phases, relative permeability models are used to approximate preferential ow between phases. Figure 2.15 shows the relative permeability curves used in these EcoShale simulations. In the relative permeability curves on the left, the endpoints do not extend to zero or 1.0 liquid saturation. This means that below a liquid saturation of about 0.20, liquids no longer ow. This is the irreducible liquid saturation. In the relative permeability curves on the right, these endpoint have been extended to zero and 1.0 meaning that any liquid oil that is generated is technically able to ow out of the rock pieces. However, at low liquid phase saturations, gases ow preferentially. Figure 2.16 shows the sensitivity of the result for oil production when these endpoints are extended to zero and 1.0. Extending the relative permeability endpoints seems to have a signi cant e ect on the recovery of oil after one year. Changing this parameter alone compared to the previous case more than doubles the amount of oil produced in the simulation. Inside the capsule and for in situ situations, pressure maintenance is the predominant factor for driving transport or ow. Figure 2.17 shows that when the production pressure is reduced at 300 days, the production rate of oil increases signi cantly. Also, this suggests that there is still oil present that is not being produced because the forces driving the ow of products are not e ective. Heating, phase changes, gravity, and production pipes all a ect the pressure gradients in the system. Managing pressure in ways that drive the ow toward the production pipes is ideal. Finally, if instead of conductive, and subsequent natural convection heating, hot gases were injected into the capsule which heats the capsule, more oil could be produced. This strategy heats the rubble bed more e ciently by forced convection, and the injected gases drive out the oil by displacing it. The simulated results of this approach are shown in Figure 2.18. Again, pressure management seems to be crucial for driving out the oil. A simple mass balance for this best case scenario is shown in Figure 2.19. At the end of one year the total hydrocarbon contents of the capsule are distributed into the following categories: unreacted kerogen, produced oil, unproduced oil, hydrocarbon gases, and residual carbon. With more simulation time, it was shown that more of the unproduced oil can in fact 25 Figure 2.15. Relative permeability curves used in these simulations. The curves on the left endpoints have irreducible uid saturations while the curves on the right have endpoints extending to 0 and 1. 26 Figure 2.16. Oil production result with relative permeability endpoints extended to 0 and 1. Figure 2.17. Oil production result with reduced production pressure at 300 days. 27 Figure 2.18. Oil production result where inert gas is injected through heaters to drive oil to producers. Figure 2.19. Best case scenario mass balance at the end of one year. 28 be produced, as expected since the production rates were increasing or comparatively high at the end of one year in most these simulations. Also, the unreacted kerogen is converted to pyrolysis products with more time and heating. If every pound of initial kerogen in place were converted to a medium API gravity oil were produced, the volume of produced oil would be about 27 bbls in the simulated section of the capsule. This approach to simulations has several limitations. For example, these one at a time simulations do not provide information on interactions between the parameters that were changed from simulation to simulation. All of these simulations showed results for one year of production, and additional time may have additional implications. Again, nite di erence representation does not account for geometric complexity, so properties given to the rubble blocks, and to the gap blocks are estimated, and those property dynamics are not taken into consideration. Also, these production results can be compared to eld data where available, and evidence of further limitations can be used for model improvement. On the other hand, this type of simulations can be useful for displaying trends, parameter sensitivities, and for providing intuition and optimization tools for engineering design. This exercise of simulating a rubble bed process illustrates several challenges with complex (multiphysical, dynamic, nonlinear) simulations. Any of these parameters that were explored in this series of simulations can be massaged to t available data, and this is often an approach used. Various strategies of history matching are examples of improving models by messaging parameters to historically collected data. History matching can be useful for cali- brating models, but the uncertainty in future predictions with those models can still be very large, even with a good t of historical data. Another approach to model improvement is by isolating physical phenomena from one another, conducting very controlled experiments, and then building a physical model valid for the experimental conditions. Later several physical models are coupled, and results from the coupled models should be consistent with data from the controlled experiments. Limitations with this approach include synergistic or diminutive relationships between interacting physical phenomena must be accounted for by coupling alone, the relative importance of each physical phenomena is not readily apparent, collecting experimental data for every possible scenario is expensive, and physical models do not extend beyond the bounds of the experiments. Also, where the models fail to make accurate predictions, little can be done to improve these models. So messaging parameters to t data does not necessarily increase predictive capability, and physics isolation followed by physics coupling can be ine cient and may not capture interactions between physical 29 phenomena appropriately. 2.4 Conduction/Convection Heating Process Modeling One of the potentially prohibitive challenges with in situ oil shale production is the required number of wells to supply adequate heating coverage in rich layers of rock. AMSO is developing technology to address this issue for oil shale similar to steam assisted gravity drainage (SAGD) for oil sands. Advances in horizontal drilling have allowed drillers to cover more resource area with fewer wells. The process initially heats the reservoir by conduction. As liquids and gases are generated, a re uxing zone is created where uids near the heater rise and condense when they reach cooler rock regions. Condensing vapors fall due to gravity into a producer at the bottom of the re uxing zone. The big question for this process is if permeable pathways will develop as the process unfolds, and that they will develop in a way that is expected by the process design. Permeability dynamics are a crucial part of modeling in situ oil shale pyrolysis, especially when heating is conductive. This type of process is particularly dependent on widespread permeability generation by conductive heating, so an understanding of in situ geomechanics due to thermal stresses is important for modeling the performance of this process. General observations about solid mechanics can be applied to oil shale systems. First, brittle mate- rials, like shales, may fracture when heated due to anisotropy and nonuniform dimensional changes. Polymers in general can expand or deform a great deal during heating, and they have low thermal conductivities. Kerogen should share some of these same properties of polymers. Increased porosity reduces heat conduction e ciency and convection within pores is ine ective. This suggests that as kerogen converts to uids and porosity increases, conduction becomes less e ective. Convection potential improves if uids are able to ow because of connectivity development between pores. Free expansion of solids is stress free. Constrained expansion in solids leads to compressive stress. Temperature gradients within a solid cause di erential dimensional changes with associated compressive and tensile stresses. Rocks are generally much stronger under compression than they are under tension. Pores, or a ductile phase, impede propagation of thermally induced cracks [17]. Any fracturing would be limited due to pores and kerogen phase in the solid. Mechanical behavior of oil shale under thermal treatment should is expected to vary signi cantly depending on the con ning conditions, the richness, the porosity, and the heterogeneity of the rock. Each of these parameters can signi cantly a ect the mechanical behavior, but each of these 30 parameters also varies signi cantly spatially and temporally in such in situ processes. Experimentalists have studied the mechanical behavior of oil shale. Tisot and Sohns [46] studied oil shale response to heat and stress, and the induced permeability. They heated oil shale plugs under a variety of conditions, and concluded that \... kerogen ... is the predominant contributor to [rich oil shales] properties and to their response to heat and stress." Still, the crucial need is to generate permeability. In this study, the researchers heated oil shale fragments under compression. Regarding permeability in the compressed fragments they conclude, \In most instances the induced permeability in the column of fragments was reduced to zero." After the kerogen support had converted, the rocks could not support any stress. Another key conclusion is, \This investigation shows that structural deformation in rich oil shales can be expected to occur ahead of the retorting zone ." Thomas studied structural behavior of oil shale with overburden pressure [45]. He concluded in this study that signi cant thermal fracturing does not occur when heating oil shale in an overburden environment. However, there was an increase in permeability that he attributes primarily to removal of oil water, and the decomposition of carbonates. Finally, Shell reported results from mechanical testing of oil shale in conjunction with eld trials in the 1970s [31]. In these early eld trials, nahcolite, a carbonate mineral found together with Green River oil shale, was solution mined with hot water, and the oil shale was later heated by steam injection. Tests showed that the heated shale releases thermal stress and pressure at the open faces the nahcolite previously occupied. Actual production of oil from the oil shale encountered further di culties in these tests, but the experience led to nahcolite mining industry improvements in the Piceance Basin. American Soda company found that thermomechanical fragmentation of oil shale propagated out more than 100 feet from the well [33]. These results, however, include a solution mining of nahcolite step, and not simply conductive heating. AMSO, ExxonMobil, and Shell have each commented on what they expect for oil shale mechanical response to heating. AMSO says, \The shale ... will want to expand as it is heated, but since it is con ned by the cool shale, it undergoes compressive failure and lls the high permeability conduit with rubble ... the thermomechanical fragmentation process is expected to propagate out to retort diameters of 100 or more feet." [16] Shell states, \... injection of hot water to leach the nahcolite and other salts ... was successful ... in generating the required permeability and porosity ... it was hypothesized that bulk heating with thermal conduction would generate permeability and that the gases generated during 31 retorting will drive liquid oil from the pores of the shale." [34] Finally, ExxonMobil boldly states, \... hydrocarbons will escape from heated oil shale even under in situ stress ... [Our] set of experiments clearly indicates that, even under conditions of overburden stress, the kerogen conversion and expulsion process creates porosity and permeability that was not present in the original oil shale." [40] Chevron, however, still plans to generate permeability with explosives or other methods [27]. The general consensus seems to believe that there is enough evidence that permeability can be generated by conductive heating alone causing thermal stress and compressive failure in the rock. An empirical model for relating permeability to uid porosity is available in STARS as follows in Equation 2.22. k = koe_kmul( 􀀀 o 1􀀀 o ) (2.22) The permeability in this model simply is dependent on the uid porosity and on an empirical multiplier. As solid kerogen degrades to uids, uid porosity increases, and consequently permeability increases. The opposite is true with coking reactions where uids convert to solids, decreasing uid porosity. Another method for calculating permeability changes due to heating in STARS, is their geomechanics module [1]. Stresses are calculated from the rock expansion while heating. This expansion can lead to surface heaving. When the kerogen and the uid support within the rocks disappear when products are produced, subsidence should occur, especially where data show that kerogen contributes signi cantly to the compressive strength of oil shale [46]. The conductive/convective heating process was modeled in STARS, especially for ex- ploring how the empirical model in Equation 2.22 and geomechanics module a ect the production results. In these simulations the horizontal heating well is located 50 feet above the production well. This is likely a greater distance between wells than the AMSO design, but is done to more e ciently heat the entire 100 foot thickness by conduction, since the extent and e ect of thermal fracturing or microfracturing is not well known. The richness of the 100 foot section was 25 gal/ton Fisher Assay. Again, all the pore space was lled with kerogen. Figure 2.20 shows the temperature pro les of the resource after 7300 days of heating. The kerogen concentration pro les correspond to the temperature pro les as expected in Figure 2.21. The cumulative liquid oil production and liquid oil production rates are shown for several cases in Figure 2.22. In this gure, the factor kmul in Equation 2.22 was varied as 1, 5, and 8. Where kmul = 5, simulations were run with and without 32 Figure 2.20. Temperature pro le of the rock after 7300 days (20 years) of heating. 33 Figure 2.21. Best case scenario mass balance at the end of one year. 34 Figure 2.22. Oil cumulative production and rates results with varied kmul and with and without use of the STARS geomechanics module. 35 using the geomechanics module in STARS. Mechanical and thermal data for Green River oil shale were supplied to the geomechanics module to make stress calculations. It can be seen that each of the simulations shows a similar trend. However, there are some di erences. With greater permeability generation, the production improved as could be expected. Using the geomechanics module is less clear what results can be expected. It appears that when kmul = 5, the production of oil is slightly better when the geomechanics module is not used. With heating, expansion of the rock is expected, but loss of kerogen and uid support could cause subsidence at the same time. The stress calculations from the geomechanics module a ect pore pressure calculations as well. The permeability evolution can be visualized as the process progresses in Figure 2.23 through Figure 2.27. The permeability increases in the hot zones due to heating expansion, kerogen conversion, and decreases due to coking. The maximum permeability achieved is about 10.7 md. The heaving and subsidence, or deformation of the rock is visualized in Figure 2.28 through Figure 2.32. The maximum heaving deformation appears to be approximately 0.41 ft as estimated by the geomechanics module. With these simulations and discussion it can be concluded that understanding geome- chanics in oil shale should be crucial for most in situ heating strategies. Permeability pathways may develop due to mechanical failure, or by some other mechanism during heating to retort temperatures in an in situ environment. Permeability dynamics can have a signi cant impact on simulated results when using an empirical model, or when calculating stresses and deformation with the STARS geomechanics module. 36 Figure 2.23. Permeability after 300 days. Figure 2.24. Permeability after 900 days. 37 Figure 2.25. Permeability after 1800 days. Figure 2.26. Permeability after 6000 days. 38 Figure 2.27. Permeability after 7300 days. Figure 2.28. Deformation after 300 days. 39 Figure 2.29. Deformation after 900 days. Figure 2.30. Deformation after 1800 days. 40 Figure 2.31. Deformation after 6000 days. Figure 2.32. Deformation after 7300 days. CHAPTER 3 PARAMETER SPACE REDUCTION AND PROXY MODELING 3.1 Background Subsurface thermal production strategies are complex as has been described in the introduction of this thesis. The sensitivity of simulated results to variations in physical models or parameters can vary signi cantly throughout the parameter hyperspace. This is because several physical processes occur simultaneously, are highly coupled, and occur with widely varying time and length scales. When the parameters the calculated results are most sensitive to can be exposed, surrogate models can be created to approximate simulated results with much simpler and much more e cient models. There are several methods for creating surrogate models for reservoir simulation. There are pros and cons associated with any of these methods [51]. Proxy-modeling techniques show a strong dependence on the complexity of the full model, the size of the parameter space, and the quality of the input data. It was found in that study that the type of proxy model did not have a major impact on the performance of the surrogate model, if each were supplied adequate data. There were challenges nding global optima using surrogate models. The surrogate models also performed better with problems less characterized by nonlinearity. Surrogate models, however, can still be used to e ciently show trends, to analyze parameter sensitivities, and to explore portions of interest in a parameter space. Modelers should use caution to understand the limitations of surrogate models for making predictions, or for optimization. Simulations for understanding the Shell ICP were built to develop this methodology. Details about the simulations are described in Chapter 2, with some di erences. The heating in these simulations was not allowed to cause excessive temperatures near the heating wells as it was in Chapter 2. Factorial experimental design methods are a simple way for generating polynomial surrogate models. Factorial experimental designs give experimenters and analysts e cient 42 tools to understand the impact parameters have on a response in a process. Unlike \one at a time" experiments, factorial designs allow the researcher to estimate interactions between parameters with fewer experiments, or full simulations in this research. These factorial designs allow researchers to evaluate many factors together. Experimental design methods primarily were developed for quality assurance purposes, but have been used in a wide variety of applications [11, 26]. These experimental design tools have also been applied to various oil reservoir studies [30]. A common experimental design is the 2k full factorial design. These designs test k factors at two levels for each factor, high levels and low levels. Each combination of high and low values of each factor is called a run. Full factorial designs require 2k runs to test every possible combination of high levels and low levels for k factors. When the number of factors is excessive or runs are expensive, fractional factorial designs are used for e ciency. In fractional factorial designs runs are selectively eliminated from full factorial designs with the assumption that higher order interactions are much less signi cant than individual factors without interactions. These designs are represented as 2k􀀀p fractional factorial designs. Fewer runs are required, but information about the signi cance of higher order interactions is confounded with information about individual parameters. 3.1.1 Parameter Space Description Experimental design and analysis methods are useful for comparing the sensitivity of a response due to variable input parameters, including their possibly signi cant interactions. Full simulations, on the other hand usually require single values to be assigned to parameters rather than a range of possible values for each parameter. These techniques were used to study certain areas in the parameter space in oil shale simulations. The reservoir description in these simulations resembles the Shell ICP as described in Chapter 2. The parameters of particular interest for study are: molecular weight of kerogen, activation energy for kerogen cracking in reaction 1, activation energy for heavy oil cracking in reaction 2, activation energy for light oil cracking in reaction 3, relative permeability representation, and reaction enthalpy. Each of these parameters is required for calculating the mass, energy, and momentum balances solved by the simulator. Activation energies are required for calculating the reaction rate term in the mass balance equation in Equations 2.11 and 2.12. Relative permeability is used to calculate the ow term with Darcy's law in Equation 2.13. Reaction enthalpy is incorporated in the energy balance equation in Equation 2.12. Ranges for each of these parameters were estimated from various literature data, inherent uncertainty, or are estimated to explore sensitivities. 43 The molecular structure of kerogen in largely unknown, and variable. The molecular weight of kerogen has been reported in ranges from about 3,000 [50] to 27,000 [10]. The sto- ichiometry in the chemical reactions and the initial concentration of kerogen are dependent on the choice for molecular weight of kerogen to conserve volume and mass for all simulation runs. Consequently, when the molecular weight of kerogen is changed between simulation runs, the stoichiometry of the reactions and the initial molar concentration of kerogen in the pore space must also be adjusted for mass and volume consistency. The range of molecular weight, and associated stoichiometry for reactions 1 and 2 are shown in Table 3.1. The values in Table 3.1 demonstrate that stoichiometry for such a kinetic mechanism depends on the properties (molecular weight and H/C ratio) of the pseudo components and is therefore \nonunique." Mass and elemental balances are important in these representations. Ranges for appropriate activation energies have been reported [38], and can vary signif- icantly depending on experimental methods, and kinetic analysis techniques. These ranges for activation energy are shown in Table 3.2. Studies have reported that activation energy for kerogen pyrolysis is most appropriately modeled with some distribution [38, 39], but it is uncertain how much impact di erent representations of activation energy have on the simulation results at large scales, and with multiphysical phenomena occurring. A normal distribution with 5 kJ/mole standard deviation is shown in Figure 3.1. Distribution of activation energies for kerogen pyrolysis is a complex function, but is sometimes repre- sented by the normal distribution [15]. A perfect activation energy distribution cannot be represented exactly in STARS, so discrete quantities, determined by integrating under the distribution curve, represent kerogen reacting with a speci ed activation energy according to the distribution. Relative permeability representations are often approximated in simulation, but such approximations may have signi cant implications. The range of relative permeability curves in this study are shown in Figure 3.2, the low level being more linear and the high level being curved. The shape of the relative permeability curves depends on the resource, the wetting characteristics of the rock, and the constituents present in the pore space. Finally, heat of reaction could play and important role in the heat transfer e ciency depending on the characteristics of the associated reactions. E cient heat transfer through an oil shale reservoir is crucial to any successful operation. Heat of reaction for oil shale pyrolysis has been studied [18], but it is not certain how much heat is lost to reaction compared to heat required to raise the rock temperature. 44 Table 3.1. High and low molecular weights and associated stoichiometry for reactions 1 and 2. kerogen heavy oil light oil gas methane char coke MW(+) 20000.55 424.49 152.03 52.01 16.04 12.60 14.55 MW(-) 2974.84 424.61 151.99 51.95 16.04 12.55 14.55 Formula(+) C1479H2220 C31.75H42.82 C11.19H17.51 C3.35H11.63 CH4 CH0.55 C1.19H0.32 Formula(-) C220H330 C31.76H42.81 C11.19H17.50 C3.35H11.62 CH4 CH0.53 C1.19H0.32 Stoic rxn 1(+) -1 37.29 13.86 25.03 17.06 38.71 0 Stoic rxn 1(-) -1 5.55 2.06 3.72 2.54 5.8 0 Stoic rxn 2(+) -1 2.18 0.06 0.03 7.13 0 Stoic rxn 2(-) -1 2.18 0.06 0.03 7.13 0 45 Table 3.2. Range of activation energies. Reaction Low Eact (kJ/mol) High Eact (kJ/mol) Kerogen cracking(1) 195 225 Heavy oil cracking(2) 208 260 Light oil cracking(3) 208 260/233 Gas cracking(4) 235 270 Figure 3.1. Normal distribution of activation energies for reaction 1. 46 Figure 3.2. Oil/water and liquid/gas relative permeability curves. Sw is water saturation, and Sl is liquid saturation. Gas (krg), water (krw), oil/water (krow), and oil/gas (krog) relative permeability parameters are used in Stone's II relative permeability model. 47 3.2 Results and Discussion The initial experimental design was a 27􀀀4 fractional factorial design. The eight run design for the initial 7 factors (excluding (8) heat of reaction) is shown in Figure 3.3. Each row represents a simulation run at the parameter levels speci ed in the run. The un-coded levels for these parameters are shown in Table 3.3. Details about each of these parameters have been described. The numerical performance for each of these runs was not equal, in terms of run time and convergence. Some of the runs had excessive time step cuts due to rapid change in gas saturation. The response chosen for these runs was simulation time in order to pinpoint the possible causes of these excessive time step reductions. Figure 3.4 is a Pareto chart displaying the impact each of these parameters have on the simulation time. Activation energy for reaction 3, or factor X4, had the greatest impact on the simulation time. After investigation, it appeared that simulation time increased signi cantly when the activation energy for reaction 3 was greater than the activation energy for reaction 4. This could be due to the combination of rapid gas creation coupled with high gas mobility causing rapid gas saturation changes. The high value for factor X4 was lowered to 233 kJ/mole as shown in Table 3.2, and no major di erences in simulation time were observed in subsequent runs. This exercise illustrates the value of and application of these experimental designs. Using the lowered value for factor X4, runs in the fractional factorial design were completed with ultimate recovery of oil as the output response. None of the simulations produced acceptable amounts of oil for such a process. Upon inspection it was found that oil generated from kerogen had inadequate mobility in lower temperature zones far from the heaters to ow to the producer. As a result, oil components had large residence times in the reservoir, and eventually converted further to gas and residual components. This result gives insight into the design of such a process, speci cally the spacing needed between wells for successful operation. If heating wells are drilled too far from producing wells the residence time of the oil in hot zones of the reservoir will be excessive, and these oils will convert to gasses or residual solids, signi cantly reducing or even prohibiting production of liquid oil. However, capital and operating costs increase with the number of wells drilled. Well spacing is a crucial design consideration for this process since excessive residence time of products in the reservoir and the cost of drilling wells are competing considerations for optimal process design. The initial dimensions of the simulated domain were changed to resolve this issue of 48 Figure 3.3. Eight run fractional factorial experimental design. Table 3.3. Uncoded parameters for screening design. Factor Physical parameter Low level (-) High level (+) X1 Molecular Weight 3000 20,000 X2 Eact rxn 1 195 225 X3 Eact rxn 2 208 260 X4 Eact rxn 3 208 260 X5 Eact rxn 4 235 270 X6 Eact distibution rxn 1 without with X7 Relative permeability linear curved 49 Figure 3.4. Pareto chart for parameter e ects on simulation time. excessive residence time of the oil. Reducing the spacing between these wells assures the whole reservoir was at a high enough temperature for adequate oil mobility. Figure 3.5 shows the modi ed dimension of the simulated wedge, changing the distance between heaters from 53 ft to 26.5 ft. The same fractional factorial design was used with ultimate recovery of oil as the response. The normal probability plot in Figure 3.6 illustrates the results of the runs. Normal probability plots, like Pareto charts, are useful for visualizing the signi cance of the e ects for each factor. Dominating e ects will appear as outlier points on a normal probability plot. The points of the e ects on this plot in Figure 3.6 are linear without outliers indicating that there is no evidence from these runs that any factors are dominant or insigni cant. With the 27􀀀4 fractional factorial design used, single factor e ects are confounded with pair interaction e ects and higher order interactions. Additional runs are necessary to isolate the e ects of individual parameter contributions from confounding with the e ects of higher order interactions. Further runs were done with a 16 run fractional factorial design for 6 to 8 factors. All 8 factors, including heat of reaction were tested with this design. The design used is shown in Figure 3.7 where factors E1 - E7 represent possible interactions between parameters, but individual factors are isolated from possible confounding with interactions. The results from these runs are displayed in a Pareto chart in Figure 3.8. It appears that 50 Figure 3.5. Aerial view of simulated wedge. The distance between heating wells was reduced to 26.5 ft. Figure 3.6. Normal probability plot of the e ects on ultimate recovery of oil from 27􀀀4 fractional factorial design. 51 Figure 3.7. Experimental design for 6 to 8 factors without confounding of individual parameters. Figure 3.8. Pareto chart from 16 run fractional factorial design for 8 factors. 52 the most signi cant factors are X2, X4, X6, and X7 along with higher order interactions between parameters, likely between these most signi cant factors. These factors are acti- vation energy for reaction 1, activation energy for reaction 3, activation energy distribution representation for kerogen conversion, and relative permeability representation. It appears that activation energy for reaction 4 and heat of reaction have the least impact on ultimate recover of oil. The results from these runs can be used in a 24 full factorial design without any additional runs. The data were regressed with the polynomial model shown in Equation 3.1, y = 0 + ixi + i;jxixj + i;j;kxixjxk + i;j;k;lxixjxkxl (3.1) where 0 = the intercept (global mean), = single and higher order interaction linear coe cients, and x = input variables. This polynomial model forms a multivariate surface called a response surface. The e ects are calculated by taking the di erence of the averages of the responses at high and at low levels of each factor, and for interactions between factors, and the coe cients are half of those e ects. The coe cients are summarized in Table 3.4. The model t the experimental output exactly because the response surface interpolated the data in this case. Although this is not a theoretical model and may have little physical signi cance, insight about the signi cance of each parameter in the explored ranges can be garnered. In the context of the process geometry, and process design, a modeler can better understand how kinetics should be represented. Typically higher order linear interaction e ects are assumed to be negligible and can be used to estimate error [30]. Expert opinion and knowledge is advantageous for estimating error, and elimination of terms in this model perhaps are not justi ed since this knowledge is unknown [30]. Three random validation simulations within the experimental space were run to estimate the quality of the response surface, the empirical regression model, compared to a STARS simulation. The di erence between the response surface approximations for ultimate oil recovery and STARS simulation results ranged from 3% to 15%. The quality of the response surface could be improved at the cost of more experimental runs, either by reducing the experimental space or by adding additional runs to estimate curvature due to nonlinearities when parameters are continuous. Alternative experimental designs could possibly provide more accurate response surfaces with comparable or fewer total runs, however many of these alternative designs require additional expert knowledge about the problem or unjusti ed assumptions. 53 Table 3.4. Summary of calculated e ects. Factors Intercept 294.6071 X1 -25.9846 X2 48.71381 X3 -8.88369 X4 -36.7407 X1X2 17.66644 X1X3 -19.9203 X1X4 7.082688 X2X3 14.89081 X2X4 -2.82519 X3X4 1.098062 X1X2X3 -5.35956 X1X2X4 0.575937 X1X3X4 -3.37056 X2X3X4 -0.70869 X1X2X3X4 2.191438 54 Monte Carlo sampling calculations were performed to characterize the response surface. Random values for each parameter with uniform distributions were chosen for each calcula- tion. A histogram of 80,000 Monte Carlo calculations is shown in Figure 3.9. The average value in these runs was 294.7 bbls oil with a standard deviation of 40.1 bbls oil. A normal distribution with these values is also shown in the gure for comparison. It appears the Monte Carlo results are slightly skewed to the right of a normal distribution. This exercise helps to quantify the e ects of variations in input parameters on the desired output. The shape of this distribution could be a ected by the response surface itself, the sampling locations for Monte Carlo simulation, or by the distributions assigned to each of the factors. 3.3 Key Findings Although results for oil shale simulations in this study are calculated with theoretical governing equations, the interplay within various parameters is not trivial due to competing physical phenomena. Combinations of parameters that expose possible competing phe- nomena can have signi cant numerical implications. Molecular representations for kerogen with associated stoichiometry, heat of reaction for kerogen decomposition, intermediate oil cracking reaction (reaction 2) activation energy, and continuing gas cracking (reaction 4) reaction activation energy are insigni cant in determining the ultimate recovery of oil at the scale simulated in this paper. Kerogen cracking (reaction 1) activation energy, relative permeability representation, oil cracking to gas (reaction 3) activation energy, and activation energy distribution representation have signi cant impacts on the ultimate recovery of oil in these simulations. Expert knowledge or similar studies including large scale physical experiments are important for estimating statistical error for developing validated surrogate models. Otherwise, more runs are necessary for improving these models quality for approximating simulator results. The interplay between various ow and kinetic parameters has been explored. Geome- chanical, heat transfer, and equilibrium parameters for example may also play signi cant roles at certain scales in production results for such complex reactive transport systems. Parameters from acceptable theoretical models can also be included in experimental designs to evaluate their impact on results and to include these parameters in constructing response surface approximations as illustrated in this chapter. Response surfaces can be characterized to quantify risk and uncertainty of simulations according to variation in input data. 55 Figure 3.9. Histogram of Monte Carlo calculations of response surface. CHAPTER 4 COMBINED PYROLYSIS, IN SITU COMBUSTION, AND CO2 STORAGE PROCESS Several in situ heating methods have been conceptualized and developed in applications for production of heavy oils, oil sands, underground coal, and oil shale. Some of these include variations of in situ pyrolysis, cyclic steam injection, steam assisted gravity drainage (SAGD), in situ combustion, underground gasi cation, and microwave heating. Reservoir system complexity makes each of these heating methods challenging. E cient and uniform heating is di cult to control. Reservoir geological characteristics and boundaries limit the applicability of any heating approach. Some in situ heating methods that may be appropriate for some resource may not be appropriate for another resource due to di erences in the depositional environment, even if the resources are chemically similar. Sometimes local environmental changes within a targeted resource can cause a production technology to become unsuitable as the process unfolds. This challenge can be addressed with combination processes that take advantage of the strengths of more than one process. For example, in some SAGD operations pressure maintenance and heating e ciency become economically prohibitive challenges as the process unfolds over time and mature steam chambers begin interacting with immature chambers. Researchers have conceptualized and conducted ex- periments combining the bene ts of SAGD and in situ combustion. Following SAGD steam chamber development with in situ combustion in experiments increased oil recovery 20% by mobilizing residual oil, and possibly would isolate a mature chamber from the rest of the reservoir [29]. Oil shale resources typically are characterized with very low initial permeability. Kero- gen, the primary organic component of oil shale, is an insoluble solid. Temperatures of 650 F to 1100 F are required for kerogen pyrolysis. In situ pyrolysis by conductive heating is ine cient because heat transfer by conduction through reservoir rock is slow. Some exper- 57 iments have shown that as the initially impermeable oil shale rock is heated by conduction and kerogen is converted to uids, permeable pathways are generated [47]. Greenhouse gas emissions associated with oil shale processing and shale oil use provide another challenge for oil shale development. The emissions associated with oil shale processing depend on the source of the input heat energy. Any combustion operations are tied to the resulting emissions. Flooding with CO2 for enhanced oil recovery (EOR) strategies have proven to signi cantly increase the pro tability of several production operations [14]. Flooding with CO2 provides a drive mechanism as the reservoir and residual oils interact with the injected CO2. Another advantage of CO2 EOR is underground storage and reduction of greenhouse gas emissions into the atmosphere. An attempt to combine the bene ts of in situ pyrolysis, in situ combustion, and CO2 EOR to address the challenges with fuel production from oil shale is examined in this paper. Initially the oil shale is heated by in situ pyrolysis to convert kerogen near the heat source to uids and to create some permeability pathways. When adequate permeable pathways have developed, conductive heating is terminated and air is injected to fuel in situ combustion. The mobile heating front provides heat more e ciently to the reservoir by coke and residual oil combustion without the cost of external heating. when su cient heat and permeable pathways have been developed, CO2 is injected to drive out the remaining oil and for storage. Oil shale development in the United States has gone through cycles in response to market conditions. Development projects were active in the 1970s and early 1980s, and interest has again surged in the last decade. In 1980 the United States O ce of Technology Assessment examined the potential of oil shale technologies for development. This report only examined two developing in situ technologies called true in situ (TIS) and modi ed in situ (MIS) strategies. With TIS strategies, explosives were placed down hole, or pumped into a hydraulic fracture, to further fracture the resource followed by some heating method, like in situ combustion. With MIS strategies, a portion of the resource was mined, and subsequently explosives were placed to fracture the remaining portion [28]. A U.S. Department of Energy funded eld project terminated in 1978 concluded that no TIS technology existed at the time "that [would] allow shale oil recovery on an economically feasible basis." The report concluded that a 10% void volume is required to sustain in situ combustion and provide adequate permeability for generated oil ow. The use of explosives down hole did not provide the required void volume, especially because of formation swelling due to heating [24]. The major technological concerns expressed in these reports relate to challenges in 58 generating adequate permeability for uid injection and oil mobility. 4.1 Combined Process Simulation Background Many of the characteristics for modeling a combined process are similar to the models discussed in Chapter 3. Studies from this research have explored general sensitivity of pre- dicted oil recovery to various simulation parameters for the Shell ICP [9]. The same pyrolysis kinetics mechanism was used in this study. A sequential reaction mechanism was used where kerogen cracks to lighter products, and those products were allowed to crack further at the high temperatures. The millions of chemical species and reactions were lumped into these representative components based on molecular weight. Other kerogen pyrolysis mechanisms and species lumping schemes with a wide range of representative complexity may be used to empirically t experimental data. In simulation the representative complexity of reactions and component lumping schemes should nd a balance with computational cost. It is important to use a reaction mechanism representative enough to make results useful, but excessive complexity can make numerical simulation noninformative and even prohibitive. Kerogen ! 37:29HeavyOil + 13:86LightOil + 25:03Gas + 17:06CH4 + 38:71Char HeavyOil ! 2:18LightOil + 0:059Gas + 0:026CH4 + 7:13Char LightOil ! 0:0024Gas + 3:30CH4 + 7:86Char Gas ! 2:84CH4 + 0:51Char Char ! 0:0097CH4 + 0:023Gas + 0:78Coke The representative kerogen species was further subdivided into seven fractions with distributed reaction rates in a manner similar to other studies [41, 9]. Studies have shown that kerogen pyrolysis rates are most appropriately modeled with some distribution [38, 39]. The molecular weights for the associated stoichiometry in each pyrolysis reaction, which must balance mass and elements, are shown in Table 4.1. Combustion reactions included all hydrocarbon species reacting to CO2 and H2O. The values for activation energy, frequency factor, and heat of combustion for these reactions are shown in Table 4.2. Gas combustion reactions are not included because gas ows ahead of the combustion front with the given conditions. The same sensitivity study [9] summarized in Chapter 3 showed that with the reported simulations there was a maximum well spacing distance where oil could be produced. The geometry used in this study is similar to the Shell ICP where six vertical heating/injection wells surround a vertical production well. Only portions of two heating/injection wells and one producer were simulated. The distance between the heaters/injectors was 26.5 feet. The thickness of the simulated section was 50 feet. 59 Table 4.1. Representative component molecular weight Component Molecular Weight Kerogen 20000.55 HeavyOil 424.49 LightOil 152.03 Gas 52.01 CH4 16.042 Char 12.60 Coke 14.55 Table 4.2. Combustion kinetic parameters Combustion Reactant EAct (BTU/lbmol) Frequency Factor Hc (BTU/lbmol) Kerogen 59,450 3.02e10 1.2525e7 Heavy Oil 59,450 3.02e10 1.2525e7 Light Oil 59,450 3.02e10 2.9075e6 Char 25,200 416.7 2.25e5 Coke 25,200 416.7 2.25e5 60 The purpose of the simulations in this study is to explore the fundamental characteristics and implications of a conceptual combination process including in situ oil shale pyrolysis, in situ combustion, and CO2 enhanced oil recovery and storage. The process begins with in situ pyrolysis to heat the rock and kerogen near the heating well, converting kerogen to products and coke, and generating permeability in the heated zone. At some point pyrolysis heating is replaced by air injection for in situ coke combustion to utilize the energy from residual oil and coke, provide a moving combustion front that more e ciently heats the rest of the reservoir, and allow the heating requirement to be self sustaining. Figure 4.1 through Figure 4.3 illustrate the physical characteristics of the process that make the combination of in situ pyrolysis and in situ combustion appropriate. These gures show a two-dimensional aerial view of horizontal permeability, kerogen concentration, and coke concentration in a middle layer in an oil shale reservoir after 400 days of pyrolysis heating. It can be observed that the permeability correlates with the concentrations of unpyrolyzed kerogen and coke concentrations due to these immobile solids occupying pore space and blocking permeable pathways. It is di cult to depend on conduction to heat zones far from the heating wells. Also the coke near the heaters is a fuel that can be utilized for supplying heat to the process, especially since su cient permeability generated due to kerogen pyrolysis near the heating wells allows for air injection. 4.2 Results and Observations The rst set of simulations in this study used data from the mahogany zone in the U059 well in the Uinta Basin in Utah to estimate richness shown in Figure 2.2. The depth of this rich section of the resource is from about 665 feet to 715 feet where this well is located. In the rst step of the process the heating wells rapidly brought the temperature of the rock near the well to target pyrolysis temperatures. After 600 days the external heating was terminated and air was injected to initiate in situ coke combustion. At 2000 days air injection was replaced with CO2 injection for enhanced oil recovery and partial storage of the CO2. A comparison simulation without in situ combustion is shown to illustrate the di erences of the combined process. At 600 days in the comparison simulation, conductive heating pyrolysis continued with a lower heating rate to ensure the temperatures near the heating well did not become unrealistically high. The results presented in this study compare fuels production, energy input, and CO2 emissions for several cases. Because all these values are relative to the size of the simulated 61 Figure 4.1. Horizontal Permeability after 400 days of in situ pyrolysis heating. Figure 4.2. Kerogen concentration after 400 days of in situ pyrolysis heating. 62 Figure 4.3. Coke concentration after 400 days of in situ pyrolysis heating. 63 domain, or the e ciency of the process, the results can be normalized by some standard. For example, production of oil numbers can be normalized by 0.3645 acre-ft, giving units of bbls per acre-ft. Energy input and CO2 emissions can be appropriately normalized by the number of barrels of oil produced giving units of MBTU per bbl oil produced and ft3 per bbl oil produced respectively, for example. The input energy savings for the combined pyrolysis/in situ combustion process when compared to the pyrolysis only process were 112 MBTU, or about 25% of the total energy required as shown in Figure 4.4. Although hydrocarbon products besides coke can be consumed by combustion, 157 bbls more oil were produced with the combination process. Oil and hydrocarbon gas production results are shown in Figure 4.5 and Figure 4.6 respectively. Signi cant additional hydrocarbon gases were also produced with the combined process. This is somewhat counterintuitive since gases would readily be consumed by combustion; however, as hydrocarbon gases are generated they may ow toward the production well ahead of the combustion front. The pyrolysis only process did not produce as much oil as was expected. Upon investigation it was discovered that a signi cant portion of the generated oil pooled at the bottom of the reservoir without being produced as shown in Figure 4.7. An explanation for why this may have happened is that not enough permeability was generated far from the heaters. Consequently, oil was only mobile in hot zones of the reservoir, and further experienced secondary cracking transformations in those hot zones. This problem is remedied when the heat transfer through the reservoir is su cient to establish permeable networks between generated hydrocarbon uids and production wells. This was achieved in these simulations by increasing the heating rates from the heating wells. However, the cost of heating can become excessive. With a higher heating rate in the pyrolysis only process, shown in Figure 4.8, 93 more barrels of oil were produced in the same time period in simulation. There is a tradeo between oil recovery and heating requirement in the pyrolysis only process. Another option is reducing the distance between wells. Costs increase signi cantly with additional wells, so reduction of well spacing distance could be prohibitive for such a conductive heating process. The energy savings and oil yield are dependent on the time where the pyrolysis stage is switched to the in situ combustion stage. Figure 4.9 and Figure 4.10 show the energy requirement for the heaters and the cumulative oil production respectively where the py- rolysis stage is switched to the in situ combustion stage on day 100, day 400, day 600, and day 800. It can be seen where in situ combustion began at 100 days, oil production was 64 Figure 4.4. Cumulative energy input for combined and pyrolysis only processes. Figure 4.5. Oil production comparison. 65 Figure 4.6. Gas production comparison. Figure 4.7. Oil pools in hot zones of the reservoir if heating and permeability generation is insu cient. 66 Figure 4.8. Cumulative energy input showing results with additional increased pyrolysis heating rate case. 67 Figure 4.9. Cumulative energy input where pyrolysis is switched to in situ combustion at di erent times. 68 Figure 4.10. Cumulative oil production where pyrolysis is switched to in situ combustion at di erent times. 69 severely reduced compared to the other cases. This is because the pyrolysis stage did not have su cient time to generate permeability a signi cant distance into the reservoir. The pyrolysis stage did not produce enough coke to signi cantly fuel in situ combustion, and therefore the combustion stage did not generate su cient heat to sustain kerogen conversion to uids reactions. When the pyrolysis stage is switched to in situ combustion at 400 days, 600 days, and 800 days, the oil production rates are slightly a ected, but the nal yield is the same for each case. This suggests that there is an optimum switching period between 100 days and 400 days where energy savings can be maximized while maximizing oil yield. The controlling factors for this optimization seem to be that the pyrolysis stage must be long enough to generate su cient permeability and coke to allow the in situ combustion stage to sustain the process, and the characteristics of pyrolysis compared to in situ combustion at some point have a signi cant energy cost without any bene t of increased production. Although there are signi cant energy savings with a combination pyrolysis/in situ com- bustion process when compared with an exclusively pyrolysis process, the combination process could emit signi cantly more CO2. All simulations began a CO2 injection stage for EOR beginning at 2000 days (refer to Figure 4.4 through Figure 4.6 and Figure 4.8). It was found in these simulations that additional oil recovery was marginal and most the CO2 injected was eventually produced. This is probably due to the very small distance between wells. The pyrolysis process followed by CO2 injection left a net 8040 ft3 in the reservoir section. The combination pyrolysis/in situ combustion process followed by CO2 injection produced a net 597,000 ft3 CO2 from the combustion. It must be noted that the true CO2 emissions for pyrolysis are highly dependent on the source of energy for heating [12]. At the high temperatures involved with in situ pyrolysis and in situ combustion, it is possible carbonate decomposition reactions play a signi cant role in permeability generation and additional CO2 emissions. Additional simulations were designed to incorporate miner- alogical data and carbonate decomposition kinetics published in studies with Green River oil shale [44, 19]. These simulations assumed an average 22 gal/ton Fischer Assay oil shale with 31% dolomite and 17% calcite composition. No signi cant permeability evolution or CO2 generation due to carbonate decomposition was observed with the simulated conditions. It is uncertain how much the geological heterogeneity a ects simulated results at these scales. Simulations estimating resource richness according to data from the U059 Uinta Basin well, an average richness in the 50 ft section, and disconnected strips of organic rich layers were compared. Figure 4.11 through Figure 4.13 show the three cases for kerogen 70 Figure 4.11. Layers based on Fischer Assay U059 core data. Figure 4.12. Average richness uniformly dispersed throughout section. 71 Figure 4.13. Disconnected kerogen rich layers. 72 distribution. All simulations were tested with in situ pyrolysis heating. It is observed that the results predicting cumulative oil produced in 600 days for each of the three scenarios are very similar as shown in Figure 4.14. At these scales, and with the boundary conditions used for the wells, the variations in layer properties (porosity, initial permeability, richness) do not signi cantly change oil recovery predictions. 4.3 Key Findings Oil shale development has potential because of the massive, secure, and accessible resources. However, signi cant challenges are associated with oil shale processing. In situ processing is preferred to avoid mining, access resources inaccessible to mining, and reduce land surface footprint. In situ processing is energy intensive, oil shale is typically impermeable, and oil shale processing can have signi cant CO2 emissions in addition to emissions from the produced fuels. This study explored the possibility of a sequentially combined in situ pyrolysis, in situ combustion, and CO2 EOR process to address these speci c challenges. A thermal reservoir simulator was used to evaluate the process and compare the results with a simulated process employing in situ pyrolysis alone for supplying heat. It was found that with the dimensions of this oil shale reservoir simulation, the combined in situ pyrolysis and in situ combustion process signi cantly reduces the energy requirement for heating. However, CO2 emissions from in situ combustion may be signi cantly larger than an exclusively pyrolysis heating process depending on the source of heat energy used for the pyrolysis heating requirement. With an exclusively pyrolysis heating process it is essential that permeable pathways are developed permitting oil to ow to a production well. This requires adequate heat input and su ciently rapid heat transfer through the reservoir rock. Injection of CO2 for EOR and storage did not improve oil recovery signi cantly. The net CO2 remaining in the reservoir was insigni cant compared to the CO2 produced during the process. Carbonate mineral decomposition did not contribute signi cantly to permeability generation or CO2 emissions. The level of detail for representing layer richness and vertical permeability does not a ect the predicted oil production with the simulated conditions in this study. 73 Figure 4.14. Cumulative oil production results for three kerogen layer richness represen- tations. CHAPTER 5 THERMAL RESERVOIR MODELING DEVELOPMENT Several aspects associated with thermal and reactive reservoir simulations have been highlighted in Chapters 2 through 4. Heat transfer modes like conduction of heat through reservoir rock, convection of heat due to owing hot uids, heating due to reactions, and latent heat transfer by phase changes of reservoir components have been illustrated with the simulator applications. Multiphase uid ow through porous reservoir rock has been explored by varying parameters in Darcy's Law equations and relative permeability relation- ships. Flow through rock fractures has not directly been the focus of these applications. The importance of geomechanics modeling has been illustrated, mostly because permeability and ow paths for uids production is crucial. Phase behavior of uids a ects heat and mass balance results. Chemical reactions have been represented with varying complexity. In fact, all of these physical processes have been modeled with a limited variety of complexity and scales in these applications. The following list are key points that can be learned from the simulated results presented in those chapters. Some of these points have previously been mentioned in Chapter 1. Several physical processes may have signi cant impacts on results at di erent times during the process. All important physical processes must be represented. Practicality and e ciency require justi ed assumptions to be made. Important physical processes occur at a wide range of time and length scales. Accurate physical insight into a complex process is important; not only accurate calculated results. Speci c key points related to the simulated in situ oil shale applications are as follows. 75 Heating reservoir rock by conduction is quite slow. Convective heating requires permeable pathways. Secondary reactions cracking producible liquid uids to residual or less valuable com- ponents occur when the residence time of products in hot reservoir zones is excessive. There may be limits for the maximum distance between wells due to heat transfer and permeability limitations. High temperatures may generate permeability by thermal fracturing or by dis occu- pation of pore volume as solid kerogen converts to uids. More e cient methods for in situ oil shale processes can be designed with technological advances and improved reservoir dynamics understanding. Beyond understanding in situ oil shale reservoirs through simulation, a part of this research was to improve and develop computational approaches to more e ciently under- stand highly coupled nonlinear reservoir modeling problems with widely varying time and length scales. It is illustrated through these exercises that as models become increasingly complex due to simultaneous interacting physical processes and due to large variations in time and length scales, the relevant causes (model param- eters, physical models, physical processes, assumptions, data incorporation, and uncertainty are examples) of calculated results become less apparent. However, it is unlikely that all parameters and all models are relevant at all scales throughout the course of such a process. Greater e ciency and physical insight from modeling would be possible where only essential physical processes are modeled whenever or wherever they are appropriate. A number of ways to approach this type of problem have been studied, including adaptive mesh re nement [20] and parameterization [25] for examples. Adaptive mesh re nement involves methods to create a ne computational mesh where process dynamics occur at small length scales, and to create a coarse computational mesh elsewhere. Examples in thermal reservoir simulation where these methods are appropriate include in situ combustion appli- cations where dynamics are better modeled with high resolution at the moving combustion front. These adaptive mesh re nement methods help to increase computational e ciency with only the length scale issues discussed throughout this thesis. Parameterization involves one way or another tabularizing highly dynamical physical models to avoid having to make 76 calculations with these models each time step. The cited example uses tabularized tie lines to avoid equation of state computations each time step for calculating phase behavior. 5.1 Parameter Space Reduction Using Experimental Designs In Chapter 3 this research demonstrated a method for reducing the parameters space required for making thermal reservoir system calculations using experimental designs. The methods showed how experimental designs can be used to identify the cause of numerical performance di erences in complex simulations. The methods also showed how selected kinetics and ow parameters can be used to predict the ultimate recovery of oil using response surfaces to approximate full simulations. Random sampling of the response surfaces can give estimates for uncertainty quanti cation due to variations in the parameters in the response surface. The following steps outline how to use these methods for complex simulations. Choose parameters of greatest importance for predicting an outcome based on expert knowledge. Choose the expected range or variations in the parameters selected. Choose an appropriate experimental design. The choice should minimize the number of runs required while providing the most information possible. Insigni cant parameters can be eliminated, and further runs using more informative experimental designs can be calculated. A response surface model is constructed by regressing the calculated data. Monte Carlo sampling of the response surface parameters is used to perform several calculations on the response surface. One of the bene ts of this approach is that very complex, high dimensional problems can be reduced to simple regression models with relatively few parameters as inputs. The parameters of greatest impact show researchers and developers what data are the most important. Reducing variation in these data will reduce the uncertainty in predictions the most. This allows data collectors to know where to make investments for collecting data where they are sparse and expensive. In summary, this method for evaluating complex 77 modeling e orts allows you to simplify the models by focussing on only the most essential parameters for predicting an outcome. Signi cant parameters and processes are identi ed. A method for estimating uncertainty on a simple response surface is e cient and based on the uncertainty in the input data. This method also has several limitations. First, the choice of parameters is based on expert knowledge. The limits of expert knowledge can signi cantly a ect the quality of the results. Expert knowledge is also required for choosing the limits of the chosen parameters of interest. The reason parameters of interest must be methodically chosen is because the number of runs required by experimental designs increases exponentially with the number of parameters. Several styles of experimental designs exist, but each run costs a full simulation, and runs can become excessive with more than 10 or so parameters. Expert knowledge must choose the parameters wisely. Second, the regression response surface will only approximate the full simulation models. It cannot improve the quality of the physical representation. The quality of the approximation depends on the size of the experimental space, the choice of parameters to focus on, and the form of the regression model. It is di cult to estimate the quality of these approximations without more full simulations. These issues tend to drive these methods toward over simplistic models. Third, calculations from response surfaces are black box calculations. Trends and dynamics during a process can not be observed. Fourth, lurking variables are treated as constants while developing the response surfaces, but may be more signi cant than judged during some possible stage of the process. Finally, and importantly, deterministic simulations give no real estimate of variance in the output data. Variance must be estimated or assumed with judgement. These limitations make it di cult to put as much credence into response surfaces as one would in full deterministic simulations. Experimental design methods, and response surface methodology are typically applied to physical systems experiments. They can be used for quality improvement, and for design. Physical experiments always have some variation within the experiments, which give estimates for variance, and a t-test can be applied to determine whether a signal due to variations in a test parameter is signi cant above the noise in the data. These methods are not so clearly applied to deterministic computer simulations because physical models already exist a priori. These methods can better be compared to interpolation or regression methods for high dimensional calculations. 78 5.2 General Parameter Space Reduction With Hyperspace Regression One of the limitations of the experimental design methods discussed is a limited number of parameters and their limits must be chosen based on expert knowledge. Can the approach be generalized to approximate the whole parameter space? Often industrial process models are large scale and multiphysical processes. The feeling in this research is that a top down approach should be taken when modeling such processes. Only physics relevant at the scales of interest need to be calculated. Figure 5.1 illustrates this concept. The top of this pyramid represents highly coupled physical models with large scale applications. An example could be a eld scale reservoir simulation model. More generally, the top of the pyramid would include general (mass and energy) conservation equations. As you move down the pyramid, the models become uncoupled. For example, an energy conservation equation at the top of the pyramid can be broken into separate conduction, convection, and radiation models. Continuing with the example, if you move further down the pyramid from the conduction model, the scale of the material particles with di erent conduction properties decreases. Continuing down the pyramid gives models for calculating heat conduction coe cients for all materials, even down to molecular models that give molecular explanations for material properties for calculating heat conduction. Note that if the system is very large, heat conduction may not be a useful heat transfer mode, so calculating everything beneath it in the pyramid would be useless. The challenge in building useful models is to know which paths, and how far down the pyramid one must go to accurately model the system at the scale of interest. Using the response surface methodology discussed in this chapter more generally may provide a way to do this. Again, with deterministic simulations as the runs in experimental designs, this can be interpreted as a multivariate interpolation or regression method. To keep the number of parameters manageable, the parameters are chosen at the highest levels of the pyramid in Figure 5.1. This allows fewer parameters to be studied at a time. In regions where these parameters are signi cant, additional models moving down the pyramid can be investigated again using experimental designs to develop multi dimensional regression or interpolation surfaces. A necessary part of using experimental designs in this manner is choosing the range for each parameter. The choice of the range for a parameter is based on expert knowledge, but should cover the entire region of interest. However, the regression model will be less 79 Figure 5.1. Models at the top of the pyramid are large scale highly coupled problems. Models at the bottom of the pyramid include small scale uncoupled problems. 80 accurate with lower resolution. Figure 5.2 illustrates how parameter space can be subdivided to achieve the appropriate resolution of the regression response surfaces. This is one of the key bene ts of these methods. The regression models should not have greater resolution than needed to improve the quality of the output. This is a characteristic in full simulations that is rarely known, and assumptions must be made to generate solutions that balance accuracy with computational e ciency. In Figure 5.2 it is illustrated how two parameters can be divided into several response surface models. The models cover the range of the two parameters. Each rectangle in the gure is a di erent shape, representing that di erent submodels cover di erent areas in the parameter space, based on their appropriateness or t, rather than somewhat arbitrary assumptions about resolution. Not only does resolution refer to time and space discretization, but to the scale of any other parameters or models within a simulation. These two concepts, top down physical approach and multivariate response surface regression or interpolation, can be generally applied to cover the entire simulation parameter space of interest. This idea was developed in Visual Basic Applications with Microsoft Excel. The code and Excel worksheets perform the following sequential tasks. Choose maximum limits for parameters in the high level equations that cover entire expected parameter space. Complete fractional factorial design calculation runs. Identify four parameters with the greatest e ect on the calculated results. Convert fractional factorial runs, or calculate additional runs to complete full factorial experimental design runs. Regress or interpolate results from full factorial runs. Check the center of the response surface versus an original calculation with all original parameters in the center. Decide if the error at the center is acceptable. If unacceptable, reduce the range of the variable with the most signi cant e ect in half and repeat all steps. If acceptable, store response surface and move to next range of variables. Then repeat all steps until entire parameter space has been covered. 81 Figure 5.2. Visualization of two parameter ranges subdivided into response surface subregions. 82 This concept has been demonstrated with an arbitrary nonlinear three-dimensional surface shown in Figure 5.3. A full factorial experimental design was used to calculate the interpolation points at the corners of the surface. Figure 5.4 shows the original surface and one response surface that covers the entire parameter space. It can be seen that the response surface does not approximate the original surface very accurately, especially in the center. This is a major limitation of surrogate modeling where the desire is to cover the entire parameter space, but the approximation may not be adequate. Using the approach discussed in this section, the parameter space can be subdivided into a set of response surfaces. Figure 5.5 shows a set of four response surfaces for approximating the same original surface. These four response surfaces cost ve more full calculations of the original function, but approximate the original surface much more accurately. For the very complex subsurface reactive ow simulations focussed on in this dissertation, experimental designs can be used to generate the set of response surfaces to cover the entire parameter space in a methodical and e cient way. Unjusti ed assumptions are not necessary. Some limitations may include the number of response surfaces required for good approximations may become excessive, the number of full simulations required to generate an adequate set my be expensive, the quality of the response surface results will never be better than the original calculations, an e cient method for sorting and calling the appropriate response surfaces from the generated set depending on the location in the parameter space during a calculation is required, only one outcome can be calculated with a set of response surfaces, and dynamics driving that outcome cannot be calculated with the response surfaces alone. 5.3 Monte Carlo Sampling With Central Limit Theorem Uncertainty based on the uncertainty in input data and uncertainty that propagates from successive calculations is often a missing component of deterministic model calculations. Monte Carlo methods are often used for uncertainty quanti cation in static calculations, but can be more generally applied by combining the Central Limit Theorem with the Monte Carlo calculations for dynamic system calculations. The Central Limit Theorem basically states that the distribution of sample means is normally distributed, with a mean equal to the population mean, and a standard deviation equal to the population standard deviation divided by the square root of t