APPENDIX H: modeling in situ oil shale extraction; a subpart of Project Oil Shale Pyrolysis and In Situ Modeling. Final Project Report. Reporting period: June 21, 2006-October 21, 2009

The objective of this work is to develop a highly-coupled, multi-phenomenon modeling capability based upon the Comsol® multi-physics software package for in situ oil shale extraction using various heating technologies, including conduction, DC resistive heating and radio frequency heating, various pushing fluids and geometries for their delivery, and production of both liquid and gas products from the deposit.

APPENDIX H Modeling In Situ Oil Shale Extraction A Subpart of Project Oil Shale Pyrolysis and In Situ Modeling Final Project Report Reporting period: June 21, 2006 to October 21, 2009 Terry A. Ring, Jon Wilkey and Mikhail Skliar Chemical Engineering Department University of Utah June 17, 2009 DOE Award Number: DE-FC26-06NT15569 Submitted by: Institute for Clean & Secure Energy 155 South 1452 East Room 380 Salt Lake City, UT 84112 H-1 Project Objectives The objective of this work is to develop a highly-coupled, multi-phenomenon modeling capability based upon the Comsol® multi-physics software package for in situ oil shale extraction using various heating technologies, including conduction, DC resistive heating and radio frequency heating, various pushing fluids and geometries for their delivery, and production of both liquid and gas products from the deposit. Project Outcomes A multi-physics model of in situ extraction of oil shale has been developed which couples fluid flow, mass transfer of multiple species, heat transfer and AC (RF) and DC heating of the deposit. All physical properties used in these model equations are functions of the local chemistry of the deposit and of local temperature. After overcoming significant numerical difficulties, a 2D slice consisting of a heating and a production well located 25 feet (7.62 m) apart has been simulated for up to 5 years. The 2D slice is a right triangle consisting of the smallest repeating unit of a hexagonal drill pattern. The model calculates the concentrations of kerogen, bitumen, oil and gas at all locations in the deposit; physical properties such as viscosity, permeability, heat capacity, thermal conductivity, electrical conductivity, dielectric constant, and loss tangent; and pressure, temperature, and thermal and pressure stresses in the deposit. The results show that a pusher fluid, a gas in this work, is necessary to move the oil to the production well, that thermallyinduced stresses do not induce fracture of the deposit, and that more uniform heating of the deposit by RF heating is beneficial to oil extraction. Presentations and Papers Jon Wilkey, “Multiphysics Modeling in situ Oil Shale Extraction,” Undergraduate Research Oportunities Presentation at the University of Utah, summer 2008. H-2 Introduction A conservative estimate of the total world in-place oil shale resources is 2.9 trillion barrels (Dyni, 2003) of which 2.0 trillion barrels, constituting the bulk of the oil shale resource worldwide in both quantity and quality, are in the western United States encompassing the Piceance Basin and the Uinta Basins. The Rand report (Bartis et al., 2005) puts the range of recovery at 500 billion to 1.1 trillion barrels depending on the percent recoverable and accessibility. There are six BLM oil shale Research Development and Demonstration (RD&D) leases in Colorado and Utah, five of which are in situ production. In situ processes are also being vigorously pursued by all the major energy companies (see Table 1). However, fundamental issues related to the kinetics of kerogen conversion to natural gas and light oil products and the production of the resulting oil require further multi-physics analysis to aid in situ extraction. In situ processing is a highly energyintensive process. Better energy utilization and efficiency will be necessary to make the extraction of this resource cost effective. Since water is scarce in this part of the United States, the use of large amounts of steam as a pushing fluid will be difficult, making water conservation another important aspect of in situ processing. Table 1. In situ Processing Methods Under Investigation by Major Oil Companies in Colorado’s Piceance Basin and Utah’s Uinta Basin. Data from (Parkinson, 2006). Company Raytheon and Technologies, Inc Chevron Shale Company Exxon Mobil Heating Method Pushing Fluid CF Radio Frequency heating Supercritical CO2 In situ Containment ? Oil Hot CO2 CO2 ? Either hot fluids or electric current directed into deposit Electrical resistance heaters and 3 phase AC electric heaters Pressurized hot water Hot fluids ? In-situ methane Ice-wall Steam Ice-wall Steam/water ? In-situ methane Gas In situ methane ? ? ? Shell Frontier Oil & Gas, Inc. Shell Frontier Oil & Gas, Inc EGL Resources, Inc. Phoenix Wyoming Petro Probe IEP Superheated steam or heat transfer fluid Microwaves Hot gases Waste heat from solid oxide fuel cell; fuel cell will use produced gases to make electricity for sale ? unknown The basic concept of in situ oil shale extraction is to 1) isolate the oil shale structure to be extracted, 2) sink various production wells to allow liquid and gaseous products to be removed from the oil shale structure, 3) sink various wells to provide access to the oil shale structure for various methods of heating and/or pushing fluids and 4) turn on the heating method and /or pushing fluids and remove product from the producing wells. Reservoir isolation may occur naturally due to impervious formations surrounding the oil shale deposit. Other options include the installation of isolation walls made from cement using naturally occurring faults or of ice walls (Parkinson, 2006). The production wells are placed in the appropriate locations so that heating profiles and pushing fluids drive the oil shale extraction products to the production wells. H-3 The placement geometry is different depending upon the liquid/gas product distribution and the heating methods used. For 3-phase AC reactance heating, a Texas 5-spot pattern or a triangular pattern is often used for the heating wells (Shell, 2006). With today’s directional drilling, however, even horizontal production wells are possible, e.g. Steam Assisted Gravity Drainage (SAGD) production of oils sands in western Canada. The modeling effort allows the various reservoir and well geometries to be modeled for heating by various methods, pushing with various fluids, and production of both liquid and gas products. The model for this work is a highly coupled multi-physics model run on Comsol® multi-physics software. Model Description The multi-physics model consists of the simultaneous solution of several equations with the coupling of the physical properties used in these equations (e.g. temperature, pressure and concentration) to the results computed from these equations. Table 2 lists the equations solved and shows the extent of equation coupling in the model. The most important equations used to model in situ extraction of oil shale are the chemical reactions that take kerogen to oil and gas. A simplified overall rate equation taken from Hubbard and Robinson1 k1 k2 Kerogen Bitumen Oil + gas [1] has been used in this model. This model assumes a sequential series mechanism with first order irreversible rate constants given by: k1 = 0.0706 [1/s] exp(-1.696e5[J/mol]/RgT) k2 = 115.673[1/s] exp(-5.677e4[J/mol]/RgT) The overall reaction has a temperature-dependent heat of reaction that, based on the mass of kerogen reacted, is given by2: Hrxn = (320.07*T[1/K]-114093)[J/kg] For all modeling conditions used, this reaction heat is endothermic above 357 K. Governing Equations The governing equations for the multi-physics simulation include: Darcy’s Law [2] where u is the fluid velocity vector, is the permeability, f is the fluid viscosity, is the porosity of the deposit, f is the fluid density, g is the acceleration due to gravity and ez is the unit vector in the direction of the gravity force. The gradient of pressure is determined from the methane gas 1 2 Hubbard, A.B. and Robinson, W.E., USBM Rpt. Inv. 4744(1950) ibid H-4 pressure at any location in the deposit. The presence of methane gas is due to kerogen pyrolysis (equation [1]) or to the flow of the pusher fluid. Table 2. Modeled Multi-Physics Phenomena and Equation Couplings. Equation Phenomenon/Medium Fluid flow Oil D’arcy’s Law Gas D’arcy’s Law Pusher Fluid D’arcy’s Law Mass Transfer Kerogen/Bitumen Oil Gas Heat Transfer Deposit Coupling 1) Viscosity as a function of temperature 2) Pressure driving flow is due to decomposition of bitumen to gas 3) Porosity is function of reaction conversion & convection 1) Viscosity as a function of temperature 2) Pressure driving flow is due to decomposition of bitumen to gas 3) Porosity is function of reaction conversion & convection 1) Viscosity as a function of temperature 2) Pressure driving flow is due to applied pressure & decomposition of bitumen to gas 3) Porosity is function of reaction conversion & convection Convective Mass Transfer Convective Mass Transfer 1) Reaction is function of temperature 2) Porosity is a function of reaction & convection 1) Velocity from D’arcy’s Law for convection 2) Reaction is function of temperature 3) Diffusion coefficient is function of temperature 1) Velocity from D’arcy’s Law for convection 2) Reaction is function of temperature 3) Diffusion coefficient is function of temperature Convective Mass Transfer Conduction/Convection Radiation Boundary Condition DC Resistance Heating Deposit RF Heating Deposit 1) Chemical reactions are heat source 2) Velocity from D’arcy’s Law (oil & gas) for convection 3) Thermal conductivity is function of local chemistry & temperature 4) Density is function of local chemistry & temperature 5) Heat capacity is function of local chemistry & temperature Electrical Conduction equation Electrical conductivity is function of local chemistry & temperature Wave equation Dielectric constant & electrical conductivity are functions of local chemistry & temperature H-5 Convective Mass Transfer [3] where Ci is the molar concentration of species i, Di is the diffusion coefficient for species i, and Ri is the reaction rate for species i. Convective/Conductive Heat Transfer [4] where l is the local density of the deposit, Cpl is the local heat capacity of the deposit, kl is the local thermal conductivity of the deposit, Qrxn is the heat of reaction. QRF, the RF heating rate, is given by QRF =2 f E2 o r tan [5] where f is the frequency, E is the electric field, o is the permittivity of free space, r is the relative dielectric constant and tan is the loss tangent. If a DC current is applied to the deposit, Joule heating replaces RF heating, given by QJoule= E2/ where is the electrical conductivity of the deposit. Electrical Conduction Equation [6] where V is the voltage and is the electrical conductivity of the deposit. RF Equation (Wave Equation) [7] where E is the amplitude of the wave’s electrical potential and k, the wave propagation constant, is given by [8] k= where j is the imaginary number, is the frequency angular frequency [ = 2 f] and is the electrical conductivity of the deposit. A coupled orthogonal magnetic field also governed by the wave equation is solved in this model. The relative magnetic permeability of the deposit is assumed to be 1.0 for this model since no other data was available. Physical Properties Most of the local physical properties are determined from concentrations of rock, kerogen, bitumen, oil and gas that occur locally in the deposit at any given time using a simple molar mixing rule. Each individual property is given by a temperature dependent function that is determined from either first principles or from a best fit of temperature dependent data (e.g. H-6 viscosity of oil) for the material3. Since the pressures may be large, the compressibility of oil is considered and a virial coefficient equation of state is used for the gas (assumed to be methane). Electrical properties for the deposit, including electrical conductivity, relative dielectric constant and tan , are available as a function of frequency and temperature but not in a clean form as the data does not account for the temperature-dependent chemistry of the deposit. As a result, average properties are used for the electrical properties of the deposit except for the electrical conductivity, which is assumed to be a function of carbonaceous residue left behind after the kerogen decomposition. Mechanical properties of the oil shale as a function of temperature are not readily available. As a result, the mechanical properties at room temperature are used in this model to determine the failure strength, Young’s modulus ( ), Poisson’s ratio ( ) and the thermal expansion coefficient ( ). Thermally induced stress is generated elastically due to the temperature difference from the initial temperature of the deposit, To, as calculated by4: thermal= (To –T)/(1- ) [9] Model Geometry The system investigated consists of one-half of a 2D slice of a triangular drill pattern as shown in Figure 1. The model geometry consists of a heating well at one end and a production well at the other end of the 25 ft (7.62 m) hypotenuse of a right triangle. The production wells consist of concave, rounded off surfaces through which the heat flux and the production fluxes flow. To facilitate numerical stability, the right angle has also been rounded off to a convex surface. A view of the computational mesh with 9393 elements for this model geometry is shown in Figure 2. 3 Rajeswar, K., Nottenburg, R. and Dubow, J., Review Thermophysical Properties of Oil Shale,” J. Mat. Sci. 14, 2025-52(1979). 4 Ring, T.A., “Fundamentals of Ceramic Powder Processing and Synthesis,” Academic Press, 1999. H-7 Hexagon heating drill pattern with production well at center used as representative extraction geometry symmetry simplifies geometry Boundary Conditions Symmetry Constant Temp. Convective Flux Figure 1. Model geometry. Figure 2. Mesh for model geometry (scale in meters). Initial and Boundary Conditions Initially the oil shale deposit is filled with kerogen at 9% weight (12 gal/ton Fischer assay) with no bitumen, oil or gas present. Kerogen content varies in deposits from low values to high values; values typical for economically viable commercial operations are in the ~25 gal/ton range and H-8 higher. We have used 12 gal/ton as an initial condition for the deposit because higher values gave numerical instabilities in the model. When the kerogen decomposes, it is assumed that 63% becomes oil and 24% becomes gas with the balance being a carbonaceous deposit left behind. Initially, the deposit has a temperature of 400K and a down hole pressure of 10 atm. The boundary conditions are symmetry conditions for all of the boundary walls except the heating (lower left corner in Figure 2) and production (upper corner in Figure 2) wells. At the heating well (when used), a time dependent heat flux is used to soften the boundary conditions; the heat flux increases to 500 W/m2 over a time period of 10 hrs using an exponentially rising function. In other cases, a constant wall temperature at the heating well is used as a boundary condition. In these cases, the wall temperature increases from 400K to 1000K over a time period of 10 hours using an exponentially rising function to soften the boundary condition. When a DC electrical current is used for heating, it is passed between the heating and production wells. The voltage on the heating well surface is increased to 1,000 Volts over a time period of 10 hours using an exponentially rising function, and the production well is grounded. When an AC (RF) electrical signal is used for heating, the antenna is a dipole consisting of the heating and production wells with an external current of 2,200 A/m2. When a pusher fluid is used, the pusher fluid, methane, is applied at the heating well, with pressure increasing to 100 atm over a time period of 10 hours using an exponentially rising function. Computational Platform Simulations were performed on a Windows XP64 quad core computer with four 18.6 GHz Intel Xeon E532 CPUs and 10 GB RAM. Calculations on the 9393 element grid with 76,536 degrees of freedom were performed with a relative tolerance of 0.01 and an absolute tolerance of 0.001. The computational times for simulating 5 years of heating were typically from several days to a week. Results and Discussion Model results are given for a series of cases where the heating is done by thermal conduction, electrical conduction or electrical induction both with and without a pusher fluid. Conductive Heating with and without Pusher Fluid With the wall of the heating well being heated by a 500 W/m2 heat flux, the temperature of the deposit comes up to temperature slowly over a period that exceeds 5 years (43,800 hrs). The temperature (surface color), pressure (contours) and velocity (vectors) at the end of a 5-year heating period without pusher fluid are shown in Figure 3. The temperature has reached almost 1000K near the heating well and progressed through the deposit, but temperatures above 600K have not arrived at the production well. Pressure contours are barely visible and cover a very narrow range around that of the initial deposit pressure. The velocity vectors are very small. H-9 Figure 3. Temperature, pressure and velocity profiles at the end of a 5-year heating cycle with heating due to thermal conduction. Without pusher fluid, the amount of oil that has left the deposit is negligible as shown in Figure 4. With pusher fluid, the temperature profile is not significantly changed from that shown in Figure 3, but the amount of oil that has left the deposit is significantly larger as seen in Figure 4. This result suggests that convection is not significant in moving the heat around the deposit as the velocity vectors are small for both cases. In the case with the pusher fluid, the oil is being pushed ahead of the pusher fluid, albeit slowly. The amount of gas that has left the deposit is also negligible over the 5 year time period with and without pusher fluid; both cases provide only 10-8 weight fraction of gas at the production well (not shown). The deposit in the region of the production well in both cases is essentially plugged against gas flow by kerogen and oil filling the pores in the deposit, thus severely reducing gas permeability. H-10 Figure 4. Oil concentration (weight %) at the producer well as a function of time (hrs) covering a period of 5 yrs (43,800 h). Red line – without pusher fluid or with constant heating well temperature boundary condition, Blue line – with pusher fluid. An example of the thermally induced stress on the deposit is shown in Figure 5. The figure shows results for the same case as Figure 3 but after only 2,016 hrs (84 days) of heating. The temperature distribution (not shown) shows that only a small portion of the deposit next to the heating well has become hot; most of the balance of the deposit remains cold. The thermal stress is negative or compressive and is as high as 2.2 MPa (deep blue color at heating well). According to Pariseau5, the compressive strength of the deposit is 10.1 MPa, so the formation will not fracture due to the thermal stresses developed during heating. 5 Pariseau, W.G., “Rock Mechanics,” H-11 Figure 5. Thermal stress (color) and pressure (contours) on the deposit, both in units of Pa. Using a constant 1000K wall temperature boundary condition at the heating well, the initial heat flux is 6 times greater than the 500 W/m2 flux applied in the previous case, resulting in faster heating of the deposit. The temperature (surface color), pressure (contours) and velocity (vectors) after a little more than 1 year of heating are shown in Figure 6. In contrast to Figure 3, Figure 6 shows a zone near the heating well encompassing approximately one-half of the deposit volume that is near a uniform temperature of 1000K. The other half of the deposit volume near the production well is at a low temperature. The flow of oil to the production well is again blocked by the kerogen that has not decomposed near the production well as shown by the velocity vector arrows in Figure 6 and by the production well concentration of oil shown in Figure 4. Figure 6. Temperature, pressure and velocity profiles at the end of a 403 d (9672 hr) heating cycle with heating by the application of a 1000K heating well temperature. H-12 DC Joule Heating of Deposit With the application of electrical current to the deposit, the nascent conductivity of the deposit should act like a resistor and heat the deposit via Joule heating. After the kerogen has decomposed, a carbonaceous char with a much higher conductivity is left behind. This higher conductivity should then move the resistive zone to where the kerogen has not yet decomposed. After 5 years of applying 1,000 volts DC to the deposit, the temperature profile, shown in Figure 7, indicates that there is essentially no heating taking place. Consequently, there is no oil or gas production from the deposit (not shown). Figure 7. Temperature, pressure and velocity profiles at the end of a 5-year heating cycle with heating by the application of 1,000 volt DC to the deposit. Combination of Conductive Heating and DC Joule Heating of Deposit The heating profile obtained by combining conductive heating (Figure 3) with DC heating (Figure 7) is shown in Figure 8. There is no substantial difference between the heating profile in Figure 8 and that shown in Figure 6, indicating that Joule heating in combination with thermal conduction heating of the deposit is not significant, even when the kerogen residue is significantly more electrically conductive than the kerogen. H-13 Figure 8. Temperature, pressure and velocity profile at the end of a 403 d (9,672 hr) heating cycle with heating by the application of a 1,000K heating well temperature and 1,000 volts DC at the heating well. RF Heating of the Deposit Results from simulations of RF heating using the quasi-static approach with a frequency of 628 MHz and an external current of 2200 A/m2 are shown in Figure 9. A very different heating profile is observed with heating in the center of the triangular zone where the electric field is highest due to the dipole antenna consisting of the two well casings. While the maximum temperature is only 665K, more of the deposit is heated to this temperature, causing the oil to move toward the production well with larger velocity vectors than with any other method of heating. Again, the area near the production well is cooler and will restrict the flow of oil, but not as severely as with conductive heating (Figures 3 and 6). From Figure 9, it is clear that oil flows to both wells so that both could be used as producer wells with RF heating. The thermal stress after 5 years of RF heating is shown in Figure 10. The maximum thermal stress is a compressive 0.97 MPa, which is well below the compressive strength of the oil shale deposit (10.1 MPa). Hence, micro-cracking is not predicted to occur in the deposit. H-14 Figure 9. Temperature, pressure and velocity profiles at the end of a 5-year heating cycle with heating by the application of 628 MHz RF with 2,200 A/m2 external current to the deposit. Figure 10. Thermal stress (color) and pressure (contours) on the deposit, both in units of Pa, at the end of a 5-year heating cycle with heating by the application of 628 MHz RF with 2,200 A/m2 external current to the deposit. Conclusions A multi-physics model of in situ extraction of oil shale has been developed which couples fluid flow, mass transfer of multiple chemical species, heat transfer, and AC (RF) and DC heating of the deposit. All physical properties used in these model equations are functions of changing local chemistry and temperature. The results show that a pusher fluid is necessary to move the oil to the production well for conductive heating of the deposit. Thermally induced stresses do not induce fracture of the deposit with any forms of heating considered in this work, but higher heat H-15 loads, especially in the case of conductive heating, will cause micro-cracking of the deposit. RF heating promotes the more uniform heating of the deposit, which is beneficial to oil extraction and lowers the thermal stresses in the deposit for a given overall heating rate. H-16