CO2, Methane, and Brine Leakage through Subsurface Pathways: Exploring Modeling, Measurement, and Policy Options Mary Kang A Dissertation Presented to the Faculty of Princeton University in Candidacy for the Degree of Doctor of Philosophy Recommended for Acceptance by the Department of Civil and Environmental Engineering Adviser: Michael A. Celia June 2014 Copyright by Mary Kang, 2014. All rights reserved. Abstract Subsurface pathways, such as abandoned oil and gas wells and faults, can serve as leakage pathways for CO2 , methane, brine, and other fluids. These pathways allow fluids from deep subsurface formations to migrate into shallow groundwater aquifers or to the atmosphere. To estimate leakage rates and the associated pressure effects on adjacent aquifers, analytical models representing fluid flow in the vicinity of leaky faults are developed in Chapter 2. The incorporation of this kind of fault model in larger basin-wide multiscale models allows sub-grid-scale effects due to leakage through faults to be captured with improved efficiency. The corresponding multi-scale framework that accounts for vertical leakage to the overlying aquifer, and horizontal flows perpendicular and parallel to a fault within a grid block, is presented in Chapter 3. In Chapter 4, first-of-a-kind direct measurements of methane fluxes from abandoned oil and gas (AOG) wells in Pennsylvania are reported. These measurements may bridge the current gap between top-down and bottom-up methane emission estimates. The mean methane flux at the 19 wells for which fluxes were measured is 0.27 kg/day/well, while the mean methane flux at the control locations near the measured wells is 4.5×10−6 kg/day/location. All measured wells showed positive methane leakage. The presence of ethane, propane, and n-butane, along with the methane isotopic composition, indicate that the methane emitted from the measured wells is predominantly of thermogenic origin. In Chapter 5, the number of AOG wells in Pennsylvania is estimated, based on historical records, and found to be in the range of 280,000 to 970,000. When the mean flux rate from the measured wells is applied to these estimated total number of wells in Pennsylvania, methane emissions are 4 to 13% of currently estimated annual statewide anthropogenic methane emissions. Three of the 19 measured wells are high emitters. Because these high emitters govern the average flux, more field iii measurements are needed. Such measurement plans should be aimed at identifying attributes that aid in finding these high emitters. Leakage was found to occur at both plugged and unplugged wells. As such, existing well abandonment regulations in Pennsylvania do not appear to be effective in controlling methane emissions from AOG wells. As a mitigation strategy, inclusion of gases emitted from AOG wells in Pennsylvania’s Alternative Energy Portfolio Standard may be valuable for both promoting capture and possible use of the gas as well as for reporting and monitoring of these wells. iv Acknowledgements There are a lot of people I have to thank for getting me to this point and would like to apologize for anyone that I missed. First, I would really like to thank my advisor, Michael Celia, who has encouraged me to develop my own research ideas and pursue them. I also like to thank him for his support and allowing me to travel to conferences around the world, all of which have shaped my views on research. I truly learned a lot from Mike and I thank him for always making time to help with everything from research questions to career advice. He has been a great mentor and has helped me enormously as I begin transitioning to a career beyond Princeton. I could not have finished this dissertation without Mike. I must thank my PEI-STEP advisor, Denise Mauzerall, for her support and positive energy. She has really taught me the difference between technical and policy research, and showed me how my work could be used to have real impacts in the world. Her stories about her own life and people she has come across in life have been very inspiring and I thank her for being a mentor and a role model. During my first few years at Princeton, I had the pleasure of working with Jan Nordbotten, who inspired much of my early work and allowed me to spend a semester in Norway. I really got a new perspective on mathematics during my time at Bergen, in addition to great cultural experiences. I must thank everyone including the faculty, students, and staff in the Math Department at the University of Bergen (especially Trine Mykkeltvedt) for making my time in Bergen easy, comfortable and memorable. I would also like to thank my dissertation committee, Jan Nordbotten, Elie BouZeid, and Peter Jaffe, and my research collaborators, Matthew Reid, Yuheng Chen, Xin Zhang, Florian Doster, Cynthia Kanno and T.C. Onstott, for their guidance, help, and support. I especially have to thank Matt and Yuheng for being available 24/7 to help me with my laboratory work. v The Princeton Environmental Institute has provided me many opportunities through a fellowship and participation in the Princeton Energy and Climate Scholars (PECS). I thank Robert Socolow, Pascale Poussart, and Holly Welles for their support through my time during and after PECS. I especially would like to thank Rob for being a wonderful and supportive mentor who was always there to give me advice and guidance on career, research, and leadership. I thank the current and past members of my research group, especially Cynthia Kanno, Benjamin Court, Florian Doster, and Ryan Edwards. Cynthia has been a wonderful co-researcher and office mate, who I have to thank for her organization skills, hard work, and positive energy. I would also like to thank the staff and faculty in the Civil and Environmental Engineering department (Jane Soohoo, Melinda Matlack, and Islam Elnaggar) for solving all administrative and computer-related problems with ease. I thank all my friends who have been there for me through good and bad. I would like to thank my life-long friends, especially Betty Kim and Hyon-Ju (Christina) Chung, and new life-long friends I met in Princeton, Kyle Keller, Katie Hansen, Carole Dalin, Brianne Smith, Darshana Narayanan, Minjin Lee, and Zhong Zheng. I also thank my FG friends, especially Lina Lee, Yuri Hsu, Boas Nahm, Nancy Sung, Lisa Yuen, and Jennie Byun, for providing me the emotional and spiritual support and making me feel at home in Princeton. Last but most importantly, I would like to thank my grandmother, my parents, my brothers, my in-laws, and my husband for all their incredible love and support. No words can express the respect and awe I have for my grandmother, someone who truly understands what hard work and sacrifice are. I can’t thank my parents, Anna and Ray Kang, enough for all their love and support in everything I’ve done, including moving far away from Toronto. I have to thank my brother, Johnny, for his patience and dependability, and my brother, Sam, for teaching me to be thankful. I thank vi my in-laws, Bun Sun and Kyung-Koo Cho, for accepting me into their family and supporting my career. And I must thank my husband, Sang, for his love, support, patience, and understanding, without which I could not have finished this dissertation. The work in Chapters 2 and 3 was supported by a National Sciences and Engineering Research Council of Canada Postgraduate Scholarship awarded to Mary Kang. This work was also supported in part by the Environmental Protection Agency under Cooperative Agreement RD-83438501 as well as the National Science Foundation under Grant EAR-0934722; the Department of Energy under Award No. DEFE0001161, CFDA No. 81,089; and the Carbon Mitigation Initiative at Princeton University. The work in Chapter 4 was supported by the Princeton Environmental Institute for the Science, Technology and Environmental Policy Fellowship and the National Sciences and Engineering Research Council of Canada for the Postgraduate Scholarship-Doctoral Program awarded to Mary Kang. This work was also supported by the Yale Center for Environmental Law and Policy, which provided funding for equipment and travel. For the work in Chapter 4, I have many people to thank in addition to my coauthors. I would especially like to thank Cheryl and Joe Thomas, who opened their home to us and were always there with a helping hand and a delicious meal. I would like to Laurie Barr and Les for taking the many hours and days in difficult terrain and weather to help us. I also thank Joann Parrick for always being welcoming and supportive of our work. I thank Save Our Streams PA for aiding in site identification and access; Joseph Vocaturo for assisting with chamber construction; Tsering W. Shawa for helping with map creation and geospatial analysis; Professor Anthony Ingraffea for providing contacts in Pennsylvania. I must thank Ryan Edwards, Evan Leister, Levi Golston, and David Pal for assisting in the field and with chamber construction; Peter Jaffe for allowing the use of laboratory equipment/facilities and David Pal for assistvii ing with laboratory equipment set-up and usage; and Levi Golston, Mark Zondlo, Kang Sun, and David Miller for providing instantaneous methane measurements. The work in Chapter 5 was supported by the Princeton Environmental Institute for the Science, Technology and Environmental Policy Fellowship. viii iX To my grandmother, Tae?Sook Lee. Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 1 Introduction 1.1 1 CO2 and methane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 Contribution to climate change and mitigation potential . . . 2 1.1.2 Fluid properties . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Subsurface leakage pathways . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Abandoned oil and gas wells . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Overview of chapters . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Analytical solutions for two-phase subsurface flow to a leaky fault considering vertical flow effects and fault properties 9 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Problem definition: leakage through faults . . . . . . . . . . . . . . . 12 2.3 General equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.1 Reduction of dimensionality . . . . . . . . . . . . . . . . . . . 18 2.3.2 Dimensionless equations . . . . . . . . . . . . . . . . . . . . . 19 2.3.3 Vertical flow effects . . . . . . . . . . . . . . . . . . . . . . . . 21 x 2.4 2.5 2.6 Solutions for leakage through faults . . . . . . . . . . . . . . . . . . . 23 2.4.1 Single-phase leakage . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.2 Two-phase leakage . . . . . . . . . . . . . . . . . . . . . . . . 27 Fault representation and leakage . . . . . . . . . . . . . . . . . . . . . 33 2.5.1 A fault model . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.5.2 Effect of fault properties . . . . . . . . . . . . . . . . . . . . . 35 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3 A multi-scale model for leakage through faults using combined analytical-numerical methods 40 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 Coarse-scale numerical model . . . . . . . . . . . . . . . . . . . . . . 41 3.3 Fine-scale model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3.1 Horizontal flow perpendicular to fault . . . . . . . . . . . . . . 44 3.3.2 Vertical leakage through fault . . . . . . . . . . . . . . . . . . 46 3.3.3 Pressure correction . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.4 Horizontal flow parallel to fault . . . . . . . . . . . . . . . . . 49 Multi-scale model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.4.1 Numerical example . . . . . . . . . . . . . . . . . . . . . . . . 52 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4 3.5 4 Direct measurements of methane emissions from abandoned oil and gas wells in Pennsylvania 55 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.2.1 Site selection . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.2.2 Flux chambers and sampling . . . . . . . . . . . . . . . . . . . 59 4.2.3 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 xi 4.2.4 Flux calculation . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.3.1 Methane fluxes . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.3.2 Presence of ethane, propane, and n-butane . . . . . . . . . . . 69 4.3.3 Carbon isotopes of methane . . . . . . . . . . . . . . . . . . . 69 4.4 Uncertainties in methane fluxes . . . . . . . . . . . . . . . . . . . . . 72 4.5 Methane emissions from abandoned wells in PA . . . . . . . . . . . . 75 4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.3 5 Estimating and mitigating methane emissions from abandoned oil and gas wells - Policy analysis and implications 79 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.2 Estimating methane emissions from AOG wells . . . . . . . . . . . . 82 5.2.1 Number of AOG wells . . . . . . . . . . . . . . . . . . . . . . 82 5.2.2 Emission estimates per AOG well . . . . . . . . . . . . . . . . 89 5.2.3 Reporting and monitoring . . . . . . . . . . . . . . . . . . . . 91 Mitigation options . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3.1 Plugging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.3.2 Usage as an alternative energy source? . . . . . . . . . . . . . 96 5.3 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6 Concluding Remarks 6.1 105 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.1.1 Model development . . . . . . . . . . . . . . . . . . . . . . . . 107 6.1.2 Measurements and policies . . . . . . . . . . . . . . . . . . . . 107 Bibliography 109 xii List of Tables 1.1 Properties of CO2 , methane, and water/brine . . . . . . . . . . . . . 4 2.1 Aquifer and fluid properties . . . . . . . . . . . . . . . . . . . . . . . 13 4.1 Summary of sampling campaigns . . . . . . . . . . . . . . . . . . . . 59 4.2 Estimates of the number of abandoned oil and gas wells in Pennsylvania. 5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The number of abandoned and orphaned oil and gas wells on Pennsylvania DEP’s database 5.2 76 . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Estimate of the number of AOG wells in Pennsylvania based on changes in technology and regulatory environment . . . . . . . . . . . . . . . 86 5.3 Well attributes and their impact on well bore leakage . . . . . . . . . 90 5.4 Measured methane fluxes by attributes . . . . . . . . . . . . . . . . . 91 5.5 Summary of plugging requirements in Pennsylvania Code . . . . . . . 102 5.6 Summary of alternative energy resource types included in Pennsylvania’s Alternative Energy Portfolio Standard (AEPS) . . . . . . . . . . 103 5.7 Alternative energy sources for Pennsylvania in PJM-EIS Generation Attribute Tracking System (GATS) . . . . . . . . . . . . . . . . . . . 104 xiii List of Figures 1.1 Schematic of CO2 injection, migration, and interaction with faults and existing oil and gas wells . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1 A schematic of a vertical cross-section of an aquifer containing a fault 13 2.2 ˆ f ault = −0.05 Analytical solutions for single-phase leakage given Q m2 /day with (a) RL = RL2 = 1.0, and k = 26.7, and (b) RL = 0.25, RL2 = 0.0625, and k = 426 . . . . . . . . . . . . . . . . . . . . . . . . 2.3 27 Analytical solutions for two-phase leakage with (a) Qd = 0.23 and (b) ˆ f ault = −0.2 m2 /day, RL = 5, R2 = 25, and k = 1.0 31 Qd = 0.24, given Q L 2.4 Analytical solutions for two-phase leakage with (a) Qd = 0.22 and (b) ˆ f ault = −0.2 m2 /day, RL = 1.14, R2 = 1.31, and k Qd = 0.24, given Q L = 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Comparison of Qd values given by the numerical simulation, representing various T and m values, and the fault model (equation (2.42)) . . 2.6 35 Numerical and analytical solutions for single-phase leakage with (a) ˆ f ault = −0.05 m2 /day and m = 1 . T = 0.04 and (b) T = 0.2, given Q 2.7 31 36 Numerical and analytical solutions for two-phase leakage with (a) T = ˆ f ault = 0.2 (and Qd = 0.209) and (b) T = 0.8 (and Qd = 0.217), given Q −0.2 m2 /day and m = 1 . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 37 Numerical and analytical solutions for single-phase leakage with (a) ˆ f ault = −0.05 m2 /day and T = 1.0 . . m = 1 and (b) m = 10, given Q xiv 37 2.9 Numerical and analytical solutions for two-phase leakage with (a) m = ˆ f ault = 1 (and Qd = 0.218), and (b) m = 10 (and Qd = 0.219), given Q −0.2 m2 /day and T = 1.0 . . . . . . . . . . . . . . . . . . . . . . . . 3.1 38 The empirical relationship between Qd , Qf ault and h0 determined using the fault model and numerical simulations . . . . . . . . . . . . . . . 45 3.2 A schematic of grid blocks containing faults and adjacent grid blocks 50 3.3 Pressures (left) and saturations (right) at 7 years of the numerical example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 The leakage rate, Qf ault , (left) and the proportion of the leakage that is brine, Qd , (right) with respect to time for the numerical example . 4.1 57 Schematic of flux chamber collar and frame enclosing surface protrusions of an abandoned oil and gas well . . . . . . . . . . . . . . . . . 4.3 53 Map of abandoned oil and gas wells measured and/or on Pennsylvania Department of Environmental Protection list . . . . . . . . . . . . . . 4.2 53 60 The R2 value for the linear regression and the corresponding methane fluxes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.4 Methane fluxes measured at well and control locations . . . . . . . . 68 4.5 Average alkane ratios (A) and proportion of samples with alkane ratios greater than 0.01 at control and well location with detectable ethane, propane, and n-butane concentrations (B) . . . . . . . . . . . . . . . 4.6 70 Composition of carbon isotope of methane for select samples collected at well and control locations . . . . . . . . . . . . . . . . . . . . . . . 71 5.1 Distribution of AOG wells on DEP’s list by county in Pennsylvania . 84 5.2 Historic data on number of wells drilled in Pennsylvania . . . . . . . 85 xv Chapter 1 Introduction Subsurface pathways, such as abandoned oil and gas wells and faults, can serve as leakage pathways for CO2 , methane, brine, and other fluids. These pathways allow fluids from deep subsurface formations to migrate into shallow groundwater aquifers or to the atmosphere. Abandoned oil and gas wells and conductive faults in geologic formations serve as leakage pathways for fluids which would otherwise remain trapped beneath low-permeability layers. Fluids of concern in this dissertation are the greenhouse gases (GHGs), CO2 and methane, and brine which is present in deep subsurface formations. Leaky faults and abandoned wells must be considered in many engineering activities such as geologic sequestration of CO2 [51, 46, 68], deep well injection of hazardous and radioactive wastes [82, 83, 10], underground natural gas storage [27, 14], and oil and gas exploration and production [42, 29]. To understand the leakage potential and to reduce leakage, model development, field measurements, and governmental and regulatory policies are needed. 1 1.1 1.1.1 CO2 and methane Contribution to climate change and mitigation potential CO2 and methane are well-mixed greenhouse gases (GHGs) that greatly contribute to positive radiative forcing, which leads to global surface warming and climate change [35]. Radiative forcing represents the change in energy fluxes caused by changes in anthropogenic and natural drivers and is one of the metrics used to determine the relative contribution of an emitted compound to surface warming. The contributions to radiative forcing from anthropogenic emissions of CO2 and methane are 1.68 W m−2 and 0.97 W m−2 respectively for 2011 relative to 1750 [35]. CO2 is the dominant greenhouse gas and mitigation strategies aimed at CO2 emission reductions are needed to limit climate change and its impacts [66, 34, 35]. Methane is a strong GHG with a global warming potential (GWP) 86 times greater than CO2 in a 20-year time frame [35]. As a result, mitigating methane emissions can significantly limit climate change in the near term [79], while we wait for CO2 emissions reduction strategies to be implemented. Available technologies for reducing CO2 emissions can be broadly categorized as CO2 capture and storage (CCS), energy efficiency and conservation, renewable electricity and fuels (e.g. solar and wind), nuclear fission, and forests and agricultural soils [66]. To stabilize atmospheric concentrations of CO2 at current levels of 399 ppmv (as of March 2014), we need to deploy mitigation strategies that can reduce emissions of carbon by more than 9 billion tons per year. This requires the deployment of multiple available mitigation strategies, which include CCS. The International Energy Agency (IEA)’s technology roadmap specifies a goal of at least 30 successful CCS demonstration projects by 2020 [36]. There are 4 projects operational as of 2013 [36]. Additional policies and incentives are needed to reach IEA’s goal for CCS. 2 Given the many co-benefits of mitigating methane emissions, including air quality and agricultural productivity and the value of methane as an energy source, there are many low-cost opportunities to reduce methane emissions [79]. In fact, five out of the eight low-cost measures for reducing short-lived climate forcers (e.g. methane, tropospheric ozone and black carbon) are related to methane [84]. Three of these measures are related to the oil and gas industry: recovery and utilization of vented gas in oil production, reduced leakage during gas pipeline transmission, and recovery and utilization of vented gas in oil production. CO2 and methane emissions are driven mainly by anthropogenic activities, primarily from the fossil fuel industry [35]. World energy consumption is expected to continue to grow, reaching over 800 quadrillion Btu by 2040 [85], with much of the demand being met by fossil fuels. With current regulations and policies, CO2 emissions will grow by 45% by 2040, relative to 2010 [85]. Meanwhile, the need and urgency to reduce GHG emissions are becoming increasingly important. With the continued growth and development of the oil and gas industry, it is important to look to this industry for additional low-cost mitigation strategies. 1.1.2 Fluid properties The fluids of primary concern in this dissertation are carbon dioxide, methane, and brine. The relevant properties of these fluids are listed in Table 1.1. The properties are presented in terms of pressure and temperature conditions found at the ground surface and in “shallow formations” (1 km deep) following ref. [57]. Shallow formations are considered in the context of two different kinds of basins: cold and warm basins. A cold basin is defined to have a surface temperature of 10◦ C and a geothermal gradient of 25◦ C/km. A warm basin is defined to have a surface temperature of 20◦ C and a geothermal gradient of 45◦ C/km. The pressure gradient is assumed to be 10.5 MPa/km. At the ground surface, water is assumed to be fresh; in formations 1 km 3 below ground, the resident water is likely to be saline. Table 1.1 shows that in all three cases, CO2 and methane are significantly lighter and less viscous than water and brine, with methane being lighter and less viscous than CO2 . Table 1.1: Properties of CO2 , methane, and water/brine at the ground surface and in “shallow formations” (1 km deep). Two types of shallow formations are considered: cold and warm basins. Data are from ref. [57] and www.peacesoftware.de. Densities (kg/m3) Viscosities (mPa s) Compound Methane CO2 Water Brine* Methane CO2 Water Brine* Ground Surface 0.66 1.8 1000 0.010 0.015 1.002 Cold Basin 76 714 Warm Basin 66 307 1012 - 1230 0.014 0.0577 998 - 1210 0.015 0.023 0.795 - 1.58 0.491 - 0.883 *Salinities from 0.02 to 0.3, measured in mole fractions Given their fluid properties, methane and CO2 will migrate upward from geologic formations saturated with water or brine due to buoyancy (Table 1.1). In contrast, brine is denser than freshwater and will require significant pressure gradients to promote its leakage into shallower formations. Nonetheless, leakage rates of methane and CO2 are influenced by the presence and flow of brine and its properties. 1.2 Subsurface leakage pathways Abandoned oil and gas wells and conductive faults can extend over multiple sedimentary layers connecting deeper formations to groundwater aquifers and to the surface (Figure 1.1). Sedimentary basins are typically composed of lower permeability layers such as shales and higher permeability layers such as sandstone [70, 57]. Oil and gas reservoirs are the results of oil and gas migrating upwards from the source rock and being collected under low-permeability layers, also known as caprocks (hatched layers in Figure 1.1). Oil and gas wells are drilled through these caprocks and if improperly 4 abandoned, can compromise the ability of the caprocks to keep fluids from migrating upwards. Faults are the results of movements in the earth, which change the properties of the rock in the affected zone. The resulting properties of faults are variable and faults can act as a barrier, a conduit, or a barrier-conduit system for subsurface flow. As shown in Figure 1.1, faults can extend through multiple layers and can also compromise the ability of caprocks to keep fluids from migrating upwards. Depending on the permeability and extent/geometry of abandoned oil and gas wells and faults, these wells and faults can be a significant pathway for vertical migration of fluids. This is especially true for methane and CO2 , which are generally much less dense than the resident fluids (i.e. brine/water). lt Fau Upward Flow Figure 1.1: Schematic of CO2 injection, migration, and interaction with faults and existing (abandoned and producing) oil and gas wells (adapted from ref. [25]). 5 Within the context of CCS, specifically geologic storage of CO2 , the presence of abandoned wells and faults poses a risk of leakage. To understand this risk, basinscale models of CO2 and brine flow that cover areas in the order of 105 km2 are needed [90]. In these basin-scale models, abandoned wells and faults are relatively small-scale features and their representation requires large computational resources or an alternate approach to represent the local-scale effects. An alternate approach is to use a multi-scale model in which the local-scale effects are captured using an analytical model within a coarser-scale numerical model. This strategy is taken in ref. [24] to represent leaky abandoned wells. A similar strategy can be employed for faults but requires the additional consideration of along-fault flows. Accurate representation of these local-scale effects are important since they govern the leakage rate, which is a function of the properties of the aquifer, the leakage pathway, and the fluids and the boundary conditions driven by the coarse-scale model. 1.3 Abandoned oil and gas wells The number of oil and gas wells in the U.S. and in other oil and gas producing regions is highly uncertain. For the U.S., ref. [11] estimates that there are 3 million abandoned wells. In Pennsylvania, where the first commercial oil well was drilled in 1859, estimates range from 300,000 to 500,000 but may be up to 970,000 (as discussed in Chapter 5). In Alberta, there is an estimated 315,000 abandoned wells [87]. Similar numbers of wells are expected in states with long histories of oil and gas production such as California, Colorado, and Texas. Although abandoned oil and gas wells have been studied as a leakage risk for geologic storage of CO2 [58, 56, 57], there has been limited research on the potential of these wells to emit methane. Oil and gas wells are abandoned when they no longer produce economical amounts of oil and/or gas. Oil and gas, including methane 6 and other gases, continue to exist in the formations, even after the wells have been abandoned. Depending on how the wells were abandoned, these fluids, especially buoyant gases, will migrate from these formations to overlying groundwater aquifers and to the surface. Previous to the work presented in Chapter 4, there were no empirical studies on the amount of leakage along these wells. Abandonment procedures for oil and gas wells have been motivated mainly by oil and gas resource conservation and protection, and groundwater protection. Therefore, the main abandonment strategy is plugging, which typically involves zonal isolation, whereby formations from which leakage can occur (i.e. oil and gas reservoirs) are isolated from formations that should be protected (i.e. groundwater aquifers). While regulations exist for abandonment procedures and protocols, there is no regulation to address methane emissions from abandoned oil and gas wells. This is not surprising since methane emissions from these wells are not included in any emissions inventories [11], and the implied assumption in abandonment regulations is that leakage will not occur. The result is a lack of information to quantify methane emissions from these wells, which depend on the number of wells and the methane emissions per well. Both of these are uncertain, which makes evaluating the effectiveness of any emissions reduction strategies difficult. 1.4 Overview of chapters The modeling, measurement, and policy options for subsurface leakage pathways, specifically abandoned oil and gas wells and faults, are presented in this dissertation. These models, empirical field studies, and policies provide the necessary scientific information to determine the significance of these leakage pathways and to develop and implement effective carbon mitigation strategies and policies. 7 Chapters 2, 4, and 5 have been published or will be submitted for publication in a peer-reviewed journal. Chapter 3 will be published once additional model verification and validation are performed. In Chapter 2, a model for single- and two-phase leakage through faults considering vertical flow effects and fault properties is presented in the context of geologic storage of CO2 . This work was done in conjunction with Jan M. Nordbotten, Florian Doster, and Michael A. Celia, who are the co-authors of the corresponding paper published in Water Resources Research. In Chapter 3, a multi-scale framework is presented. This multi-scale model uses solutions derived in Chapter 2 as fine-scale representations for single- and two-phase leakage through faults. This work was done in conjunction with Jan M. Nordbotten and Michael A. Celia, and was presented at the Society of Industrial and Applied Mathematics Geoscience Conference, Padova, Italy, June 17-20, 2013. In Chapter 4, first direct field measurements of methane leakage from abandoned oil and gas wells in Pennsylvania are presented, along with measurements of other hydrocarbons and carbon isotopes of methane to gain insight into the source of the emitted methane. This work was done with 4 different research groups at Princeton University and the following collaborators are co-authors of the corresponding paper: Cynthia Kanno, Matthew C. Reid, Xin Zhang, Denise L Mauzerall, Michael A. Celia, Yuheng Chen, and Tullis C. Onstott. Some of the material in this chapter has been presented at the American Geophysical Union, Fall Meeting, San Francisco, CA, December 9-13, 2013. In Chapter 5, policies for abandoned oil and gas wells, specifically relating to improving methane emission estimates and developing and implementing mitigation strategies, are presented in the context of Pennsylvania. The co-authors of this work are Denise L. Mauzerall and Michael A. Celia. Concluding remarks are presented in Chapter 6. 8 Chapter 2 Analytical solutions for two-phase subsurface flow to a leaky fault considering vertical flow effects and fault properties 2.1 Introduction Conductive faults in geologic formations serve as leakage pathways for fluids which would otherwise remain trapped beneath low-permeability layers. To estimate leakage and the associated pressure effects on adjacent aquifers, flow in and around faults must be modeled. Computational efficiency of the modeling tool is important, especially if we consider leaky faults in the context of basin-scale representations [24, 90] or the need to perform risk analysis [12, 55]. One approach to achieve computational efficiency in basin-scale models is to employ analytical models locally, in the vicinity of the fault, within a multi-scale model. The multi-scale model in consideration here consists 9 of a two-dimensional (2-D) vertically-integrated numerical model at the coarse-scale and various analytical or empirical representations at the fine-scale. Steady-state analytical methods have been successfully employed as sub-grid-scale representations of producing oil and gas wells [67] and CO2 and brine leakage through abandoned wells [24] in reservoir or basin-scale numerical simulations. Similar to wells, faults are relatively small-scale features in basin-scale numerical models. However, in a 2-D vertically-integrated space, faults are linear features that can extend over multiple grid blocks in the horizontal direction. This is in contrast to abandoned wells that can be viewed as point features in a 2-D domain. Therefore, a multi-scale model in basins containing faults requires methods to represent sub-grid-scale (local-scale) processes and to handle inter-grid-block fluid transport. Within a numerical domain of rectangular grid blocks, flows associated with faults, both within a grid block and between grid blocks, can be viewed as a combination of flows perpendicular and parallel to the fault. Accurate and efficient representation of the local scale flow within a grid block, perpendicular to the fault, is the primary motivation for the analytical solutions and fault representations developed here. Analytical models designed for leakage through fractures, faults, and other planar features have been developed assuming vertical equilibrium [71, 86]. Vertical equilibrium conditions exist when vertical dynamics occur quickly such that vertical flows are equilibrated and horizontal flows dominate. The assumption of vertical equilibrium, also referred to as the Dupuit assumption, is commonly used to reduce the dimensions of problems describing flow through porous media [9, 33, 57]. Applicability of the vertical equilibrium assumption has been explored based on numerical simulations [15, 16] as well as analysis of the governing equations [59]. Perturbation analysis has also been used to develop criteria to determine when the vertical equilibrium assumption is valid [88, 44, 19]. According to these criteria, vertical equilibrium solutions are 10 valid for isotropic or mildly anisotropic aquifers with characteristic vertical lengths that are small relative to characteristic horizontal lengths. Aquifers in sedimentary basins tend to be anisotropic with vertical permeabilities that are smaller than horizontal permeabilities. Scenarios where vertical equilibrium is not achieved are possible, especially adjacent to a leakage pathway, where vertical flows associated with leakage can drive significant vertical flows in the vicinity of the leakage feature. To address this, ref. [56] improved the analytical solution for upconing around a well by imposing a linear structure to the vertical flow field, thereby relaxing the vertical equilibrium assumption. This solution is not directly applicable to faults due to the difference in geometry (i.e. radial versus Cartesian). Furthermore, only aquifer properties were explicitly considered and variations due to different wellbore properties were not considered. Fault properties, such as fault permeability, fault width, and anisotropy in the fault permeability, have been shown to affect groundwater (single-phase) flow field in the adjacent aquifer [28, 78, 4, 5, 89]. In petroleum reservoir modeling, faults are commonly defined either implicitly by modifying properties of the surrounding grid blocks, or explicitly as a volume of porous medium. The conventional approach in production simulation models is to implicitly represent the effect of faults on fluid flow using transmissibility multipliers based on empirical estimates of fault permeability and thickness [50, 49, 23]. In this paper, fault properties and vertical flow effects are incorporated into a vertically-integrated two-phase flow model and steady-state analytical solutions are derived for single- and two-phase flow adjacent to a conductive fault. Numerical simulations are performed to establish a relationship between fault properties and flow to a fault in an aquifer. The steady-state solutions and the fault model are mainly designed for future integration into a basin-wide multi-scale model, in which the coarse scale is a 2-D vertically-integrated numerical model. Note that these 211 D vertically-integrated numerical models can be stacked to represent the vertical sequence of aquifer-aquitard systems [60]. In section 2.2, the problem of leakage through a fault and the associated flow in a single confined aquifer is defined. In section 2.3, general equations for two-phase flow and the representation of vertical flow effects are presented. The time derivative is included in these general equations for completeness. In section 2.4, steady-state analytical solutions are derived for single- and two-phase leakage through a fault. In section 2.5, a fault model that relates fault properties to inputs of the analytical solutions is proposed and the effects of fault properties are explored through comparisons of the analytical and numerical solutions. In sections 2.4 and 2.5, we study the role of vertical flow effects and fault properties by comparing plume shapes. Finally, concluding remarks are provided in section 2.6. 2.2 Problem definition: leakage through faults Although leakage through faults is a general problem with applications in many different engineering disciplines, we develop our problem here in the context of geologic storage of CO2 in saline aquifers. Therefore, the two fluids of concern are CO2 and brine. Values for parameters in illustrative examples throughout the paper are based on ref. [8] and are summarized in Table 2.1. CO2 and brine properties are calculated from pressures, temperatures, and salinities at conditions typical of CO2 injection in North America. We consider a relatively deep (3000 m) aquifer, with a surface temperature of 10◦ C, a geothermal gradient of 25◦ C per kilometer, and a salinity of about 0.2 mole fraction. A single conductive fault zone in a horizontal confined homogeneous and anisotropic aquifer is conceptualized as shown in Figure 2.1. Leaky faults can be conceptualized as consisting of a low-conductivity (i.e. low-permeability) core that 12 Table 2.1: Aquifer and fluid properties. Parameter CO2 density, ρˆl Brine density, ρˆd CO2 viscosity, µ ˆl Brine viscosity, µ ˆd ˆ Aquifer thickness, H ˆ0 CO2 thickness at the outer boundary, h Value 733 kg/m3 1099 kg/m3 0.0611 mPa s 0.511 mPa s 5m 1.25 m ˆ f ault Q No Flow No Flow ˆ inner h ˆf W ˆ x) h(ˆ x ˆ CO2 ˆ0 h zˆ ˆ H Brine No Flow ˆa = x L ˆouter − x ˆinner Figure 2.1: A schematic of a vertical cross-section of an aquifer (white area) containing a fault (gray area). A gray line is used to represent the fluid-fluid interface given leakage through the fault. Note that leakage leaving the aquifer through the fault is defined to be negative. prevents flow horizontally across the fault, surrounded by two high-conductivity (i.e. high-permeability) damage zones that permit flow vertically along the fault [21]. This highly conductive damage zone poses a risk of leakage when it extends into overlying aquifers. This damage zone is referred to here as “the fault zone” or simply “the fault” and is depicted as the gray area in Figure 2.1. The low conductivity core is represented as a no flow boundary. We seek to represent one fault zone and the adjacent aquifer in this paper. Within a multi-scale framework, a different grid block can be used to represent the fault zone and the corresponding aquifer region lying on the opposite side of the fault core. In addition, the no-flow assumption for the fault core can be relaxed in a multi-scale model. Given typical aquifer thicknesses that are small relative to the vertical extent of the fault, the fault is viewed as fully-penetrating and vertical within the aquifer. The 13 ˆ f , and finite permeabilities fully-penetrating vertical fault zone has a finite width, W parallel, kˆzf , and perpendicular, kˆxf , to the fault. (We assume here the kˆyf is in the same order of magnitude as kˆxf and kˆzf .) (Note that dimensional variables are denoted by a hat.) Permeabilities of conductive faults can be several orders of magnitude larger than the aquifer [78]. In this paper, we focus on the xˆ- and zˆdirections and conditions when min(kˆxf , kˆzf ) ≥ max(kˆx , kˆz ), where kˆx and kˆz are the horizontal and vertical permeabilities of the aquifer. We assume that changes in fault or aquifer properties and the impact of boundary conditions in the yˆ-direction (into the page in Figure 2.1) are small for the domain in consideration. The presence of a conductive fault zone naturally orients the flow in the aquifer to be perpendicular to the fault, following the Tangent Law [57], and thus, the predominant horizontal flow direction in the aquifer is assumed to be in the xˆ-direction. The permeabilities in the fault are anisotropic with kˆzf equaling or exceeding kˆxf , and the fault anisotropy ratio, m = kˆzf /kˆxf , is equal to or greater than one. This is in contrast to aquifer permeabilities, which are typically larger in the horizontal direction than in the vertical direction. The presence of a low conductivity core and large and anisotropic fault permeabilities promotes predominantly vertical flow in the fault. In the confined aquifer adjacent to the fault, the two fluids are assumed to be separated by a sharp interface, as shown in Figure 2.1. This corresponds to strong buoyancy and a negligible capillary transition zone. For many problems related to geologic sequestration of CO2 , fluid segregation due to buoyancy will occur in the order of years to tens of years [16]. Given a 50-year injection period, which is typical for geologic sequestration, density segregation can generally be assumed in post-injection leakage scenarios but also during injection in many cases. The sharp interface sepaˆ x, yˆ, tˆ), where tˆ is time. The interface rating the fluid regions is defined to be at zˆ = h(ˆ shape and location will evolve over the post-injection period. Local steady-state so14 lutions, adjacent to the fault, as derived in this paper, can be used as derived or can be used to capture transient behavior using the method of successive steady-states [9]. Within a multi-scale model, a new “steady-state” solution would be employed at every time step of the coarse-scale numerical model, driven by transient changes in the coarse-scale model. ˆ and extends over a finite The confined aquifer has a uniform finite thickness, H, ˆ a . Typical thicknesses of formations in which CO2 injection is horizontal length, L being considered, such as the Mt. Simon formation [90], can vary from a few meters to hundreds of meters. In this paper, the length of the aquifer is chosen to be 5.2 times greater than the aquifer thickness of 5 m. Using this same ratio for a 100 m thick aquifer, the horizontal extent would be 520 m, which is in the range of typical lengths for grid block sides in basin-scale models. Within a multi-scale model, the ˆ a , which should be sufficiently dimensions of the grid block can be used to define L large to capture the flow field in the region significantly affected by a leaky fault and more importantly, the region where vertical flows are significant. The domain in the yˆ-direction is assumed to be sufficiently small and CO2 injection is assumed to occur at a sufficient distance away from the fault such that variations in the yˆdirection caused by the radial nature of the injection are negligible. Therefore, the fault and aquifer permeabilities in the yˆ-direction are not explicitly considered in the local solution (although they would be included in the coarse-scale model). The main objective here is to represent local-scale horizontal flow perpendicular to the fault (ˆ x-direction) in the aquifer for inclusion in a multi-scale model. The idea is to have variations along the fault in the yˆ-direction captured in the multi-scale model as a combination of sub-grid-scale and inter-grid-scale representations [20]. With ˆ ∂h ∂ yˆ = 0 in the local-scale model, vertical integration of horizontal flows throughout the thickness of the aquifer yields a flow system described by one-dimensional partial ˆ = h(ˆ ˆ x, tˆ). differential equations, which can be solved for h 15 The location where the fault meets the aquifer is at xˆ = xˆinner . At this location, ˆ inner , is expected to be at a minimum given leakage of CO2 or CO2 the CO2 thickness, h ˆ inner is allowed to vary and brine. For both single- and two-phase leakage conditions, h from 0 to the thickness of the CO2 plume at the outer boundary. This is in contrast ˆ inner is assumed to be zero for all two-phase leakage conditions. to ref. [56] where h ˆ inner is found In this work, we found this assumption to be invalid, and instead, h to be a function of aquifer and fault properties. At the outer boundary located at xˆ = xˆouter , CO2 occupies the top portion of the formation corresponding to a thickness ˆ 0 = h(ˆ ˆ xouter ), while the remaining H ˆ 0 is saturated with brine. Within a multiˆ −h of h scale model, this outer boundary would be defined by the fluid distribution in the adjacent grid block. Vertical equilibrium would be assumed in this adjacent grid block since the coarse-scale model is a 2-D vertically-integrated model. Furthermore, a sharp interface between the two fluids is assumed to be established in the adjacent grid block. Leakage of one or both fluid phases is defined as any fluid that arrives at xˆinner and enters the fault. Since the local-scale domain does not extend above the top of the aquifer, fluid(s) leaving the top of the fault are assumed to continue flowing upward in the fault model. The influence of overlying aquifers can be captured if the aquifer-fault system modeled here exists in a multi-layer system of aquifers and aquitards [60]. Specifically, the pressure difference in the fault between the bottom ˆ f ault of the overlying aquifer and the top of the source aquifer can be used to obtain Q [L2 T−1 ], defined as the total leakage out of the fault zone at the top of the fault per ˆ f ault is comprised of one or two fluid(s) unit horizontal length in the yˆ-direction [24]. Q and in the case of two-phase leakage, it is important to determine the proportion of ˆ f ault represented by each fluid. At steady-state conditions, there is no accumulation Q 16 in the fault zone and the proportion of leakage for each fluid, α, is given as Qα = ˆ H 0 qˆα dˆ z ˆ f ault Q = ˆf W f ault qˆα,z top ˆ f ault Q dˆ x = ˆα Q ˆα Q f ault where qˆα [LT−1 ] is the horizontal flux of fluid α in the aquifer, qˆα,z (2.1) top [LT−1 ] is the ˆ α [L2 T−1 ] is the vertical flux of fluid α in the fault at the top of the fault (ˆ z = 0), Q total leakage of fluid α, and α is equal to d for the denser fluid and l for the lighter fluid. The horizontal flux, qˆα , describes flows perpendicular to the fault and parallel to the horizontal top and bottom of the aquifer. Positive qˆα values indicate flux from the fault to the right and into the aquifer (i.e. injection) and negative qˆα values ˆ f ault correspond to flux into the fault from the aquifer (i.e. leakage/extraction). Q ˆ f ault , leaving follows the same sign convention as qˆα . For example, the leakage, Q the system as shown in Figure 2.1 is a negative number. In this paper, we consider ˆ0 < H ˆ with (1) single-phase leakage of the lighter fluid (CO2 ) situations when 0 < h only (Qd = 0), and (2) two-phase leakage of both the lighter (CO2 ) and denser (brine) fluids (0 < Qd < 1). 2.3 General equations Although steady-state solutions will be based on one-dimensional representations, we begin with the full three-dimensional equations to review the assumptions made. Two-phase flow through porous media assuming negligible capillary pressure can be represented using the multiphase extension of Darcy’s Law, defined as ˆαk ˆ · (∇ˆ q ˆα (ˆ x, yˆ, zˆ, tˆ) = −λ p + ρˆα gˆ∇ˆ z ), (2.2) ˆ [L2 ] is the intrinsic permeability where q ˆα [ML−3 ] is the Darcy flux for each fluid, k ˆ α = krel,α /ˆ tensor, λ µα [M−1 LT] is the mobility tensor of each fluid, krel,α [-] is the 17 relative permeability tensor of each fluid, µ ˆα [ML−1 T−1 ] is the viscosity of each fluid, pˆ [ML−1 T−2 ] is the pressure of both fluids, ρˆα [ML−3 ] is the density of each fluid, and gˆ [LT−2 ] is the gravity acceleration. Mass conservation of each fluid in the porous medium is given by ∂(φˆ ρα sα ) + ∇ · (ˆ ρα q ˆα ) = ρˆα Fˆα , ∂ tˆ (2.3) where φ [-] is the porosity, sα [-] is the saturation, and Fˆα [T −1 ] is the source/sink term. In section 2.3.1, the spatial dimensions of equations (2.2) and (2.3) are reduced using vertical integration and neglecting variations in the y-direction (see discussion in section 2.2). Next, characteristic scales are defined and the governing equations are expressed in dimensionless form in section 2.3.2. In section 2.3.3, vertical flow effects and a structured representation of vertical flow, introduced in ref. [56], is reviewed and included in the pressure and mass balance equations. The result is a set of governing equations that include vertical flow effects. The equations are especially valuable in leakage scenarios but are applicable to any situation where vertical flow effects are significant and a linearly structured representation of the vertical flow field is suitable. 2.3.1 Reduction of dimensionality To reduce the spatial dimensions of equations (2.2) and (2.3), vertical equilibrium and a sharp interface can be assumed and vertical integration can be applied [57]. The resulting 2-D partial differential equations (PDEs) can be solved for the interface ˆ = h(ˆ ˆ x, yˆ, tˆ). We neglect variations in the y-direction (i.e. location, h ˆ ∂h ∂ yˆ = 0), ˆ = h(ˆ ˆ x, tˆ). For for reasons stated in section 2.2, to obtain a set of 1-D PDEs for h ˆ the integrated equations horizontal homogeneous confined aquifers with thickness, H, 18 have the following form, ˆ ∂h ∂ φ + ˆ ∂ xˆ ∂t ˆ ∂h ∂ −φ + ∂ tˆ ∂ xˆ ˆ h qˆl dˆ z = 0, (2.4) 0 ˆ H ˆ h qˆd dˆ z =0 (2.5) for the lighter, l, and denser, d, fluids, respectively. In equations (2.4) and (2.5), incompressible fluids and solids, no residual saturation of either fluids, and no source/sink terms are assumed. The fluid viscosities and densities and the aquifer properties are assumed to be uniform in space and constant in time. We assume that the principal axes of the permeability tensor are in line with the coordinate axes. Thus, the horizontal flux, qˆα , away from the fault in the x-direction is ˆ α kˆx qˆα (ˆ x, zˆ, tˆ) = −λ ∂ pˆ ∂ xˆ , (2.6) α ˆ α = krel,α /ˆ where λ µα [M−1 LT−1 ], kˆx is the horizontal permeability [L2 ], and ∂ pˆ ∂x ˆ α is the horizontal pressure gradient in each fluid region. We assume here that the xˆ- and zˆ-axes are aligned with the principle directions. Given a sharp interface, krel,α = 1 since we assume that fluids have been segregated by gravity into two zones that are fully saturated with one fluid [9]. 2.3.2 Dimensionless equations To rewrite equations (2.4), (2.5), and (2.6) in dimensionless form, the independent variables are scaled as follows x= xˆ , ˆa L z= zˆ , ˆ H t= 19 tˆ ˆL ˆ a /Q ˆ f ault ) (φH . (2.7) The horizontal characteristic length scale is determined by the extent of the aquifer, ˆ a . The vertical characteristic length scale is chosen to be H. ˆ The dimensionless L dependent variables are defined as p(x, z, t) = pˆ(ˆ x, zˆ, tˆ) ˆ (ˆ ρd − ρˆl )ˆ gH , h(x, z, t) = ˆ x, zˆ, tˆ) h(ˆ , ˆ H q(x, z, t) = ˆ qˆα H . ˆ f ault Q (2.8) For pressure, the characteristic scale is defined based on buoyancy effects. The characˆ The remaining teristic scale for the thickness of the lighter fluid, h, is chosen to be H. dimensionless groups are λ= ˆl λ , ˆd λ ρα = ρˆα , (ˆ ρd − ρˆl ) Γ=δ ˆ d (ˆ ˆ kˆx λ ρd − ρˆl )ˆ gH , ˆ f ault Q (2.9) ˆ L ˆ a , λ is the mobility ratio, ρα is the density of fluid α normalized by the where δ = H/ density difference between the two fluids, and Γ is the ratio of buoyancy to viscous forces scaled by δ. Note that Γ is defined here with respect to the denser fluid. Based on equations (2.7), (2.8), and (2.9), the dimensionless forms of equations (2.4) and (2.5) are ∂h ∂ = − ∂t ∂x ∂h ∂ = ∂t ∂x h ql dz , (2.10) 0 1 qd dz , (2.11) h where ∂p , ∂x d ∂p ql = −Γλ . ∂x l qd = −Γ 20 (2.12) (2.13) 2.3.3 Vertical flow effects The vertical equilibrium assumption has been shown to be reasonable when RL > √ −1 10 [44] and/or when RL2 1 [88], where RL = δ k and k = kˆx /kˆz . This implies that vertical flow effects can be important if the horizontal length scale is not sufficiently larger than the vertical length scale and/or if kˆx kˆz . Sedimentary basins are typically horizontally stratified with horizontal permeabilities that are 1 or 2 orders of magnitude greater than vertical permeabilities. In addition, since a finite-length aquifer representing a portion of a plume is considered, the horizontal length scale may not be sufficiently larger than the vertical length scale depending on the discretization of the domain. Therefore, the vertical equilibrium assumption can be invalid and vertical flow effects may need to be considered. This is especially true in the vicinity of leakage pathways. If RL is sufficiently small, vertical flows may be important and the vertical equilibrium assumption inherent in the derivation of the 2-D vertically-integrated equations needs to be relaxed. Numerical simulations have shown that a linear structure with maximum vertical fluxes at the interface and no vertical flux at the top is a good approximation for qz,α [56]. In that work, the vertical flux structure becomes slightly nonlinear closer to the well in the denser fluid only. However, this is where the magnitude of the vertical flux, qz,d h(x,t) , is the greatest and any difference between the linear approximation and the numerically-obtained vertical flow structure is relatively small. A linear structure is the lowest-order interpolation between the two vertical fluxes. More complex structures can be employed but typically require the definition of additional parameters, which are difficult to quantify. Therefore, to account for vertical flow effects, we employ the lowest-order interpolation between no vertical flux at the top and bottom boundaries and the vertical flux at the interface to represent the vertical flux field in the z-component of equation (2.2) as a piecewise linear 21 interpolation: qz,d = qz,d qz,l = qz,l 1−z for z > h(x, t), 1 − h(x, t) z for z < h(x, t). h(x,t) h(x, t) h(x,t) (2.14) (2.15) At the interface, a relationship between horizontal and vertical fluxes can be derived based on the substantial derivative method or the mass conservation method [9, ] to get qz,α h(x,t) = qα ∂h h(x,t) ∂x + ∂h . ∂t (2.16) The first term on the right side of equation (2.16) reflects the effects of horizontal fluxes, while the second term represents the change in interface location over time. In other words, any vertical flux can be translated into a horizontal flux and changes in interface location. Substitution of equations (2.14), (2.15), and (2.16) into the z- component of equation (2.2) and integration over the thickness of the aquifer gives  k (1 − z)2 ∂h ∂h   + ρd z + qd h   Γ 2(1 − h) ∂x ∂t         k 1−h ∂h ∂h p(x, z, t) = p(x, 1, t) − ρd + qd h + h + ρl z +  Γ 2 ∂x ∂t          k h2 − z 2 ∂h ∂h   + ql h + Γλ 2h ∂x ∂t 22 for z > h, (2.17) for z < h. Vertical flow effects can be included in equations (2.12) and (2.13) by incorporating equation (2.17) to get ∂h ∂ ∂p k ∂ ∂ 1−h ∂h ∂h ∂h = h + h qd h + +h ∂t ∂x ∂x z=1 ∂x Γ ∂x ∂x 2 ∂x ∂t 2 2 ∂h ∂h 2h ∂ ∂h ∂h k ∂ h ∂h ql h + + ql h + + Γλ ∂x 3 ∂x ∂x ∂t 3 ∂x ∂x ∂t (2.18) and ∂h ∂p ∂ h = ∂t ∂x ∂x z=1 k ∂ (1 − h)2 ∂ + Γ ∂x 6 ∂x qd ∂h h ∂x + ∂h ∂t + 1 − h ∂h 6 ∂x qd ∂h h ∂x + ∂h ∂t (2.19) . Equations (2.18) and (2.19) can be solved numerically to obtain h = h(x, t). However, we choose to focus on solving the steady-state solutions to provide efficient solutions for inclusion in a basin-scale multi-scale model. Further, the steady-state solutions can be employed with the method of successive steady states [9] to capture changes in time. Equations (2.18) and (2.19) are presented in this section on general equations for completeness. 2.4 Solutions for leakage through faults Steady-state analytical solutions including vertical flow effects are derived for the problem of subsurface flow due to leakage through a fault. Two different steady-state leakage scenarios are considered: single-phase leakage of the lighter fluid only in section 2.4.1, and two-phase leakage in section 2.4.2. Single-phase leakage occurs when the leakage rate is small. In this case, Qd = 0. Two-phase leakage occurs when the leakage rate is larger. In this case, 0 < Qd < 1. For these two cases, we derive ordinary differential equations (ODEs) for the stationary solution accounting for vertical 23 flow effects. Simplified solutions are obtained by assuming vertical equilibrium and employing a small perturbation approach (Muskat’s approximation) [9]. 2.4.1 Single-phase leakage For sufficiently low leakage rates, only the lighter fluid will flow into the fault zone (Qd = 0) such that h qd = 0 and ql (x, z)dz = 1. (2.20) 0 In other words, the dense fluid is stagnant throughout the dense fluid zone (h ≤ z ≤ 1) including at the interface such that qd h = 0. Correspondingly, the pressure gradient in the x-direction in the dense fluid zone is zero and thus ∂p ∂x z=1 = 0. However, for the lighter fluid, the pressure gradient in the x-direction is non-zero and with the inclusion of vertical flow effects (equation (2.17)), we get dh h2 + z 2 k ∂p = + ql ∂x dx 2h2 Γλ 2 dh dx h + h2 − z 2 k d 2h Γλ dx ql dh h dx . (2.21) For z = h, the last term in equation (2.21) vanishes and the resulting expression is combined with equation (2.13) to yield ql h =− Γλ dh dx 1+k dh dx 2. (2.22) Substitution of equations (2.21) and (2.22) into equation (2.20) and subsequent integration gives the interface upconing solution for the stagnant dense fluid case, a second-order nonlinear ODE: −1 dh 2h =h −k Γλ dx 3 dh dx 1+k 3 dh dx 2 24 h2 −k 3 2 d h 2 dh dx dx2 1+k dh dx 2 2. (2.23) To solve equation (2.23), the interface location and the vertical equilibrium condition are specified at the outer boundary: h(x = xouter ) = h0 , (2.24) −1 . Γλh0 (2.25) dh dx = x=xouter The solution to equation (2.23) is obtained using a numerical solution technique (e.g. Runge-Kutta methods). This solution is referred to as the full ODE solution for the single-phase leakage case. In equation (2.23), vertical flow effects are represented in the last two terms. When k → 0, these last two terms vanish and equation (2.23) simplifies to a separable equation without vertical flow effects. Upon integration with the boundary condition (2.24), we get h2 = h20 − 2 (x − xouter ). Γλ (2.26) This equation is referred to here as the “vertical equilibrium” solution and is analogous to the Dupuit solution for upconing around a well [9]. Note that δ is included in the Γ term and is not necessarily small. As δ becomes large, the slope of h(x) becomes smaller; while as δ becomes small, the slope of h(x) becomes larger. Since h is a minimum at xinner , we can determine when equation (2.26) is applicable by setting hinner to its lowest possible value of zero to get Γ ≥ 2 xinner − xouter . h20 λ 25 (2.27) Equation (2.27) is useful in determining if leakage is single-phase or two-phase. For ˆ f ault values, which corresponds example, two-phase leakage is promoted for higher Q to lower Γ values. Another simplified solution can be obtained by employing a small perturbation approach where we assume h ≈ h0 to get h = h0 − 1 (x − xouter ). Γλh0 (2.28) This solution is analogous to Muskat’s approximation for upconing around a well [9] and is referred to here as Muskat’s approximation. Figure 2.2 shows a comparison of the full ODE solution, the vertical equilibrium ˆ f ault = −0.05 m2 /day given parameters solution, and Muskat’s approximation for Q presented in Table 2.1 and δ = 0.194. Two different k values, 26.7 and 426, are considered. As expected, the vertical equilibrium solution (dash-dot line) departs from the full ODE solution (solid line) as k increases and RL decreases. This is especially true near the fault (xinner = 0). For RL = RL2 = 1.0, the vertical equilibrium solution is essentially identical to the full ODE solution. This RL of 1.0 is lower than the criteria for vertical equilibrium suggested by [44] (RL > 10) and by [88] (RL2 1). Nonetheless, decreasing values of RL from 1.0 correspond to non-vertical-equilibrium conditions as demonstrated in the case with RL2 = 0.0625 << 1, where the vertical equilibrium solution departs noticeably from the full ODE solution near the fault. Muskat’s approximation also departs from the full ODE solution as x → 0 in the same manner and is the poorest approximation for the range of parameters explored here. Given the simplicity of the vertical equilibrium solution (equation (2.26)), there is no advantage to using Muskat’s approximation but it is included here for completeness. Figure 2.2 also shows that vertical flow effects are confined to a region near the fault. To approximate the location where vertical flow effects become negligible, we 26 0 0.05 0.1 0.1 0.15 0.15 0.2 0.2 0.25 0.25 0.3 0 0.2 0.4 x 0.6 0.8 Full ODE Vertical equilibrium Muskat’s approximation 0.05 h h 0 Full ODE Vertical equilibrium Muskat’s approximation 0.3 0 1 0.2 0.4 x 0.6 0.8 1 (b) (a) ˆ f ault = −0.05 m2 /day Figure 2.2: Analytical solutions for single-phase leakage given Q with (a) RL = RL2 = 1.0, and k = 26.7, and (b) RL = 0.25, RL2 = 0.0625, and k = 426. Note that xinner = 0 and xouter = 1. can balance the first term with the second or third term on the right side of equation (2.23) since the second and third terms represent vertical flow effects. If k small, xtransition ∼ 2k /3 h; while if k dh 2 dx is large, xtransition ∼ dh 2 dx is 3k /2 h. These steady-state transition locations, xtransition , can be viewed as a relative measure of √ vertical flow effects. The dependence of xtransition on k is in line with the discussion √ on RL in the previous paragraph, recalling that RL = δ k . It is important to note that there are other factors controlling the relative importance of vertical flow effects. ˆ f ault increases For example, vertical flow effects also become more significant as Q ˆ f ault , all three solutions converge. and for very small Q 2.4.2 Two-phase leakage At higher leakage rates, both the lighter and denser fluids contribute to the total leakage rate with 0 < Qd < 1, where 1 Qd = h qd (x, z)dz and Ql = h ql (x, z)dz. 0 27 (2.29) Note that Qd + Ql = 1. In contrast to the previous case on single-phase leakage, the pressure gradients in each fluid zone are non-zero and can be obtained by taking the derivative of equation (2.17) in the x-direction. Substituting the pressure gradients for each of the two fluids into equation (2.29) gives: Qd ∂p =− Γ(1 − h) ∂x ∂p Ql =− Γλh ∂x k − 6Γ z=1 k dh + 3 − dx 6Γ z=1 dh dx 2 dh dx 2 4 − λ d dx qd h + (1 − h) qd h − 3(1 − h) 2 dh dx ql h − d dx 2h d λ dx dh qd (2.30) h dx qd dh h dx dh ql (2.31) h dx The pressure gradient at the bottom of the aquifer (z = 1) is non-zero (since the denser fluid is flowing) and is unknown. To eliminate this unknown pressure derivative, equation (2.31) is subtracted from equation (2.30) to get dh k 2 Γλ + dx 3 dh dx 2 ql h − λqd h + λ(1 − h) d dx qd dh h dx +h d dx = ql dh h dx Qd λ Q − l(2.32) . 1−h h At the interface where z = h, equations (2.12) and (2.13) can be combined with the corresponding pressure derivatives to get ql h − λqd h =− Γλ 1+k dh dx . dh 2 dx The vertical equilibrium assumption is invoked to approximate qα (2.33) h only in the terms where the interfacial fluxes are contained within the derivative. Since vertical equilibrium is assumed, the pressure derivative dp/dx is independent of z and qd h ≈ Qd , (1 − h) and ql 28 h ≈ Ql . h (2.34) Substituting equations (2.33) and (2.34) into equation (2.32), we get an expression with h(x) as the only dependent variable: dh k Γλ + dx 3 + Ql − 2Γλ 1+ d2 h 1 − dx2 h dh 3 dx 2 k dh dx 2 dh dx + λQd = d2 h 1 + 2 dx 1−h dh dx 2 (2.35) Qd λ Q − l. 1−h h Using numerical methods, equation (2.35) can be solved for a given Qd and two outer boundary conditions, which are the interface height, h(x = xouter ) = h0 , and vertical equilibrium. It is important to note that the interface location at the inner boundary h(xinner ) = hinner can be specified instead of Qd since there is a unique Qd for each hinner . However, there is a range of possible Qd (or hinner ) values given a set of aquifer properties and thus, more information is needed to constrain the problem. This is in contrast to previous work of ref. [56] where hinner is assumed to always be zero in twophase leakage conditions. To further constrain the problem, numerical simulations are performed and the role of fault properties on Qd are explored in section 2.5. As with the single-phase leakage case, equation (2.35) can be simplified by assuming k → 0: dh Qd Q = − l . dx Γ(1 − h) Γλh (2.36) Equation (2.36) is referred to here as the vertical equilibrium equation. Alternatively, a simplified ODE can also be obtained if the interface is assumed to be sufficiently smooth with dh/dx being small so that higher-order derivatives are negligible. Then all second-order derivative terms, as well as products of first-order derivatives, are omitted to obtain a separable equation, which is identical to the vertical equilibrium expression presented as equation (2.36). The solution to equation (2.36) is the 29 following transcendental equation: F (h, x) = F (h0 , xouter ), (2.37) where −x h2 − Γλ 2(Qd λ + (1 − Qd )) Qd λ(1 − Qd ) hQd λ + ln ((1 − Qd )(1 − h) − Qd λh) . + 2 (Qd λ + (1 − Qd )) (Qd λ + (1 − Qd ))3 F (h, x) = The solution to equation (2.37) can be obtained by specifying a Qd value and applying the same two outer boundary conditions used to solve the full ODE. When La → ∞, the first term on the right side of equation (2.37) goes to negative infinity and thus, one of the terms on the left side must go the negative infinity. Therefore, the expression contained in the logarithm in the last term of F (h, x) must go to zero, which gives lim Qd = La →∞ 1 − h0 . 1 + h0 (λ − 1) (2.38) This limit also represents the upper value for Qd and is referred to here as Qd,max . The lower limit of Qd , for a given set of fluid and aquifer properties, can be obtained using equation (2.37) by setting hinner to its minimum value of 0. Finally, an analogous expression to Muskat’s approximation is found by setting k to zero and replacing h with h0 in equation (2.35) to get h = h0 + Qd 1 − Qd − Γ(1 − h0 ) Γλh0 (x − xouter ). (2.39) This equation is referred to here as Muskat’s approximation. In the two-phase leakage case, the advantage of Muskat’s approximation is that it provides a simple algebraic equation, which is useful in a multi-scale framework. 30 0 0.05 0.1 0.1 0.15 0.15 0.2 0.2 0.25 0.25 0.3 0 0.2 0.4 x 0.6 0.8 Full ODE Vertical equilibrium Muskat’s approximation 0.05 h h 0 Full ODE Vertical equilibrium Muskat’s approximation 0.3 0 1 0.2 (a) 0.4 x 0.6 0.8 1 (b) Figure 2.3: Analytical solutions for two-phase leakage with (a) Qd = 0.23 and (b) ˆ f ault = −0.2 m2 /day, RL = 5, R2 = 25, and k = 1.0. Note that Qd = 0.24, given Q L xinner = 0 and xouter = 1. 0 0.05 0.1 0.1 0.15 0.15 0.2 0.2 0.25 0.25 0.3 0 0.2 0.4 x 0.6 0.8 Full ODE Vertical equilibrium Muskat’s approximation 0.05 h h 0 Full ODE Vertical equilibrium Muskat’s approximation 0.3 0 1 0.2 0.4 (b) (a) x 0.6 0.8 1 Figure 2.4: Analytical solutions for two-phase leakage with (a) Qd = 0.22 and (b) ˆ f ault = −0.2 m2 /day, RL = 1.14, R2 = 1.31, and k = 20. Note Qd = 0.24, given Q L that the vertical equilibrium solution can give non-physical results. Also, note that xinner = 0 and xouter = 1. Figures 2.3 and 2.4 show comparisons of the full ODE solution to the two simplified ˆ f ault of 0.2 m2 /day given δ = 0.195 and parameters presented in Table solutions for a Q 2.1. The maximum Qd for both cases shown in these two figures is 0.24; while the minimum Qd values for Figures 2.3 and 2.4 are 0.23 and 0.22 respectively. In general, ˆ f ault . However, the range of Qd values for a given problem decreases with increasing Q the range in hinner values remain the same (between 0 and h0 ). 31 Comparing the full ODE solution (solid line) and the vertical equilibrium solution (dash-dot line) in Figure 2.3 shows that vertical flow effects are not significant for k = 1, which corresponds to RL = 5 and RL2 = 25. This is generally in line with the criteria by ref. [44] and ref. [88] for vertical equilibrium. As with the case for single-phase leakage, Muskat’s approximation (dashed line) provides a poor approximation when compared to the full ODE solution, which is especially true near the fault. For larger k values (Figure 2.4), the vertical equilibrium solution does not provide an improved approximation to the full ODE solution compared to Muskat’s approximation. In fact, when vertical flows are very significant (Figure 2.4a), the vertical equilibrium solution gives negative values of hinner . For leakage of CO2 to occur, 0 ≤ hinner ≤ 1 when Qd < 1 (since a sharp interface between the fluids is assumed). Therefore, the vertical equilibrium solution can give non-physical results if vertical flows are sufficiently important. In contrast to the single-phase leakage case, RL ≈ 1 is not sufficiently large for the vertical equilibrium assumption to be valid, as shown in Figure 2.4 for k = 20, RL = 1.1, and RL2 = 1.3. For all k values, vertical flow effects become less significant as Qd approaches Qd,max and hinner → h0 . However, the influence of Qd is greater when vertical flow effects are more significant. Furthermore, similar to the single-phase leakage case, vertical flow effects become more significant ˆ f ault increases, while all three solutions converge for very small Q ˆ f ault . as Q Figures 2.3 and 2.4 also show that mismatches among the solutions are greatest in the vicinity of the fault. Scaling relationships similar to those determined for the single-phase leakage case can be used to define a relative measure for the significance of vertical flow effects for two-phase leakage. The first and second terms in equation (2.35) can be balanced for small values of k dh 2 dx to give xtransition ∼ (3k /2) h. 32 (2.40) The first term in equation (2.35) can be balanced with the third and fourth terms in (2.35) to give xtransition ∼ k qd h k ql h h and xtransition ∼ h. 3Γ 3Γλ In the two-phase leakage case, xtransition is proportional to k , as opposed to (2.41) √ k in the single-phase leakage case. This implies that in the case of two-phase leakage, the anisotropy in aquifer permeabilities leads to greater vertical flow effects. This is in line with the discussion above where a larger RL value is required for the vertical equilibrium assumption to be valid in the two-phase leakage case. 2.5 Fault representation and leakage The analytical solutions for leakage through faults derived above do not explicitly include fault properties. They are defined given a Qd value (or hinner ) and two outer boundary conditions. Recall from equation (2.1) that at steady-state, Qd represents both the integrated dense fluid flow out of the aquifer into the fault and the integrated ˆ f ault . Therefore, Q dense fluid flow out of the top of the fault as a portion of Q d represents both aquifer and fault properties. To explore the combined effect of fault and aquifer properties, a 2-D numerical two-phase flow model is used to simulate vertical (ˆ z -direction) and horizontal (ˆ xdirection) flows in the aquifer and the fault zone. The 2-D model is solved using a fully-coupled fully-implicit approach. For the aquifer and fluid properties, values given in Table 2.1 are used along with k = 20 (kˆx = 400 mD and kˆz = 20 mD), which corresponds to RL = 1.1 and RL2 = 1.2. The outer boundary conditions are hydrostatic pressures and a fixed interface location, h0 = 0.25, which are identical to the boundary conditions applied to the analytical solutions. The conductive fault zone of finite width is explicitly represented by a region with kˆzf 33 kˆz and kˆxf ≥ kˆx . The ˆ f ault , is specified as the flow leaving from the top of the fault zone. The leakage rate, Q transient simulations with an initially flat interface are run to steady-state conditions, which is when Qd is constant. The numerical simulations are used to determine a ˆ f ault ); this is relationship between Qd and fault properties for a given leakage rate (Q presented in section 2.5.1. The leakage rate is considered as an input value here since ˆ f ault will be solved using information about pressure in a multi-scale framework, Q conditions in the overlying aquifer. In section 2.5.2, the numerical simulations are compared to analytical solutions to explore the effects of fault properties. 2.5.1 A fault model ˆ f kˆzf , and the Ref. [4] and ref. [5] identified the vertical fault transmissivity, Tˆzf = W anisotropy ratio, m, as key fault properties governing flow in the fault. Based on this and equation (2.1), the following fault model is proposed, Qd = Qd,0 + aT ln (T ) + am m, (2.42) ˆ f ault and aquifer properties (e.g. h0 ), where Qd,0 , aT , and am are constant for fixed Q and T is the ratio of vertical and horizontal transmissibility given by T = ˆ f kˆzf W . ˆ kˆx H (2.43) Note that T can be thought of as the ratio of the horizontal and vertical pressure gradients in the fault and can be obtained by rearranging the terms in equation (2.1). Results from a series of 2-D numerical simulations with various T and m values ˆ f ault and aquifer properties. are used to empirically obtain Qd,0 , aT , and am given Q The simulations represent a range of fault properties for fixed aquifer properties. The ˆ f , and 400 mD to 10 D range of fault properties considered are 0.2 m to 1.0 m for W for kˆzf and kˆxf . The m values considered range from 1 to 10, which correspond to 34 Q'd from numerical simulations 0.220 0.215 0.210 0.205 y = 1.000x R² = 0.995 0.200 0.195 0.195 0.200 0.205 0.210 0.215 Q'd estimated using the fault model 0.220 Figure 2.5: Comparison of Qd values given by the numerical simulation, representing various T and m values, and the fault model (equation (2.42)) for Qd,0 = 0.223, aT = 0.009, and am = 8 × 10−5 . Note that xinner = 0 and xouter = 1. limits based on field data identified in ref. [80]. Figure 2.5 shows that equation (2.42) with empirical coefficients can reproduce Qd obtained from the numerical simulations. Such fault models can be used to determine the value of Qd to be used in the analytical solutions. This empirical relationship is useful in a multi-scale framework since the relationship can be determined a priori with only a few numerical simulations with a domain representing only a single grid-block. While a more comprehensive fault model representing a wider range of conditions could be pursued, it is beyond the scope of the current paper. 2.5.2 Effect of fault properties The numerical simulations are compared to the corresponding analytical solutions with the same Qd value. Here, the effects of the two key variables in the fault model, T and m, are explored for both the single- and two-phase leakage cases. For converged numerical simulations, run times of 5-8 hours with a MATLAB code are required. This is in contrast to several seconds required to solve ODEs of the analytical solutions in MATLAB. In the case of single-phase leakage, Figure 2.6 shows that the full ODE and vertical equilibrium solutions match the numerical simulation results well. As expected, 35 0 0.05 0.1 0.15 0.1 0.15 0.2 0.2 0.25 0.25 0.3 0 0.2 0.4 (a) x 0.6 0.8 Full ODE Vertical equilibrium Muskat’s approximation Numerical result 0.05 h h 0 Full ODE Vertical equilibrium Muskat’s approximation Numerical result 0.3 0 1 0.2 0.4 (b) x 0.6 0.8 1 Figure 2.6: Numerical and analytical solutions for single-phase leakage with (a) T = ˆ f ault = −0.05 m2 /day and m = 1. 0.04 and (b) T = 0.2, given Q Muskat’s approximation departs significantly from the numerical simulation results, especially near the fault. The numerical results are relatively insensitive to T for ˆ f ault = −0.05 m2 /day. A slightly better match values in the range of 0.04 to 0.2 for Q is obtained with the larger T value. No changes in the numerical simulation results are observed as T is increased further and thus, the analytical solution is applicable when T > T1critical ≈ 0.2 (where the subscript “1” denotes single-phase flow). In the two-phase leakage case, Figure 2.7 shows that the full ODE solution provides a significantly better approximation to the numerical simulation results when compared to Muskat’s approximation or the vertical equilibrium solution. The full ODE solution is unable to fully capture the flow field near the fault when T values are below 0.8 ˆ f ault = −0.2 m2 /day (where the subscript “2” denotes and thus, T2critical ≈ 0.8 for Q two-phase flow). In both the single- and two-phase leakage cases, the analytical solutions appear to be relatively insensitive to m within the range of typical values of 1 to 10, as shown in Figures 2.8 and 2.9 respectively. Nonetheless, in the two-phase leakage case, higher m values correspond to a slightly improved match (Figure 2.9) and further improvements may be obtained for larger values of m, which have been shown to be 36 0 0.05 0.1 0.15 0.1 0.15 0.2 0.2 0.25 0.25 0.3 0 0.2 0.4 x 0.6 0.8 Full ODE Vertical equilibrium Muskat’s approximation Numerical result 0.05 h h 0 Full ODE Vertical equilibrium Muskat’s approximation Numerical result 0.3 0 1 0.2 0.4 x 0.6 0.8 1 (b) (a) Figure 2.7: Numerical and analytical solutions for two-phase leakage with (a) T = 0.2 ˆ f ault = −0.2 m2 /day (and Qd = 0.209) and (b) T = 0.8 (and Qd = 0.217), given Q and m = 1. Note that xinner = 0 and xouter = 1. 0 0.05 0.1 0.15 0.1 0.15 0.2 0.2 0.25 0.25 0.3 0 0.2 0.4 x 0.6 0.8 Full ODE Vertical equilibrium Muskat’s approximation Numerical result 0.05 h h 0 Full ODE Vertical equilibrium Muskat’s approximation Numerical result 0.3 0 1 0.2 0.4 x 0.6 0.8 1 (b) (a) Figure 2.8: Numerical and analytical solutions for single-phase leakage with (a) m = 1 ˆ f ault = −0.05 m2 /day and T = 1.0. Note that xinner = 0 and and (b) m = 10, given Q xouter = 1. plausible [5]. Therefore, the critical value for m, above which the analytical solution and numerical simulation results match, appears to be greater for two-phase leakage than single-phase leakage. Numerical results with higher m values do show stronger vertical flow in the fault. This is in line with ref. [5], where small m values led to more horizontal flow in the fault, while large m values led to predominantly vertical flow in the fault. 37 0 0.05 0.1 0.15 0.1 0.15 0.2 0.2 0.25 0.25 0.3 0 0.2 0.4 (a) x 0.6 0.8 Full ODE Vertical equilibrium Muskat’s approximation Numerical result 0.05 h h 0 Full ODE Vertical equilibrium Muskat’s approximation Numerical result 0.3 0 1 0.2 0.4 (b) x 0.6 0.8 1 Figure 2.9: Numerical and analytical solutions for two-phase leakage with (a) m = 1 ˆ f ault = −0.2 m2 /day (and Qd = 0.218), and (b) m = 10 (and Qd = 0.219), given Q and T = 1.0. Note that xinner = 0 and xouter = 1. 2.6 Conclusions The steady-state analytical solutions for flow in an aquifer adjacent to a fault with the capability to represent vertical flow effects, coupled with an empirical model for the fault, provide a good model for vertical leakage through faults. This is true for both single- and two-phase leakage through faults. In the single-phase leakage case, vertical flow effects are important when RL = RL2 ≥ 1. In the two-phase leakage case, anisotropic fault permeabilities lead to greater vertical flow effects and a larger RL value, of 5 or greater, is needed to assume vertical equilibrium. Fault properties, including fault width and permeabilities in both the horizontal and vertical directions, strongly govern two-phase leakage through faults including the fraction of dense fluid that leaks, Qd . In fact, fault properties are necessary to obtain a unique solution for Qd . Realistic T and m values are generally sufficiently large for the analytical solutions to be applicable. However, the critical values of T and m, above which the analytical solutions are applicable, are greater in the case of two-phase leakage. Comparisons of the analytical and numerical simulation results show that the steady-state analytical solutions along with a fault model can be used in place of fine38 scale numerical simulations in the vicinity of a leaky fault. The analytical solutions along with the fault model can be obtained in the order of seconds as opposed to hours or tens of hours for converged numerical simulations. The development of these solutions satisfies the motivation of this chapter, which is to provide an efficient approach to include sub-grid-scale leakage through faults in basin-scale models. 39 Chapter 3 A multi-scale model for leakage through faults using combined analytical-numerical methods 3.1 Introduction A major motivation and application of the work on models of leakage through faults presented in Chapter 2 is for its use as a fine-scale representation in a multi-scale model. The inclusion of the analytical solutions and fault model allow for an accurate and efficient description of the leakage through faults and its role in the surrounding formation. In this chapter, the corresponding multi-scale framework is presented. The multi-scale framework considers perpendicular flow to the fault and parallel flow through the fault (both vertical and horizontal). Abandoned wells are point features in a two-dimensional (2-D) space and the effect of fluid leakage through abandoned wells can be represented radially. In contrast, faults are linear features in a 2-D space and can have lengths in the order of kilometers. Therefore, faults often extends over multiple grid blocks in typical basin-scale models. In considering 40 a multi-scale framework for leakage through faults, new considerations for horizontal flow in the fault within a coarse-scale grid block and between coarse-scale grid blocks, which are not made for well leakage, must be made. These flows are represented in the grid block containing the fault and the aquifer adjacent to the fault using fine-scale pressure and fluid distributions. The fine-scale pressure and fluid distributions hinge on the horizontal flow field perpendicular to a leaky fault, which is the topic covered in Chapter 2. Therefore, the model is capable of representing both single-phase and two-phase leakage through faults, along with any vertical flow effects, and the effects of fault properties. The multi-scale model developed here focuses on a single formation in which two fluids (e.g. CO2 and brine) are present. This single-layer model can be used to represent one of many stacked layers that represent the sequence of geologic formations in a basin of concern [60]. This single layer is represented by a 2-D vertically-integrated numerical model at the coarse scale and by a system of nonlinear equations at the fine scale. The fine scale model represents horizontal flow perpendicular to the fault in the aquifer, vertical flow through faults, and horizontal flow through faults. In this chapter, the coarse-scale numerical model (section 3.2), the various components of the fine-scale model (section 3.3), and a brief description of how these components interact (section 3.4) are presented. In section 3.4, a numerical example of a multi-scale model implementation for one simple case in the context of geologic storage of CO2 is presented. 3.2 Coarse-scale numerical model The coarse-scale numerical model represents 2-D vertically-integrated two-phase flow through porous media [24, 57]. The assumptions made here is that vertical equilibrium conditions exist, the fluids are incompressible, capillary pressures are negligible, and 41 the properties of the storage formation (aquifer) are homogeneous. Applying these assumptions to the mass balance equations presented as equation (2.2), we get ∇ · UΣ = ΥΣ , (3.1) where ΥΣ represents the sum of vertically-integrated source/sinks of the two fluids [LT−1 ], and UΣ represents the sum of Uα , and α is l or d for the lighter and denser fluids respectively. Assuming a sharp interface between the two fluids, Uα is the vertically-integrated (horizontal) flux [L2 T−1 ] for fluid α as follows: ¯ h k ¯ (∇P ) (H − h) µd 0 H k¯ ¯ , Ul = ul dz = − h ∇P − ∆α ρg∇h µl ¯ h Ud = ud dz = − (3.2) (3.3) ¯ is the coarse-scale location of the interface between the lighter and denser where h fluids [L], k is the intrinsic permeability of the aquifer [L2 ], H is the thickness of the aquifer [L], P is the pressure at the bottom of the aquifer [ML−1 T−2 ], ∆α ρ = ρd − ρl , ρd and ρl are the densities of the denser and lighter fluids respectively [ML−3 ], g is the acceleration due to gravity [LT−2 ], and z is the vertical coordinate [L]. This 2-D model can be solved using the IMplicit Pressure Explicit Saturation (IMPES) method. Therefore, for a given time step, equation (3.1) is solved implicitly for pressure and then this information is used to explicitly solve equations (3.2) or ¯ The pressure, P , and the thickness of the lighter fluid h ¯ in the grid block (3.3) for h. containing a fault are passed on to the fine-scale model, in which this coarse-scale information is used as boundary conditions. 42 3.3 Fine-scale model We consider a fault with a single impermeable core and damage zones on both sides of the core [21]. We seek to represent flow through the damage zones (also referred to as fault zones), which can be of higher permeability than the surrounding aquifer and can extend through overlying layers of low and high permeability. Vertical leakage through faults occurs when a fault extends from a storage formation (or aquifer) through a low-permeability layer (or caprock) into an overlying the aquifer. We are especially concerned with fault zones (simply referred to here as faults) with permeabilities that are orders of magnitude greater than the aquifer permeability. Flow in the fault zone on one side of an impermeable fault core is considered implicitly within a coarse-scale grid block containing a fault since the fault width is assumed to be small relative to the coarse-scale grid block width. The fault is conceptualized to be located at an edge of a rectangular grid block, and given an impermeable fault core, this edge is defined to have no flow. Therefore, the horizontal flow perpendicular to the fault can be described using the analytical solutions and fault model presented in Chapter 2. The analytical solutions and fault model are used to determine the pressure and fluid variations in the direction perpendicular to the fault. The pressures and fluid distribution at the fault, coupled with the pressure in the overlying aquifer, drives vertical leakage from the formation and is used to determine the horizontal along-fault flow in the fault. Therefore, equations for pressure corrections and fluid distributions based on the flow field perpendicular to the fault are derived. These corrections are designed to be applied to the coarse-scale pressures and saturations, which are assumed to be an average over the coarse-scale grid block. Considerations for how these fine-scale components can be represented within and between coarse-scale grid blocks are also discussed. 43 3.3.1 Horizontal flow perpendicular to fault Analytical solutions that represent horizontal flow perpendicular to a fault are presented for single-phase and two-phase leakage conditions. The steady-state solutions describe the fluid-fluid interface location, h(x), between the two fluids and captures the effects of a leaky fault. The full solution including vertical flow effects can be solved or simplified solutions can be used when appropriate. In the case of vertical equilibrium, the vertical equilibrium solutions can be used; while in the case of small-leakage, the Muskat’s approximations can be used. Depending on the aquifer and fault properties, the simplest but accurate solution should be used. An algorithm for the selection of the solution can be designed based on the work presented in Chapter 2. The corresponding interface, h(x), can be used to determine the pressure distribution in the direction perpendicular to the fault. Fault properties can be included by an empirical relationship, whose coefficients can be defined by numerical simulations of the fine-scale domain. Numerical simulations of 29 cases with different leakage rates, Qf ault , and thicknesses of the lighter fluid at the outer boundary, h0 , are shown in Figure 3.1 using the properties outlined in Chapter 2 for a fault with T = 1.0 and m = 1. The corresponding RL2 is 4.8, which should be high enough for the vertical equilibrium solution to be valid in the singlephase leakage case and sufficiently high for the vertical equilibrium solution to be valid in the two-phase leakage case. The relationships shown in Figure 3.1 are based on numerical simulations in the x- and z-dimensions so no conditions on vertical flow are specified. Figure 3.1 shows that there are two important transition leakage rates: one is when the denser fluids begin to flow such that we have two-phase leakage and the other is when Qd reaches a maximum. These transition leakage rates may also be derived analytically; however, the changes in slope between these transition leakage rates are determined using numerical simulations. Therefore, the overall function, 44 Qd (Qf ault , h0 ), is determined based on the both numerical and analytical solutions. Considering the first transition leakage rate, Figure 3.1 shows that the denser fluid begins to leak for a smaller Qf ault when h0 is smaller, which is evident from equation (2.27). The leakage rate at which Qd transitions to maximum Qd occurs at a lower leakage rate when h0 is smaller. Furthermore, h0 is inversely related to the maximum Qd (as described in equation (2.38)), and more of the denser fluid leaks for thinner plumes. 0.6   Q'd   0.5   0.4   0.5  (empirical  model)   0.3   1.25  (empirical  model)   0.2   0.5  (numerical  model)   0.1   1.25  (numerical  model)   0.0   1.E-­‐03   1.E-­‐02   1.E-­‐01   1.E+00   Qfault  [m2/d]     1.E+01   1.E+02   Figure 3.1: The empirical relationship between Qd and Qf ault for h0 = 0.5 and 1.25 m determined using the fault model and numerical simulations for T = 1.0 and m = 1. The numerical simulations used to determine coefficients of the fault model and validate it are also shown. It is only between these two transition leakage rates that fault properties play a role in the flow field. Therefore, only a few numerical simulations are required to determine the effects of fault properties in the case presented in Figure 3.1. In many cases, it may be sufficient to interpolate between the two transition points. The need to perform numerical simulations and the number of simulations required are based on the relative values of these two transition points. 45 3.3.2 Vertical leakage through fault Vertical flux through faults is represented using the multi-phase extension of Darcy’s Law qfαault kf ault kr,α (Sfαault ) =− µα pf ault − ptop − ρα g ∆zcap (3.4) where pf ault is the pressure in the fault at the top of the aquifer [ML−1 T−2 ], ptop is the pressure in the fault at the bottom of the upper aquifer [ML−1 T−2 ], ρα is the density of fluid α [ML−3 ], kf ault is the vertical permeability of the fault [L2 ], kr,α (Sfαault ) is the relative permeability of fluid α as a function of the saturation of fluid α in the fault (Sfαault ), and ∆zcap is the thickness of the caprock separating the two aquifers [L]. If both the lighter and the denser fluids are flowing, qf ault = qfdault + qfl ault [LT−1 ], and Qd = qfdault /qf ault [-]. Sfdault and Sfl ault , which sum up to 1.0, are unknowns. Equation 3.4 describes the flow in the fault from the top of the formation being modeled and the bottom of an overlying formation. We focus on leakage scenarios only in this chapter. The leakage rate (qf ault ) is a key unknown. To determine this, vertical flow through the fault is defined by the pressures above and below the caprock, which is solved simultaneously with horizontal flow and pressure variations perpendicular to the fault. 3.3.3 Pressure correction Pressure corrections are made in grid blocks containing faults to determine the pressure variation in the direction perpendicular to the fault. This information is also used to determine the horizontal flow parallel to the fault, in the fault and in the adjacent aquifer. 46 The average grid block pressure, P or p¯coarse (ˆ z ) [ML−1 T−2 ], can be corrected to represent the pressure decrease near the leaky fault as follows: pˆcorrected (ˆ x, zˆ) = p¯coarse (ˆ z ) − ∆ˆ p(ˆ x, zˆ), (3.5) where pˆcorrected is the pressure corrected to account for effects due to leakage [ML−1 T−2 ], and ∆ˆ p is the correction term. In general, variables denoted with a hat are dimensional. Pressure corrections are derived from the analytical solutions presented in Chapter 2. The correction term is derived analytically in dimensionless form using xouter p(x1 , z)dx1 ∆p(x, z) = p¯f ine (z) − pf ine (x, z) = x xouter − x − pf ine (x, z), (3.6) where p¯f ine is the average of the fine-scale pressure within a grid block [-], pf ine is the pressure representing local-scale effects of leakage determined by the analytical solutions [-], and xouter is the location of the outer boundary [-]. The variables follow the non-dimensionalization presented in Chapter 2 (section 2.3.2). Equation (3.6) can be applied to cases where vertical flow is important and cases where vertical equilbrium exists. An algebraic equation for ∆p can be determined for both the vertical equilibrium solution and Muskat’s approximation in the single-phase leakage case, and only for Muskat’s approximation in the two-phase leakage case. The obvious advantage of deriving these algebraic solutions is that they greatly reduce computational demand and should be used when appropriate. Since pressure corrections are used to determine the leakage rate to the overlying aquifer, we derive these algebraic expressions for the pressure correction term at the fault. In the single-phase leakage case, the pressure correction term is only applicable to the region consisting of the lighter fluid only (0 < z < h) since the dense fluid 47 is stagnant. For single-phase leakage, the pressure at the top of the aquifer can be defined as pf ine (x, z = 0) = p¯ z=1 − ρd + h(x) (3.7) where p¯ at z = 1 is the same for both fine and coarse scales. (The pressure equation requires knowledge of the interface locations, h(x), which can be obtained using the equations for horizontal flow perpendicular to the fault.) Combining equations (3.7) and (3.6), the corresponding pressure correction terms are ∆pV,1 (x = xf ault , z = 0) = − ∆x , 2Γλh0 Γλ ∆pM,1 (x = xf ault , z = 0) = − 3∆x − h20 h30 (3.8) − h20 2∆x + Γλ 1/2 2∆x + Γλ 3/2 , (3.9) where the subscripts V and M represent vertical equilibrium and Muskat’s approximation respectively, the subscript 1 identifies the expression as corresponding to single-phase leakage, and ∆x = xouter − xf ault . In the two-phase leakage case, the pressure correction term is applicable for all z and the pressure correction solution can be derived analytically using ∂p Qd −(1 − Qd ) (x, z = 1) = − or ∂x Γ(1 − h(x)) Γλh(x) (3.10) following equations (2.30) and (2.31). We can analytically derive the pressure correction term assuming Muskat’s approximation to get a2 Qd (1 − h0 ) ln(1 − h0 ) Γ∆x a2 Qd − (1 − h0 + ∆x/a) ln (1 − h0 + ∆x/a) Γ∆x aQd − (1 + ln(1 − h0 + ∆x/a)) (3.11) Γ ∆pM,2 (x = xf ault , z) = − 48 where the subscript 2 denotes that the equation is for two-phase leakage and a= Qd Qd − . Γ(1 − h0 ) Γλh0 (3.12) The full ODE (including vertical flow effects) and the vertical equilibrium forms of the pressure correction solution for two-phase leakage are not derived analytically since an expression for h(x) alone are not available. However, they can be determined by solving the ODE and the transcendental equation, and then solving equation (3.6) with equation (3.10). 3.3.4 Horizontal flow parallel to fault The horizontal flow parallel to the fault within a grid block can be separated as the flow in the fault and the flow in the adjacent aquifer. For a fault that extends over multiple grid blocks, the flow within the fault should be described by pressures inside the fault (“2f” and “3f” in Figure 3.2). Similarly the flow between the aquifer portion of the grid blocks should be defined by the pressures in the aquifer (“2a” and “3a” in Figure 3.2). Weighing the flows by the proportion of the grid block volume occupied by the fault and the aquifer, the horizontal flow from grid block 2 to 3, for example, is U2,3 = U2f,3f · ∆xf ault + U2a,3a · (∆x − ∆xf ault ) ∆x (3.13) for rectangular grid blocks. Between a grid block with and without a fault at the end of a fault (e.g., grid blocks 1 and 2 in Figure 3.2), the horizontal flow from grid block 1 to 2, for example, is U1,2 = U1,2f · ∆xf ault + U1,2a · (∆x − ∆xf ault ) . ∆x 49 (3.14) 3f 3a 2f 2a y 1 x Figure 3.2: A schematic of grid blocks containing faults and adjacent grid blocks. The grey blocks represent coarse-scale grid blocks without a leaky fault and the white grid blocks represent the coarse-scale grid blocks with a leaky fault. Black dots represent the center of the grid block and the grey dots represent the center of the leaky fault within a grid block. To determine the horizontal flow in the fault and in the aquifer, the pressures and saturations of the fault and the aquifer in each grid block is needed. The pressures in the aquifer can be taken to be the average pressure of the grid block; while pressures inside the fault can be determined using the pressure corrections presented in section 3.3.3. To summarize, pia = p¯i (3.15) pif = p¯i − ∆pi (3.16) where subscript i is the grid block index, subscript a stands for aquifer, and subscript f stands for fault. We also need to determine h in the fault and in the aquifer since our coarse-scale model is vertically integrated. Therefore, corrections for h are needed to determine the portion of the flow that is through the fault. Corrections for h can be determined by ¯ f ine − hinner = ∆h = h 1 ∆x 50 ∆x 0 h(x )dx − h(xinner ). (3.17) Similar to pressure corrections, corrections for h must be derived for single-phase and two-phase leakage using vertical equilibrium and/or Muskat’s approximation. For the single-phase leakage case using Muskat’s approximation, we get ∆hM,1 Γλ =− 3∆x 3/2 h20 h20 − 2 + ∆x Γλ 3/2 − h20 + 2 ∆x. Γλ (3.18) For two-phase leakage case using Muskat’s approximation, we get ∆hM,2 = ∆x 2 Qd 1 − Qd − Γ(1 − h0 ) Γλh0 . (3.19) Therefore, the horizontal flow parallel to the fault is described at the coarse scale with corrections for the presence of the fault. 3.4 Multi-scale model The multi-scale model for the coarse-scale and fine-scale models described above requires the solution of the fine-scale model at every coarse-scale time step. The finescale components can be included as an intermediate calculation within a time step or be integrated into the coarse-scale equations. The equations are developed such that the vertical flow through faults and horizontal flow perpendicular to the fault must be solved simultaneously as an intermediate calculation. An algorithm that identifies the problem as single- and two-phase leakage and evaluates the significance of vertical flow effects is implemented to determine the appropriate equations. The solution of this nonlinear system of equations, which gives the leakage rate and the proportion of the leakage that is brine, is fed into the coarse-scale numerical model in the grid block containing the fault. In contrast, the horizontal flow parallel to faults are determined by modifications to the coarse-scale pressure and saturation 51 equations. These modifications can be viewed as acting like a relative permeability term at the coarse scale. Considerations for the appropriate outer boundary length and the ability of the fine-scale model based on steady-state solutions to represent the changes within a given time step are not considered here. 3.4.1 Numerical example The example presented here is in the context of CO2 storage in saline formations and shows CO2 injection in one grid block and a fault that provides a leakage pathway for CO2 and brine in another grid block. For the fine-scale model, only horizontal flow perpendicular to the fault and vertical flow through the fault are considered and pressure corrections derived assuming Muskat’s approximation are used. The horizontal flow perpendicular to the fault considers fault properties and uses the empirical relationship presented in Figure 3.1. Therefore, the fault and aquifer properties are as described in section 3.3.1. The 2-D domain shown in Figure 3.3 covers an area of 1500 m by 1500 m and consists of 100 coarse-scale grid blocks with dimensions of 50 m by 50 m. The outer boundary condition is constant pressure at 100 bars. The initial conditions are uniform pressure at 100 bars and no CO2 in the system. As shown in Figure 3.3, the injection location is above and to the right of the leaky fault that is 50 m long. The leaky fault extends throughout the length of the grid block and is on the left side of the coarse-scale grid block. As shown in Figure 3.3, the pressure distribution around the grid block with the leaky fault is asymmetrical because the impermeable core is represented by a no-flow condition at the face adjacent to the fault. The saline aquifer becomes saturated with CO2 near the injection site and a CO2 plume is formed. The CO2 plume is asymmetrical with larger saturations near the grid block containing a leaky fault. 52 Pressure [bar] Thickness of Lighter Fluid [m] 1400 1400 100 1000 95 800 90 600 4 1200 Horizontal Distance [m] Horizontal Distance [m] 1200 4.5 400 3.5 1000 3 800 2.5 600 2 1.5 400 1 85 200 200 500 1000 Horizontal Distance [m] 1500 0.5 500 1000 Horizontal Distance [m] 1500 Figure 3.3: Pressures (left) and saturations (right) at 7 years for the case with injection in one grid block and a leaky fault in one grid block. 1 70 0.9 60 0.7 40 Q’d Leakage Rate [m2/day] 0.8 50 0.6 0.5 30 0.4 20 0.3 10 0 0 0.2 2 4 6 Time [yr] 8 0.1 0 10 2 4 6 Time [yr] 8 10 Figure 3.4: The leakage rate, Qf ault , (left) and the proportion of the leakage that is brine, Qd , (right) with respect to time for the case with injection in one grid block and a leaky fault in one grid block. Figure 3.4 shows how the leakage rate increases rapidly in the first 2 years but remains at approximately 60 m2 /day until the end of the simulation time. Similarly, 53 Qd drops rapidly in the first two years and then stabilizes at around 0.2 for the remainder of the simulation time. 3.5 Conclusions The ability to capture the pressure variations and fluid distribution within a grid block containing a leaky fault allows for the representation of horizontal flow perpendicular to the fault, horizontal flow through the fault, and vertical flow through the fault. Fault properties are represented using an empirical model in the horizontal direction perpendicular to the fault and in the Darcy flux terms for horizontal flow parallel to the fault and vertical flow through the fault. The approaches for including these effects within a coarse-scale grid block and between coarse-scale grid blocks with and without leaky faults allow for a multi-scale model to be developed. Future work includes validation and testing of the multi-scale model under a wide range of scenarios such as faults that extend over multiple grid blocks and representation of cross-fault flow. Also, cases when faults do not extend throughout the fault or are not aligned with a side of a grid block need to be considered. The equations and concepts presented in this chapter can be used to develop a multi-scale model that represents leakage through faults in basin-scale models. Such basin-scale models can also include leakage through other subsurface pathways such as abandoned wells. In considering leakage at the basin scale, both of these pathways should be and can be represented by combining the methods presented here and in ref. [24], which provides a multi-scale framework for leakage through abandoned wells. 54 Chapter 4 Direct measurements of methane emissions from abandoned oil and gas wells in Pennsylvania 4.1 Introduction Methane emissions from abandoned oil and gas wells are estimated to be the second largest potential contribution to total U.S. methane emissions above U.S. Environmental Protection Agency estimates and are not included in any emissions inventory [11]. Further, there are no measurements that can be used to estimate the methane emission potential of these wells [11]. Methane is a greenhouse gas (GHG) and its oxidation produces ozone (O3 ) that degrades air quality and adversely impacts human health, agricultural yields, and ecosystem productivity [79]. Therefore, it is important to understand methane emission sources so that appropriate mitigation strategies can be developed and implemented. Efforts to improve estimates of methane emissions to the atmosphere from oil and gas production in the U.S. are being driven by growth in unconventional production. 55 Estimates of methane emissions from activities on producing oil and gas sites, such as well completion, routine maintenance, and equipment leaks, are used to develop bottom-up estimates [31, 1]. However, a comparison of bottom-up and top-down estimates indicates that there may be missing sources in bottom-up estimates [32, 69, 52, 11]. Here, we focus on one missing source: abandoned oil and gas wells. There is no regulatory requirement to monitor or account for methane emissions from abandoned wells in the U.S. Methane leakage through abandoned wells linked to recent growth in unconventional oil and gas production is being studied as a groundwater contamination issue [62, 26, 39, 38, 53], but no direct evidence for leakage through abandoned wells to groundwater aquifers currently exists. Abandoned wells have been connected to subsurface methane accumulations that have caused explosions, which are major concerns in urban areas with oil and gas development or natural gas storage reservoirs, as well as in coal mines [27, 14]. Therefore, any monitoring is focused on detecting large concentrations. The result is a lack of information to quantify methane emissions from abandoned oil and gas wells. To characterize abandoned oil and gas wells’ potential as a significant methane source, we made first-of-a-kind direct measurements of methane fluxes from 19 wells in various locations across McKean and Potter Counties in Pennsylvania (PA) (Figure 4.1). As shown in Figure 4.1, only 1 of the 19 wells is on the PA Department of Environmental Protections (DEP’s) list of abandoned and orphaned wells. (Orphaned wells can be defined as abandoned wells with no responsible party available, other than the state.) This is indicative of the general scarcity of available information on this class of old wells. Given the lack of records on the wells we measured, no distinction was made between oil and gas wells; the wells were simply categorized as plugged or unplugged, based on surface evidence of cementing and/or presence of a marker. With this criterion, 5 of the 19 measured wells (26%) were classified as plugged. 56 McKean County Potter County # ## # # # # # # # "# ## " # "" # ## # # # # # # # # # # # ## # # """" # # # "" # # # # ## # # # # ## # # " # # # #""# # # # # # "" " Foster # # " # # # # # # " # # # " " # # # " ## # # # ### # " """ # # # " "" # " # # " # # # " # # # " # # # # " " # # # # # " # # " # # # # # # "# #" # # # # # " # # # # # # # # # # # # # " # " # # # " " # # " Oswayo " " # " " " " "# # Ceres # Genesee # # # """"" " " " # " " " # " " # " " # # # Corydon "## # # "" # # # " # # # # #" " # " # # #"# # # ## # # "" # # ## # # # # # # Otto # # # " # #"" "## # # # # "" " # " Sharon ## # # # # # # # # # # # # # " # # # # " " # # # # " # " " " Annin # ## "# " # ### # # # " " # #"#" ## " # # Eldred # " # " "" # # # # # " # # Bradford # " # " " # "# # # # # " " # # # # "" " # # " " " # # # " " # Glade # # # # # # # # " " " # # "" # # # # Allegany " # " # # # # # " " # # " "# " # # Legend " "" "" " " " "" " # " " Keating Hebron " " # "" " # # " Measured Wells # # # # Mead # Eulalia ## Lafayette # Hamilton ## ## # " ## Liberty # # # " # # # # "# # # # " ## # DEP Abandoned List # # # # # # # "" # Coudersport ## Sweden # # ##### # # # "" # # " # ## " # # # DEP Orphan List ## " # # # # # # ### "" ## # # " # #" # # # # # # " # # # # # # # # ## # # # # # # " # ## DEP Plugged List # # # # # " Homer # # ## # # Hamlin ## # # # # ## # # # # # # # # # # # # Allegheny National Forest # # # # # # Sheffield # # ## # # # ## # # ## " # # # # # # ## # # # " # # # Summit # # # ## # # ## # # # # ## # # # "" # Keating # # # # Norwich # # PA State Game Lands # " " # #" # Sergeant " ## " Wetmore # " # # ( ! ! ( ( ! ! ( Pennsylvania ( ! ! ( ! ( ( ! ( ! ( ! # # # # " # # # # ## # # # ### ### # # ### # # # ## # " Jones 10 Kilometers ## # ## # ## # " " " 0 5 Highland Land Cover Forest Grassland River Wetland Total Number of Locations Wells Controls 9 23 5 10 1 1 4 8 19 42 #### Portage "" Portage Location Allegheny National Forest City of Bradford Otto Township Hebron Township Total Sylvania # " # # Wharton County Ownership McKean Private Public Private Private McKean McKean Potter Public Townships Counties Number of Wells Unplugged Plugged 2 1 1 7 3 14 3 0 0 2 0 5 On PA DEP List 1 0 0 0 0 1 Figure 4.1: The 19 measured wells are located in McKean County and Potter County in Pennsylvania. The map also shows the location of the 12,131 abandoned, orphaned, and plugged wells on the Pennsylvania Department of Environmental Protection (DEP)’s website as of January 17, 2014. Note that only the western portion of Potter County is shown in the detailed map. 57 In addition to methane, we also analyzed the samples for ethane, propane, and n-butane, and carbon isotopes of methane, to provide insight on the potential sources of the emitted methane. This work provides previously unavailable data on methane fluxes and other emissions from abandoned oil and gas wells. 4.2 4.2.1 Methods Site selection We selected abandoned oil and gas wells for measurement based on location information, access (legal and logistical), wellhead geometry above ground, land cover, geographic coverage, and plugging status. Site visits were made to confirm well locations and evaluate logistical access issues. Wells on both public lands and private properties were considered. Private properties were only considered if the surface landowner granted access since flux chambers were designed to enclose and not touch the well. Effort was made to ensure that wells in different land cover areas were measured and that we had broad geographic coverage. We measured 19 wells over 5 sampling rounds (Table 4.1) in McKean and Potter Counties in PA. The measured wells were located on four different land cover types: grassland, wetland, river, and forest (Figure 4.1). Five of the measured wells were located in the Allegheny National Forest. Nine were located on a 40-acre private property in Otto Township. The remaining two wells, one on a private property and one in the West Branch Tunungwant Creek, were in the City of Bradford. The three wells in Hebron Township, Potter County, were on two different private properties. At each well site, measurements of 1 to 6 controls located 0.1 to 62 meters from the measured well were taken. Wells were selected to ensure that the numbers of plugged and unplugged wells measured were representative of abandoned wells in PA (Figure 4.1). The proportion 58 Table 4.1: Summary of sampling campaigns. Wells Controls July 17 - 18, 2013 1 0 Average Of Mean Daily Temperatures over the Sampling Round (◦ C) 24 July 24 - 25, 2013 5 3 13 Otto Township July 31 August 1, 2013 8 11 17 Otto Township October 8 - 11, 2013 13 23 11 Otto Township, City of Bradford, Allegheny National Forest January 26 - 31, 2014 14 15 -12 Otto Township, City of Bradford, Allegheny National Forest, Potter County Sampling Campaigns Number of Flux Measurements Locations Otto Township of measured wells that are plugged (26%) is similar to the proportion of wells that are plugged on the PA DEP list, which is 29%. Only 1 of the 19 measured wells was on the PA DEP’s list of abandoned and orphaned oil and gas wells as of January 17, 2014. 4.2.2 Flux chambers and sampling Chamber design and construction Static chamber methodology was adapted from techniques to measure trace gas fluxes from soil-plant systems [47, 75]. For well measurements, multiple-component chambers of various geometries were designed to enclose the wellhead and measure the emissions of methane and other trace gases from the well and the immediate sur59 Side View Plan View w w Well d h3 Well Groove Rigid Frame h 2 Collar Ground Surface h1 Chamber Figure 4.2: Schematic of flux chamber collar and frame enclosing surface protrusions of an abandoned oil and gas well. All variables other than the chamber diameter, w, depend on location. rounding area (Figure 4.2). For controls, both single- and multi-component chambers were designed. A schematic of the multiple-component chamber (Figure 4.2) shows the distances, d and h3 , that separate the well from the chamber. To accommodate different geometries of surface protrusions (e.g. well casing, monument/marker) and/or cement bases, we built three multiple-component chambers with the following footprints: circular with a 29 cm diameter (small), circular with 50 cm diameter (medium 1), circular with 55 cm (medium 2), and rectangular with widths of 90 cm (large). The chamber heights were adjusted to accommodate surface protrusions and ranged from 0.23 to 2 m. The resulting flux chamber geometries were cylindrical and rectangular prismatic. 97% of measurements were performed using cylindrical chambers. In the single-component chambers, vent tubes and sampling ports were installed through holes drilled in the rigid bucket. The multiple-component chamber has three major parts: collar, rigid frame, and bag. The chamber collars were made of prefabricated cylindrical buckets of rigid plastic material that were trimmed using a jigsaw or aluminum siding with grooves in which a tight bungee cord would be used to seal the collar and the bag (Figure 4.2). Rigid frames were built to support the bag and form cylindrical or rectangular chambers. The frames were made of PVC pipe, rigid plastic, aluminum siding, and/or fence posts. Bags were made of 4 mm60 thick polyethylene drop cloths and designed to snugly fit over the rigid frame and be airtight. Single-component chambers have a height of 25 cm and enclose volumes of 14 to 16 L depending on field conditions. Multiple-component chambers are larger and enclose volumes of 37 to 1400 L. At least one set of vent tubes and sampling ports were installed on each of the bags. Vent tubes made of Tygon lab tubing with 1/8” (3.1 mm) inner diameter were installed on the top face of each chamber to transmit atmospheric pressure changes to the enclosed air volume. Tube lengths were approximately 1 m, which is more than a factor of 4 greater than minimum lengths recommended in [47]. Sampling ports with one-way Luer connections were created on the side of the chamber at approximately half the height of the chamber. The vent tubes and sampling ports were installed on the chambers using bulkhead connectors, brass hose barbs, and Teflon tape. To ensure sufficient mixing in these larger chambers, we installed closed-air circulation systems composed of one or more 12 Volt DC cooling fans each powered by 8 AA batteries in all multiple-component chambers. The fans and battery packs were installed on the rigid frame to maximize mixing inside the chamber. Chamber deployment To obtain a snug fit, minimize soil disturbance, and limit air leakage, the singlecomponent chambers and multiple-component chamber collars were placed in 1” to 2” grooves (h1 in Figure 4.2) in the soil. For multiple-component chambers, the rigid frames were placed on top of the collars, which were secured by PVC fittings, and enclosed within the chamber bag. The chambers were sealed around the collar with either a water lock or with tight bungee cords that fit into grooves around the collars. Chambers enclosing wells were installed around surface protrusions and in some cases, cement bases, and did not touch any material visibly associated with the well. 61 We also chose flux chambers based on magnitudes of the methane flux. For high emitting wells, larger chambers were deployed to ensure that the chamber could be considered static and concentration changes could be measured at 5-minute intervals. For low emitting wells, smaller chambers were used to minimize sampling duration and maximize mixing. Sampling Gas samples from chambers were collected in 20 mL Wheaton serum vials for gas chromatographic analysis and 150 mL Wheaton glass serum bottles or 0.5 to 3 L SamplePro FlexFilm air sample bags for isotopic analysis. Butyl rubber stoppers were used in all Wheaton vials and sealed with aluminum crimps. Used vials were flushed with more than 100 mL of ultra-high purity nitrogen or air zero. All vials were flushed at least twice with ambient air or air zero and evacuated with a hand or mechanical pump to pressures of <10 kPa. We took 3 to 17 gas samples in 20 mL vials mainly at 5, 10, and 20-minute intervals over durations of 20 minutes to 25 hours depending on knowledge of the methane emissions. For each sample, a 50 mL air sample was extracted from the chamber using a 60 mL syringe, which was then transferred to a vial, using a needle, or directly to a bag. To minimize contamination, syringes were flushed at least 10 times prior to use in a new chamber and at least 5 times between samples taken at different times from the same chamber. We also took one 150 mL and one 0.5 L or 3 L sample from the chambers for isotopic measurements at the end of each measurement period. Throughout the sampling period, the pressure in the chamber was maintained at atmospheric pressure through a vent tube, which was tested occasionally with a flow meter. The samples in glass bottles were analyzed within 3 (4) weeks of collection for alkane (isotope) concentrations. 62 4.2.3 Analysis Alkane concentration C1 −C4 light hydrocarbons, including CH4 , C2 H6 , C3 H8 , and n-C4 H10 , were analyzed using flame ionization gas chromatography (GC) on a Shimadzu GC-2014 instrument. The carrier gas was ultra high purity helium, and hydrocarbon gases were separated on a 10-foot packed Porapak-Q column with a temperature program that involves a 2-minute isothermal period at 100◦ C followed by a temperature ramp of 10◦ C per minute to 150◦ C. The flame ionization detector was held at 200◦ C. Air samples were extracted from the vials into a glass syringe, equilibrated to atmospheric pressure, and injected into the instrument via a sample loop. Air samples were discarded when the volume extracted into the syringe was less than 5 mL. Instrument precision based on triplicate injections of a 1.00% methane standard is 2% for the July/ August 2013 samples, 8% for October 2013 samples, and 2% for January 2014. Peak identification in the sample chromatograph was accomplished by comparing retention times to those in the chromatograms of pre-mixed gas standards from AirLiquide. These standards were used to develop calibration curves for each gas. The calibration curves were obtained using linear regression of the peak areas and the mixing ratios of the pre-mixed standards, with the intercept set to zero. For methane, we used standards at 2.04 ppmv, 5 ppmv, 50 ppmv, 122 ppmv, 2527 ppmv, 1%, and 100%. Based on this calibration information, the methane concentrations observed in the chamber ranged from 0.82 ppmv to 52%. For C2 -C4 alkanes, we used a standard with ethane, propane, and n-butane concentrations of 5 ppmv, 50 ppmv, and 1%. The 5 and 50 ppmv standards also contained i-butane and were only available for the January 2014 samples. The objective for measuring the C2 -C4 alkanes was to determine their presence and their relative concentrations. Detection is defined qualitatively by the shape of the peak, relative to the baseline, and a minimum peak 63 area of ∼100 mVolts·min. Alkane presence at a well or control location was noted when alkanes were detected in at least two samples. For these locations, the average ratio of the C2 -C4 alkane with respect to methane was calculated. Isotopic concentration To measure the C isotopic composition of CH4 , a near-IR Continuous Wave-Cavity Ring-down Spectrometer (CW-CRDS) was used. This system consists of four distributed feedback laser diodes (DFB-LD), three of which were tuned to the absorption line peaks of 12 CH4 , 13 CH4 and CH3 D at 6047 cm−1 , 6049 cm−1 and 6108 cm-1, re- spectively, and a fourth that measured the baseline at 6051 cm−1 . The multiple laser design provides long-term stability of the system and increases the data acquisition rate. The acquisition frequency was further increased by utilizing a semiconductor optical amplifier to initiate cavity ring-down events. Optical isolators and spacing and orientation of optical elements were all used to prevent any Etalon effects. The optical cavity has a length of 0.653 meters with one pair of super mirrors with reflectivities of over 99.9993% at a wavelength of 1.651 µm. This is equivalent to a 93.3 km absorption path length. A heat controller stabilized the temperature of the cavity at 30±0.05◦ C. The high repetition rate combined with the super-high reflectivity mirrors and long-term temperature stability yields isotopic measurements of CH4 that are much more precise than can be achieved on current commercial CRDS instruments. The system has a detection limit of 1.9x10−12 cm−1 corresponding to 10 pptv of CH4 at 100 Torr. CH4 gas isotopic standards that have been analyzed by isotope ratio mass spectrometers were used to calibrate the CRDS. The δ 13 C-CH4 values were reported relative to Peedee belemnite (VPDB) standard. For ambient air samples that contained 1.9 ppmv CH4 the precision of the δ 13 C of the CH4 is ±1.7‰[13]. The 25 cc samples were equilibrated with the evacuated cavity to a pressure of ∼50 Torr to reduce peak overlaps between CO2 , H2 O and 64 12 CH4 absorption peaks with the CH3 D absorption peak. The resulting isotopic precision of δ 13 C for atmospheric CH4 was > 2‰ if the concentration is over 100 ppmv. 4.2.4 Flux calculation Fluxes, in units of mass per time per well, starting at the moment of chamber deployment, were calculated using linear regression in Matlab on the concentration data, c [mass/volume], over time, which is then multiplied by the chamber volume: F = dc · Ve , dt (4.1) where dc/dt is the slope of the fitted line for c(t) and Ve is the effective chamber volume. At least three data points, each representing different times, were used in the linear regressions. For most fluxes, five to seven data points were used. We assumed that the volumes of surface protrusions were small relative to the chamber volume. For control locations, fluxes were scaled based on the land areas covered by the chamber for the control and the nearest well location as follows: Fcontrol,scaled = Fcontrol,raw · Awell Acontrol (4.2) where Fcontrol,scaled is the methane flux scaled to the area of the well [mass/time/well], Fcontrol,raw is the methane flux of the control before scaling [mass/time/control], Awell is the ground surface area of the chamber at the well [length2 ], and Acontrol is the ground surface area of the chamber at the control [length2 ]. We calculated fluxes of methane only. Fluxes of C2 -C4 alkanes are not presented here since there are insufficient numbers of ethane, propane, and n-butane detections to calculate fluxes in most locations. Also, we have low confidence in the C2 -C4 alkane concentrations since a standard in the range of most samples was not available at the time of the July, August, and October 2013 sampling campaigns. 65 In several chambers, changes in concentration declined with time during the sampling duration. This decline in methane accumulation could be because the chamber could no longer be considered static or due to altered diffusion gradients. For these flux measurements, data from the earlier time points where the relationship is linear and more representative of the flux at the time of chamber deployment, and less likely to be affected by altered concentration gradients, were used. Data from the first 10 to 81 minutes were used for 88% of the fluxes. The sampling duration was based on the assumed magnitude of the flux, with longer durations used for smaller fluxes. The goodness of fit obtained by linear regression was evaluated with the R2 value. Here, R2 values greater than 0.8 were assumed to lead to a good flux estimate. However, fluxes with low R2 values were not discarded here since this would lead to biased results that favor higher fluxes [76]. R2 values corresponding to the measured methane fluxes are correlated with the absolute value of the flux with low R2 values associated with smaller fluxes (Figure 4.3). Therefore, p-values (null hypothesis: flux = 0) were calculated and fluxes were set to zero for p-values greater than 0.2. This resulted in 2 of the measurements at wells and 12 measurements at controls to have zero fluxes. The 2 zero-flux measurements at wells are one of three measurements made at each of the two wells, which were sampled over 3 sampling campaigns; the other two flux measurements at these two locations are positive and non-zero. Flux values were discarded if there was clear evidence of contamination. Evidence for contamination is determined by considering the sequence of sampling and equipment usage, concentrations observed at the last location of equipment usage, and comparisons of fluxes to other measurements and published values. A total of 3 flux measurements out of a total of 97 flux measurements were discarded. Fluxes from controls were compared to published values. Of the measured fluxes at well locations, 90% are 1 to 8 orders of magnitude greater than methane fluxes observed in other areas with similar land cover but without wells [54, 2, 48]. In 66 (" .$" !#'" !#&" !#%" @066A" BC41=C6A" !#$" !" (#)*!+" (!+" (!," (#)-!+" (!*+" (#)*!," (!*," (#)*!(" (!*(" (#)-!(" (!(" (#)-!," /012340"5678"9:;<2=<>066?" Figure 4.3: The R2 value for the linear regression and the corresponding methane fluxes. contrast, over half of the methane fluxes from control locations are representative of fluxes found in other areas with similar land cover. However, it is important to note that the control locations may also be affected by the presence of the well. 4.3 4.3.1 Results Methane fluxes Methane fluxes from abandoned wells were found to be significantly higher than methane fluxes observed at controls (Figure 4.4). The mean flux at well locations was 11,000 mg/hr/well (0.27 kg/day/well), while the mean flux at control locations was 0.19 mg/hr/location (4.5x10−6 kg/day/location). Methane emissions at well locations were based on good linear fits with 88% of the fluxes having R2 values greater than 0.8. Positive methane fluxes were observed at all 19 wells with average values ranging from 6.3x10−1 mg/hr/well to 8.6x104 mg/hr/well. Average methane fluxes at controls ranged from -1.2x10−1 mg/hr/location to 4.2 mg/hr/location. Sources of uncertainty included flux chamber design, deployment, sampling, laboratory analysis of samples, and data selection for regression analysis. We estimate that the combined effect of the various sources of uncertainties in flux estimates will lead to errors within a factor of 2 67 !"#$%&$%'()*+,(%!-./(0%% *+,-.+/0' 95->:,5'J/3M50'' N?4G>.GO5//'+.'?4G>.G/+P:Q+,R' !"#(%)' !%)' J+.50-' 1,2/34456'75//0' 8/34456'75//0' F.:00/:,6' K@L5.' 75-/:,6' !%&' !"#(%&' !"#(%!' !%!' !"#$%!' !%$!' !% !"#$%&' !"#$%&' $!%$&' $&' $!%$!' !"#$%!' Figure 4.4: A total of 42 and 52 measurements were made at wells and at locations near the wells (controls), respectively, in forested, wetland, grassland, and river areas !"#$%&'#()&%*)+,#)+$#")#-)./#0$%1./$0$+"1#2$/$#0%3$4##!((#2$((#()&%*)+1#'%5$#%"#($%1"#)+$# in July, August, October 2013 and January 2014. &)+"/)(#0$%1./$0$+"#26"'6+#78#0#)-#"'$#2$((#()&%*)+4### )G!G!H' 9:.;'<:,4='95->:,5'#?@00@+,0'A.+?'BC:,6+,56'DEF'75//0' I' of our estimate. This error is small relative to the seven orders of magnitude variation in measured fluxes. Furthermore, most of the sources of measurement uncertainty would bias the measured fluxes to be lower than their actual value. Uncertainties in methane flux measurements are discussed further in Section 4.4. Methane fluxes at well locations appeared to be unaffected by land cover, which included forest, grassland, river, and wetland. In contrast, we found that methane fluxes at control locations were dependent on land cover. A large proportion of fluxes from controls in forests and grasslands were negative (i.e. methane sinks) and ranged from -10−2 to 10−1 mg/hr/location, while the fluxes from controls in wetlands were consistently positive and relatively high, ranging from 10−2 to 101 mg/hr/location. We found seasonal effects were present in controls, with lower methane fluxes observed in the January 2014 sampling round. While there is no evidence of significant seasonal effects in the methane fluxes from wells, additional measurements are needed to reach a firm conclusion. According to regulations on well abandonment, wells are plugged to limit vertical migration from subsurface source formations (oil and gas reservoirs and coalbeds), 68 which includes minimizing impacts on groundwater. We found that methane fluxes from plugged wells were not necessarily lower than methane fluxes at unplugged wells. For example, in the grassland area, both the largest and the second lowest methane fluxes originated from plugged wells. Evaluation of plugging status and wellbore integrity was difficult using only visual inspection at the surface and lack of additional information. 4.3.2 Presence of ethane, propane, and n-butane The presence and concentration of ethane, propane, and n-butane are useful for identifying the methane source as thermogenic or microbial. Because ethane is not coproduced during microbial methanogenesis, the presence of ethane-to-methane ratios greater than 0.01 indicates gas of largely thermogenic origin in groundwater ([53, 81]). A similar threshold is not readily available in literature for propane-to-methane and n-butane-to-methane ratios, but we expect this threshold value to be less than 0.01. Ratios of ethane, propane, and n-butane relative to methane were more frequently greater than 0.01, and at higher values, for wells than for controls (Figure 4.5). The presence of these non-methane hydrocarbons in controls indicates that there may be subsurface horizontal gas flow away from the well and subsequent emissions to the atmosphere. We did not find a consistent ratio for wells or controls with the alkane ratios ranging from 1x10−5 to 0.8. The high variability in alkane ratios may be a result of mixing between various microbial and thermogenic (deeper) sources. 4.3.3 Carbon isotopes of methane Carbon isotope information provides additional evidence suggesting that the source of methane from the wells is likely to represent a mixture of microbial and thermogenic sources. In general, methane originating from thermogenic sources is more enriched in 13 C, whereas methane originating from microbial sources is relatively depleted in 69 =#B*C07,D#EF()# 478*+,9:,()*+,#;*!=# &!=# %!=# !=# H/+(./7D# I,77D# '()*+,# -./0*+,# +123(*+,# '()*+,# -./0*+,# +123(*+,# J*K# JLK# Figure 4.5: Average alkane ratios ([C2 H6 ]/[CH4 ], [C3 H8 ]/[CH4 ], and [n-C4 H10 ]/[CH4 ]) (a) and proportion of samples with alkane ratios greater than 0.01 (b) at control and well location with detectable ethane, propane, and n-butane concentrations are calculated for samples collected in July, August, and October 2013 and January 2014. 13 C. We found that the samples collected at wells were likely to be more enriched in 13 C than those collected at controls (Figure 4.6). A comparison of the methane δ 13 C values to that of known thermogenic and microbial sources [77] indicates that most of the methane fluxes from wells were thermogenic or a mixture of microbial and thermogenic sources. Three of the 26 measurements at wells had methane δ 13 C values in the microbial range. The methane δ 13 C values of the measured wells ranged from -71‰to -21‰. This range is broader than published methane δ 13 C values of thermogenic methane in natural gas in the northern Appalachian basin, which range from -47.9‰to -30.7‰[41]. The methane δ 13 C values at controls ranged from -85‰to -56‰, indicating control sources were more likely to be of primarily microbial origin. Figure 4.6 also shows that locations with larger methane fluxes emitted methane that was more enriched in 13 C. Wells with methane fluxes that were greater than 103 mg/hr/well were likely to be emitting methane of thermogenic origin; while wells with fluxes in the order of 100 to 101 mg/hr/well emitted methane of microbial, thermogenic, or mixed thermogenic/microbial origin. Methane emitted from most control locations is in the microbial range; however, one measurement reveals that methane emitted from control sources can contain thermogenic sources of methane as well. If we consider the integrated fluxes from all these wells, the methane emitted is 70 FDE<9IDJ=' !%+' !"#7%+' FDG?H' @A?9*01) !-(' %"()' !%$!' !"#$%!' !%$&' !"#$%&' !"#$%&' $!%$&' $!%$!' !"#$%!' %"&*' ,R' %"!!' !"!' %")!' $(%' $)%' $*%' $+%' $,%' $-%' $&%' !!-.$./,'01'23456' FDE<9IDJ=' %"%%!' ".&78./,'' !%%' .9:;<9=' >?==' FDG?H' @A? -&#+,?'@' A$2#")B$&=' 69C6 C8=96C F-'GHI' F&+/&$JJ' I$*+&)J='6986: 699C !"#$%&' (#)$&*+,%)"+# 69;9:69D9 698>:6996 D8=69; 699;:699D 6998:;>>7 ;>>9:;>6C HJ)"2%)$?'T32K$&'+N'O$,,J' T+)'-PP+3#)"#/'N+&'H#S%#P$?' G&",,$?'"#'M"2$'F$&"+?J'U")S' I$P+.$&L'O$,,J 1"JJ"#/'G%)% ;E=DC6 -.$&%/$ CD=;>E -PP+3#)"#/'N+&'H#S%#P$?' I$P+.$&L'O$,,J 01"#"2324 01%5"2324 01"#"2324 01%5"2324 6<7=69> 6<7=69> <>8=D7D <>8=D7D CD=;>E D8=69; 6;C=6DD 6<;=<9; C8=96C C8=96C 6;9=;7E 6;9=;7E D=C86 8=6C; D=C86 8=6C; D=C8D 67=DED ;7=;79 ;C=9D7 ;C=9D7 ;7=;79 ;C=9D7 ;7=;79 ;E=DC6 ;E=DC6 ;E=DC6 ;E=DC6 ;9D=>C9 C6>=6DE 96C=6D4 !"#$%&'()& *+%%, 6>=9;6 -./0/// 12/0/// average of wells drilled per year from the preceding and following time periods and interpolating between these two numbers. Both methods are used to determine the number of wells drilled in the years with missing data. With these interpolations, we can use the data in Figure 5.2 to calculate the total number of wells in Pennsylvania during these time periods. From this total number, we subtract the number of active and inactive wells (available from the DEP), which is 10,921 as for March 1, 2014, to get the total number of AOG wells (Table 5.2). The total number of AOG wells, not accounting for enhanced recovery wells, range from 280,000 to 310,000. After the Bradford Oil Field peaked in 1881, the field continued to produce significant quantities of oil due to enhanced recovery methods, mainly water flooding. Ref. [3] estimate the discovery and development of water flooding techniques to be in the 1880s and 1890s. However, water flooding was illegal in Pennsylvania until 1921 86 and water flooding wells drilled prior to 1921 are unlikely to be represented in any records. Even after 1921, it does not appear that enhanced recovery wells were considered oil and gas wells and are likely to be under-reported until 1995 when the 1992 amendments of the Oil and Gas Act of 1984 came into play. In some DER progress reports, the number of oil and gas wells presented are stated to not include wells for enhanced recovery; while others are less clear. One exception is the progress report for 1952 in which [22] states that the number of oil and gas wells (including dry wells) drilled in 1952 is 1899, of which 1366 (72%) are for enhanced recovery. If we use the 1952 ratio of all wells to wells not related to enhanced recovery (1899:533 = 3.6:1) from 1859 to 1994, we obtain 910,000 to 970,000 AOG wells (Table 5.2). There is significant uncertainty in the numbers for 1859 to 1928 and thus, it is difficult to say that application of the factor of 3.6 leads to an overestimate. In the 1929-1949 time period, the factor of 3.6 is likely to be on the low side for two reasons: (1) the “fivespot” method, developed in 1928 [3], meant that 4 additional wells were drilled for enhanced recovery, and (2) we did not represent the increase in drilling that is likely to have accompanied Pennsylvania’s second peak in production in the 1930s and 1940s, which largely came from water flooding in the Bradford Oil Field. For the 1950-1994 time period, the factor of 3.6 may lead to an overestimate near the end of this time period as water flooding became less popular. Therefore, the factor of 3.6 applied to the 1849 to 1994 numbers may lead to both over and underestimates. For the 1995 to 2007 time period, our estimates are likely to be on the low side since we did not apply any factors for enhanced recovery and our interpolation methods are unlikely to have captured any peaks in wells drilled prior to 2009 (although a peak in the number of permits issued occurred in 2008). Putting all of this together, we estimate the number of AOG wells in Pennsylvania to range from 280,000 to 970,000. This increases the range in estimates of methane emissions from AOG wells to 4 to 13% 87 of current total estimated annual anthropogenic methane emissions in Pennsylvania using the average methane fluxes of the 19 measured wells. Figure 5.2 also shows the number of permits issued for oil and gas wells (DEP database). There are 18 permits issued from 1918 to 1955. It is unclear whether these permits are a result of typographic errors or are actual permits issued. Either way, the number of permits issued from 1918 to 1955 is not useful for determining the number of oil and gas wells drilled during this time period. For example, other sources report number of wells from the 1918 to 1929 to be in the thousands, however, virtually no permits were issued during this time period. In 1956, there is a dramatic increase in the number of well permits issued, which is when permitting of new drilling in Pennsylvania began. However, the number of permits are less than the total number of wells drilled including those for enhanced recovery, confirming that permitting of enhanced recovery wells was not requied. In general, trends in the number of permits issued are echoed in the number of wells drilled but with less pronounced peaks, which occur in 1984 and 2008. Pennsylvania’s Oil and Gas Act of 1984 modernized Pennsylvania’s regulations related to oil and gas development and some of the permits may be for wells drilled prior to 1984. The second peak in 2008 is likely related to the recent shale gas boom in Pennsylvania. The proportion of permits issued that end up as new wells drilled varies from 27% to 85% with an average of 49% for 1951 to 2013. Estimating the number of wells drilled based on permits issued is likely to lead to large underestimation or overestimation and thus, we do not directly use the number of permits issued to estimate the number of wells drilled here. So far, we only focused on AOG wells, which includes orphaned and lost oil and gas wells. However, there is another well category, in which wells are labeled as “inactive”. The numbers presented above do not include inactive wells that may also be leaking methane into the atmosphere. According to Chapter 78 of the Pennsylvania Code 78.102 on “Criteria for approval of inactive status”, the owner or operator must 88 monitor these wells on an annual basis, specifically the “annular vents for gas flow volumes”. If this gas flow volume exceeds 5,000 cubic feet per day (141.6 m3 /day), “the owner or operator is required to notify” the DEP and “take remedial action approved” by the DEP. This rate is significantly larger than the highest measured methane flux in Chapter 4, which is 3.2 cubic meters per day. Since we do not have any data on these “inactive” wells, we assume that they are sufficiently managed and are not included in the above calculations but note that they may also be a significant source of methane to the atmosphere. Given the 150-year history of oil and gas development in the U.S. and the corresponding changes in technology and regulatory environments, the number of AOG wells in other states is likely also uncertain and may require analysis similar to that performed here for the number of AOG wells. Figure 5.2 shows the number of wells drilled in California from 1973 to 2009 are similar to the number drilled in Pennsylvania. California, where the first commercial oil well was drilled in 1865 [6], is likely to have a number of AOG wells on the same order of magnitude as Pennsylvania. 5.2.2 Emission estimates per AOG well Measurements in Chapter 4 show that the average methane flux from AOG wells are controlled by outliers. In fact, removing the three AOG wells with fluxes of 104 mg/hr/well results in an average that is an order of magnitude lower. This is significant since the average flux is used as the “emission factor” when estimating methane emissions from AOG wells. Therefore, the distribution of these outliers and the attributes that identify them are important. Using data from Alberta’s Energy Resources Conservation Board, [87] identified and evaluated well attributes that impact well bore leakage, which are summarized in Table 5.3. This work was done in the context of understanding the role of abandoned wells and the risk of leakage from CO2 storage formations in deep saline aquifers. 89 Information on most of the attributes listed in Table 5.3 is not available for the AOG wells measured and presented in Chapter 4. As a result, we focus on geographic area and whether the AOG well is plugged (i.e. cemented). However, as more measurements and information becomes available, other attributes identified by [87], including those found to have minor or no impact, should be explored. Table 5.3: Well attributes and their impact on well bore leakage by [87] No  Impact   •  Well  age   •  Well-­‐opera4onal   mode   •  Comple4on   interval   •  H2S  or  CO2   presence   Minor  Impact   •  Licensee   •  Surface-­‐casing   depth   •  Total  depth   •  Well  density   •  Topography     Major  Impact   Geographic  area   Wellbore  devia4on   Well  type   Abandonment  method   Oil  price,  regulatory   changes,  and  SCVF/GM   tes4ng   •  Uncemented  casing/hole   annulus   •  •  •  •  •  Table 5.4 presents the average, maximum, and minimum fluxes based on well status (plugged vs. unplugged) and geographic area (county) using the measurements from the 19 wells presented in Chapter 4. Unplugged wells have a mean methane flux that is 2 orders of magnitude larger than the mean for plugged wells. Fluxes from unplugged wells range over 6 orders or magnitude; while fluxes from plugged wells range over 4 orders of magnitude. Both the largest and smallest fluxes are observed at unplugged wells. Further, the maximum methane flux from plugged wells is only 1 order of magnitude smaller than the maximum methane flux from unplugged wells. The narrower range of values for plugged wells may be because fewer plugged wells were measured. Therefore, it is difficult to conclude based on these measurements, especially given the relatively small number of samples, that plugged wells emit less methane than unplugged wells on a per well basis. In terms of geographic area, Potter County, where gas development in the Oriskany sandstone took place, appears to have higher methane fluxes per well. However, relatively few measurements were 90 Table 5.4: Measured methane fluxes by attributes Attribute Well@status Unplugged Plugged Geographic@ Potter@County area McKean@County Well@status Methane@Flux@per@Well@(mg/hr/well) Average Minimum Maximum 2E+04 6E801 9E+04 8E+02 1E+00 4E+03 6E+04 6E+03 9E+04 3E+03 6E801 5E+04 Count 14 5 3 16 made in Potter County (3 out of 19) and these measurements were confined to a relatively small region (<2 mi2 ) in the northeastern portion of the county. Additional measurements of methane fluxes from AOG wells should be made and the attributes that most influence these fluxes should be identified and evaluated. Understanding the attributes that control emissions from AOG wells can allow for the development of efficient measurement plans and mitigation strategies. 5.2.3 Reporting and monitoring To reduce uncertainties in the methane emissions potential of AOG wells and to promote mitigation, improved reporting and monitoring of AOG wells is needed. In Pennsylvania, oil and gas development and production is overseen by the DEP’s Office of Oil and Gas Management (OOGM) pursuant to the Oil and Gas Act, the Coal and Gas Resource Coordination Act, and the Oil and Gas Conservation Law. This office is also responsible for the Abandoned and Orphan Well Program. A well owner or operator is supposed to identify AOG wells on property that they own or lease (The General Assembly of Pennsylvania, House Bill No. 1950); however, there is no clear requirement for these wells to be reported to the DEP. According to the DEP website, if an AOG well is found by the general public, the DEP is supposed to be notified within 60 days of its discovery. Once an AOG wells is reported, the DEP is supposed to investigate the well to determine if it qualifies as an abandoned or orphan well and should include the well in its database. Given the fact that only 12,000 wells out of the 280,000 to 970,000 existing AOG wells are included in this database, there 91 appears to be a lack of discovery and/or reporting of AOG wells and updating of the database. In fact, DEP’s database of AOG wells grew by only 66 wells in a 10 month period (Table 5.1). At this rate, it will take over 3000 years for even the lower bound in the number of AOG wells (280,000) to be reported and included in DEP’s database. Reporting and monitoring under the current Abandoned and Orphan Well Program in Pennsylvania is not designed to quantify methane emissions. Methane emissions from AOG wells are not included in any emissions inventories. However, their inclusion in emissions inventories is likely to be important for promoting mitigation since inventories are used to evaluate emissions reduction priorities [64]. GHG emissions including emissions of carbon dioxide, methane, nitrous oxide, hydrofluorocarbons, perfluorocarbons, and sulfur hexafluoride are estimated by the Pennsylvania Climate Change Advisory Committee (CCAC) as a part of the Commonwealth’s Climate Change Action Plan (CCAP), which must be updated every three years pursuant to the Pennsylvania Climate Change Act of 2008 [64]. To estimate the GHG emissions for the Commonwealth, the CCAC uses U.S. Environmental Protection Agency’s (USEPA’s) State Inventory and Projection Tool [64] and presents the results in the CCAP in terms of CO2 equivalents, lumping together the GHGs emitted from each industrial sector. Much of the CCAP is dedicated to detailing plans for reducing GHG emissions, which include switching from coal to natural gas for electricity production, reduced fugitive emissions, and inclusion of waste-to-energy facilities in the AEPS. There is no mention of AOG wells in the latest CCAP [64]. Methane emissions are generally estimated as a part of GHG emissions inventories by the USEPA at the federal level and by states using USEPA’s guidelines. Therefore, for methane emissions from AOG wells to be included in GHG emission inventories both at the state and federal level, action must be taken by the USEPA. THe USEPA should develop and adopt a methodology to account for these emissions and guidelines 92 should be published. Current methodology should focus on measurement-based methods so that more data on methane emissions potential become available. Given the fact that methane emissions from AOG wells are spatially distributed, measurements of these wells are unlikely to be promoted in USEPA’s GHG Reporting Program, which are for facilities with emissions equal to or greater than 25,000 metric tons of CO2 e (40 CFR Part 98). It is difficult to implement a reporting program that relies on industry to report methane emissions from AOG wells since there is no industry that remains responsible for them. Therefore, an alternate approach to incentivize and finance reporting and monitoring is required. Further, state programs responsible for AOG wells should not focus solely on plugging but also on building an accurate and comprehensive database of these wells. Current Pennsylvania regulations (Pennsylvania Code, Chapter 78) require well records providing many of the attributes listed in Table 5.3 to be provided to the DEP. However, this regulation does not address the issue of missing data on older AOG wells, which require compilation of historical public and industry records. 5.3 Mitigation options Regulations for oil and gas wells in Pennsylvania, but also in most states, require wells to be plugged upon abandonment for groundwater protection, oil and gas resource conservation, and safety. Reducing gases emitted to the atmosphere is not one of the goals of plugging and in the coal areas of Pennsylvania, venting may even be required. Therefore, we review current trends in plugging of AOG wells and regulations on plugging to evaluate the potential of current regulations to reduce methane emissions from AOG wells. Given that the total methane emissions from AOG wells are governed by the high emitting wells, we explore the option of using gases emitted from these wells. Specifically, we review Pennsylvania’s AEPS to determine if gases 93 emitted from AOG wells can be considered an alternative energy source under current regulations. 5.3.1 Plugging There are many methods to plug a well but in general, a plug (often a cement plug) is placed and cemented in an abandoned well over intervals that intersect formations to be isolated. Formations that require isolation include oil and gas formations (including those that are depleted), coal layers, and groundwater aquifers. To ensure that these formations are sufficiently isolated, plug intervals are required to extend from a specified distance below the bottom of the formation to a specified distance above the top of the formation. Intervals below or between plugs may be left unfilled or may be filled. Additionally, a plug is placed and cemented near the surface and the well is filled with material to the surface. Often, a marker at the surface is required to identify the well as a plugged well. Also, production casings are typically removed before plugging. According to the Interstate Oil & Gas Compact Commission (IOGCC), the average cost to plug a well varies from $6,000 to $12,000 [37]. In Pennsylvania, the DEP is responsible for plugging orphan AOG wells, which is funded by the Abandoned Well Plugging Fund and the Orphan Well Plugging Fund. A $50 surcharge is added to new well permits for the Abandoned Well Plugging Fund. Additional surcharges of $200 for gas wells and $100 for oil wells on new well permits are used to finance the Orphan Well Plugging Fund. A total for 66 wells were plugged from April 2013 to January 2014 according to changes in the number of DEP plugged wells on the DEP’s list (Table 5.1). At this rate, it will take 180 years to plug the 6,300 orphan wells currently on the DEP’s list. (We do not include abandoned wells that are not considered to be orphans here since the costs of plugging should be covered by the owner / operator of the well.) This implies that there may not be sufficient funding available for plugging and that changes to the Orphan Well 94 Plugging Fund may be required. Currently, given the large number of AOG wells that require plugging coupled with limited funding, the DEP prioritizes wells based on the risk they pose to human health and safety and their potential to contaminate air, water, and soil [63]. It is unclear whether these priorities also promote plugging of AOG wells that emit large quantities of methane. Table 5.5 summarizes key aspects of Pennsylvania’s regulations on plugging, which vary depending on whether the well is in a coal area, whether the well penetrates a workable coal seam, the condition of the surface (or coal protective) casing, and the condition of the production casing. In coal areas, vents are generally required, which implies that methane may be emitted intentionally from these plugged wells. In noncoal areas, vents are not required but the length of the plug below the surface is reduced from 100 to 200 feet to a 50-foot minimum. The plug material required for the plug interval over stratum bearing or having borne oil, gas, or water is generally cement, which is defined as a “mixture of materials for bonding or sealing that attains a 7-day maximum permeability of 0.01 millidarcies”. However, when the production casing is cemented, the regulations offer the option of a mechanical plug, and in workable coal seams, expanding cement plugs and bentonite gels are also given as options. The top and bottom of a plug interval is defined by specified distances above the top and below the bottom of the “stratum bearing or having borne oil, gas, or water” respectively. These specified distances varies from 20 to 100 feet, with less stringent options being available in the case of cemented production casings. The interval between plugs may be left unfilled when the production casing is cemented; these intervals are required to be filled with nonporous material when the production casing is not cemented or not present. Nonporous material is defined in the Pennsylvania Code, Chapter 78, as “nontoxic earthen mud, drill cuttings, fire clay, gel, cement or equivalent materials approved by the Department that will equally retard the movement of fluids”. However, in the Coal and Gas Resources Coordinate 95 Act, the definition of nonporous material is “sand pumpings, mud or other equally nonporous material” in section 13(1) but “bentonite gel or other equally nonporous material” in section 13(3). The porosity and permeability of sand and mud differ by many orders of magnitude [70] and therefore, it is unclear how much fluid migration between the plugs is permitted. Nonetheless, these plugging requirements should be sufficient for zonal isolation of oil and gas bearing formations (or stratums) and groundwater reservoirs, which ideally should prevent vertical migration of fluids from oil and gas formations to groundwater reservoirs or to the surface. However, a key assumption here is that the depths and thicknesses of the oil and gas formations and groundwater reservoirs are known and any uncertainties are within the 20, 50 or 100 feet buffers in plug intervals. Abandonment methods and regulations along with oil and gas production technology have changed significantly over the past 150 years of oil and gas production in North America [87, 43]. Many plugged wells are likely to have been plugged differently than as detailed in the regulations in Table 5.5. Furthermore, integrity of the plug and well bore degrades over time due to geologic conditions and events, other oil and gas production activities, and presence and chemistry of fluids that are in contact with the plug and well bore. However, there are no regulations that require monitoring of plugged wells for leakage of methane and other gases. 5.3.2 Usage as an alternative energy source? According to Pennsylvania’s Plan for Addressing Problem Abandoned Wells and Orphaned Wells [63], the DEP’s plan is to prioritize plugging “problem wells” and to encourage adoption. (Problem wells are defined as “wells that threaten the health and safety of people or property or pollution of the waters of the Commonwealth”.) Encouraging adoption is aimed at the oil and gas industry. AOG wells are unlikely to produce oil and/or gas above the economic limit. However, the industry may adopt 96 a well to use it as an enhanced recovery well or to use it for other purposes such as waste injection. Given that the number of AOG wells on the DEP’s list have grown from approximately 8,000 in 2000, when the plan was published, to 12,000 today, it appears that additional strategies should be added to Pennsylvania’s plan for addressing AOG wells. We propose the strategy of considering gases emitted from AOG wells to promote its use locally by landowners for residential or commercial use and reporting and monitoring of these wells. According to the U.S. Energy Information Administration (EIA), the residential, commercial, and industrial price of natural gas as of January 2014 are $9.24, $8.09, and $5.61 per thousand cubic feet of gas. The highest emitting abandoned well emits 3.2 cubic meters of methane per day. We can use the rate of emitted methane only (although natural gas is generally 70-90% methane) to estimate the value of the emissions at $376, $330, and $220 per year for residential, commercial, and industrial users of this high emitting well. If plugging costs range from $6000 to $12,000, it would take 16 to 32 years to generate enough savings from the use of emitted methane, without accounting for interest rates and changes in the price of natural gas. Therefore, the savings from the use of the emitted methane may not be a sufficient incentive to promote use. Here, we consider the inclusion of gases emitted from AOG wells in Pennsylvania’s AEPS so that an additional source of revenue can incentivize their use. Alternative energy resources currently included in Pennsylvania’s AEPS are presented in Table 5.6. The alternative energy resources are as presented in the Pennsylvania AEPS Alternative Energy Credit Program website, the AEPS annual reports [72, 73], Pennsylvania Code Chapter 75.1, and in PJM-EIS generation attribute tracking system (GATS) (Table 5.7). We have attempted to match the definitions from these sources; however, there are some uncertainties as to whether we have correctly matched the categories with italicized defintions in Table 5.6. For example, black liquor is listed 97 as a resource type on the AEPS website and in GATS as a Tier I energy source. However, there is no explicit mention of black liquor in Chapter 75.1 and the closest possible definition is assumed to be biomass energy. Further, black liquor can qualify as both a Tier I and Tier II resource. Another arbitrary resource type is “Other gas”, which can be related to biomass energy, biologically derived methane gas, or demandside management. It appears that the actual alternative energy resource types are based on interpretations of the definitions in Pennsylvania Code Chapter 75.1 Therefore, we can consider gases, including methane but also other hydrocarbons, emitted from AOG wells as: • Biologically derived methane gas, “which includes methane from the anaerobic digestion of organic materials from yard waste, such as grass clippings and leaves, food waste, animal waste and sewage sludge. The term also includes landfill methane gas.” (Tier I) • Demand-side management, which includes “reuse of energy from exhaust gas or other manufacturing by-products that are used in the direct production of electricity at the facility of a customer” (Tier II) • Distributed generation systems, “ which means the small-scale power generation of electricity and useful thermal energy” (Tier II) It is important to note that many of the definitions use the word, “includes”, which allows for other sources not explicitly mentioned in the regulations to be eligible. The AEPS categorizes sources as Tier I and Tier II. The alternative energy requirement for 2012 was 3.5% and 6.2% of total energy sold for Tier I and II resources respectively [73]. There is a separate requirement just for solar, which is 0.0325% of total energy sold [73]. The weighted average price of alternative energy credits (AECs) in 2012 is $180.39, $5.23, and $0.17 for solar, Tier I, and Tier II resources 98 [73]. For the 2013 energy year, the AEC prices have increased for Tier I and II sources to $8.31, and $0.22. Each AEC equals 1 megawatt hour (MWh) of alternative energy produced. According the U.S. EIA, 7.96 Mcf (thousand cubic feet) of natural gas is required to generate 1 MWh of electricity, assuming a heat rate of 8,039 Btu/kWh and a fuel heat content of 1,023,000 Btu per Mcf. This translates to 5.1 AECs per year for the high emitting well (methane emission rate of 3.2 cubic meters per day), which amounts to $26.76 and $0.87 using Tier I and II average prices for 2012 and $42.52 and $1.13 using Tier I and II average prices for 2013. With an additional revenue of $42.52 assuming qualification as a Tier I source, the years needed to recover the cost of plugging range from 14 to 29 (without accounting for interest rates and changes in natural gas prices). Two of the potential alternative energy resource types that AOG wells can fall under are in Tier II. Given the low price of Tier II AECs, it is unclear whether the inclusion of AOG wells as a Tier II resource will be sufficient to promote its use. However, the AEC prices are highly volatile with Tier I prices ranging from $0.13 to $100.00 and Tier II prices ranging from $0.01 to $20.00 in 2013. Additional details on AEC prices, in terms of location, resource, and facility, are needed to evaluate the ability of the AEPS to promote the use of gases emitted from AOG wells. Further, the value of the emitted methane is quantified here relative to natural gas prices. Therefore, increases in natural gas prices can make the use of gases emitted from AOG wells more valuable. Note that other valuable gases (such as propane and butane) are also emitted from these wells and their value should also be included in any economic valuation. The potential for AOG wells to be considered as an alternative energy resource may promote reporting of these wells and monitoring of gases emitted. The reporting and monitoring promoted is likely to be biased towards high-emitting wells, which 99 will have a significant impact on overall methane emissions. A more detailed study on the economics of plugging and use, including interest rates and changes in natural gas prices, with detailed AEC prices is needed. 5.4 Conclusions With currently available data, we estimate that there are between 280,000 and 970,000 AOG wells in Pennsylvania, which translates to 4 to 13% of total estimated state-wide anthropogenic methane emissions in Pennsylvania. Methane emissions from AOG wells remain highly uncertain due to uncertainties in the number of wells drilled over the years and variability in methane emissions per well. The history of oil and gas development in the region is critical for determining the number of AOG wells. In Pennsylvania, this requires the consideration of wells drilled for enhanced recovery and the changes in regulations on permitting, registration, and reporting. In terms of emissions per well, high emitters govern total methane emission estimates and mitigation of only these high emitters (which represent 16% of measured wells in Chapter 4) leads to an order of magnitude reduction in estimated methane emissions. More field measurements are needed and measurement plans should be aimed at identifying attributes that aid in finding these high emitters. Reducing methane emissions from AOG wells does not appear to be a goal in well abandonment regulations. Well abandonment, as currently described in Pennsylvania’s regulations, focuses on well plugging and does not include monitoring. Measurements presented in Chapter 4 have shown that these plugged wells do emit methane. Although plugging should prevent vertical migration of fluids, especially from oil and gas formations uncertainties in the depths and thicknesses of formations bearing oil, gas, and/or water and the possible degradation of the plug over time makes the effectiveness of plugging for reducing methane emissions uncertain. Also, 100 in the case of wells in coal areas, gas vents designed to emit methane and other gases are required. Gases, including methane and other hydrocarbons, emitted from AOG wells may be considered as an alternative energy resource under current definitions. We believe that the consideration of gases emitted from AOG wells as an alternative energy source is valuable for both promoting its use but also for catalyzing reporting and monitoring. Based on current AEC pricing, the inclusion of gases from AOG wells will reduce the number of years to recoup the cost of plugging by 2 to 3 years. However, given the volatility in AEC and natural gas prices, further research on the economics of this option must be studied. Similar analysis of the history of oil and gas development is needed in other states to estimate the methane emissions from AOG wells for the state and for the U.S. Such analysis is necessary for including methane emissions from AOG wells into national GHG emissions inventories and to bridge the current gap in top-down and bottomup methane emission estimates. Improved methane emission estimates can aid in evaluating, developing and implementing cost-effective methane emission reduction strategies, such as addressing methane emissions from AOG wells. 101 Table 5.5: Summary of plugging requirements in Pennsylvania Code, Chapter 78, and Section 13 of the Coal and Gas Resource Coordination Act !"#$%&'()*+'" !"#$$+,$%J)N#+*)&),(/ !"#$%+,()*-'"%.-)*%/(*'(#&% 4,()*-'"%0)".3% 0)'*+,$%.*%1'-+,$%0.*,)%.+"2% ',5%0)(3)),% $'/%.*%3'()* 6"#$/ !"#$%0)".3% 9'/% (1)%/#*7'8)% -),( .*%8.'"% 6*.()8(+-)% 8'/+,$%/)'( !"#$%&:"()"+*,-",./,0"1234"+/5/)3/6"07.8,+/"*."+*,-";.*3/+32?%7))(%'0.-)% 8)&),()5%.*%,.(%6*)/),( ',5%0)".3%6"#$/ A??B7..(% 6"#$ C)/ :0;%!*.5#8(+.,%8'/+,$% 8)&),()5 E??B7..(% 6"#$ C)/ :8;%!),)(*'()/%3.*G'0")% 8.'"%/)'&/%:/)8(+.,% AH:';:A;2%:E;2%:H;%.*%:I;%.7% (1)%<.'"%',5%9'/%J)/.#*8)% <..*5+,'(+.,%K8(; @+"")5%3+(1% ,.,6.*.#/% &'()*+'" =+,+&#&%.7%E?%(.%A??%7))(% ,.%)F6"+8+(% '0.-)%',5%>?%7))(%0)".3% /('()&),( 6"#$/ <)&),(%.*% D)6'*(&),(B '66*.-)5% &)81',+8'"%6"#$ LF6',5+,$% =+,+&#&%.7%E?%(.%>?%7))(% 8)&),(%6"#$2% '0.-)%',5%0)".3%6"#$ 0),(.,+()%$)" @+"")5%3+(1% A??B7..(% 0),(.,+()%$)"%+,% 6"#$ /.&)%8'/)/ !"#$%&="()"+*,-",./,0"1234"07.8,+/"*."+*,-";.*3/+32/."*."+/5/)3 :';%!*.5#8(+.,%8'/+,$%,.(% <)&),( =+,+&#&%.7%>?%7))(%'0.-)% @+"")5%3+(1% A??B7..(% 8)&),()5%.*%,.(%6*)/),( ',5%0)".3%6"#$/ ,.,6.*.#/% 6"#$ &'()*+'" :0;%!*.5#8(+.,%8'/+,$% 8)&),()5 <)&),(%.*% =+,+&#&%.7%E?%(.%A??%7))(% D)6'*(&),(B '0.-)%',5%>?%7))(%0)".3%7.*% '66*.-)5% 6"#$/ &)81',+8'"%6"#$ :8;%!),)(*'()/%3.*G'0")% LF6',5+,$% =+,+&#&%.7%E?%(.%>?%7))(% 8.'"%/)'&/%:AH:';:E;%.*%:I;% 8)&),(%6"#$2% '0.-)%',5%0)".3%6"#$ .7%(1)%<.'"%',5%9'/% 0),(.,+()%$)" J)/.#*8)%<..*5+,'(+.,% K8(; !"#$%&?"()")*)+*,-",./,0"1234"07.8,+/"+,02)9")*3"+/5/)3/6"*.")*3";./0/)3 :';%!*.5#8(+.,%8'/+,$%,.(% <)&),( =+,+&#&%.7%>?%7))(%'0.-)% 8)&),()5%.*%,.(%6*)/),( ',5%0)".3% E??B7..(% 6"#$ C)/ @+"")5%3+(1% A??B7..(% 0),(.,+()%$)"%+,% 6"#$ /.&)%8'/)/ C)/ >?B7..(%6"#$% M. '(%&+,+&#& ,.%)F6"+8+(% /('()&),( >?B7..(%6"#$% M. '(%&+,+&#& @+"")5%3+(1% ,.,6.*.#/% &'()*+'" >?B7..(%6"#$% M. '(%&+,+&#& :0;:E;%4,%8)&),()5%6.*(+.,% ,.%)F6"+8+(% .7%6*.5#8(+.,%8'/+,$2% /('()&),( &+,+&#&%.7%A??%7))(%'0.-)% ',5%>?%7))(%0)".3%7.*% 8)&),(%6"#$/ >?B7..(%6"#$% M. '(%&+,+&#& <)&),(%.*% =+,+&#&%.7%E?%(.%A??%7))(% D)6'*(&),(B '0.-)%',5%>?%7))(%0)".3%7.*% '66*.-)5% 8)&),(%6"#$/ &)81',+8'"%6"#$ !"#$%&'"()")*)+*,-",./,0"1234"+/5/)3/6"07.8,+/"+,02)9" :';%!*.5#8(+.,%8'/+,$%,.(% <)&),( =+,+&#&%.7%>?%7))(%'0.-)% 8)&),()5%.*%,.(%6*)/),( ',5%0)".3% <)&),(%.*% D)6'*(&),(B '66*.-)5% &)81',+8'"%6"#$ C)/ @+"")5%3+(1% ,.,6.*.#/% &'()*+'" :0;%!*.5#8(+.,%8'/+,$% 8)&),()5 :0;%!*.5#8(+.,%8'/+,$% 8)&),()5 ,.%)F6"+8+(% /('()&),( C)/%+,% :A;2%:E;2% :I;% 102 Table 5.6: Summary of alternative energy resource types included in Pennsylvania’s Alternative Energy Portfolio Standard (AEPS). The first and second column are as presented in the Pennsylvania AEPS Alternative Energy Credit Program website and in the annual reports [72, 73] !"#$% !*+(),-+'.(%",()/0% &'() 1(234)5(%&06(728 B B B B E*-5F%G'H43) K4(*%9(** 93-*%M',(%M(+:-,( O0C)3 B G-,C@'**%D-2 B B R+:()%E'3N-22%D-2 R+:()%D-2 BI% $3*-)% $3*-) B B BB T',C T33CVT33C%T-2+(% $3*'C2 T33CVT33C%T-2+(% $3*'C2%W%E*-5F%G'H43) E*-5F%G'H43) BB E*-2+%K4),-5(%D-2 B BB BB BB BB BB 9:-6+()%;<=>%?(@','+'3,2A% !*+(),-+'.(%(,()/0% 234)5(2 !"##$%&#'()**%+,+-./ 7'L8%K4(*%5(**2 7L'8%93-*%N',(%N(+:-,( 7.8%G3PQ'N6-5+% :0C)363P() 7.'''8%E'3*3/'5-**0%C()'.(C% N(+:-,(%/-2 7.''8%E'3N-22%(,()/0 !"##$%&#'()**%+,+-./2%!"###$% &#'3'.#4)33/%5+-#"+5% (+67),+%.)* 7'8%$3*-)%6:3+3.3*+-'5%3)% 3+:()%23*-)%(*(5+)'5% (,()/0S%7''8%$3*-)%+:()N-*% (,()/0 7'''8%T',C%63P() !"##$%&#'()**%+,+-./ 9:-6+()%;<=>%?(@','+'3,2A%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% $34)5(2%-2% &'()%B%3)%BB%-*+(),-+'.(%(,()/0%234)5( C(@',(C%',% D!&$ 0#+-%1%!"##$%&#'()**%+,+-./ EGI%EGJ &'()%B%7.'8%K4(*%5(**2 K9 &'()%B%7.'''8%93-*%N',(%N(+:-,( 9MMI%9MD &'()%B%7'''8%G3PQ'N6-5+%:0C)363P() O#$ BB T33CVT33C%T-2+(% $3'*C2 BB T33CVT33C%T-2+(% $3'*C2%W%E*-5F%G'H43)Z &'()%B%7.''8%E'3N-22%(,()/0 RED 0#+-%1%!"$%&#'3'.#4)33/%5+-#"+5%(+67),+%.)*2% RD !"##$%&#'()**%+,+-./ $RG &'()%B%7''8%T',C%63P() 0#+-%1%!"##$%&#'()**%+,+-./ TU? T?$ !"##$%&#'()**%+,+-./ 0#+-%1%!"##$%&#'()**%+,+-./ T?$I%EGI%EGJ !"##$%&#'()**%+,+-./ &'()%BB%7.'8%D(,()-+'3,%3@%(*(5+)'5'+0%4+'*'X',/% EGI%EGJ Y0Q6)3C45+2%3@%+:(%64*6',/%6)35(22%-,C% P33C%N-,4@-5+4)',/%6)35(22 0#+-%11%!###$%9+(),5:*#5+%(),).+(+,6 EKD ?'2+)'Y4+(C%D(,()-+'3, 7L'''8%?'2+)'Y4+(C% /(,()-+'3,%202+(N2 O0C)3 7'.8%G-)/(Q25-*(% :0C)363P() M4,'5'6-*%$3*'C%T-2+( !"##$%&#'()**%+,+-./ R+:()%D-2(2 !8##$%9+(),5:*#5+% (),).+(+,6 T-2+(%93-* T-2+(%O(-+ GKD &'()%B%7'8%$3*-)%6:3+3.3*+-'5%-,C%23*-)% +:()N-*%(,()/0 !8##$%9+(),5:*#5+% (),).+(+,6 E*-2+%K4),-5(%W%R+:()% !8##$%9+(),5:*#5+% D-2(2 (),).+(+,6 BB BB &'()%B%7.8%E'3*3/'5-**0%C()'.(C%N(+:-,(%/-2 7L8%T-2+(%53-* !8##$%9+(),5:*#5+% (),).+(+,6 !"##$%&#'()**%+,+-./ !"##$%&#'()**%+,+-./ 0#+-%11%!###$%9+(),5:*#5+%(),).+(+,62%!"##$% EKDI%RD 1,6+.-)6+5%4'(;#,+5%4')3%.)*#<#4)6#',% 6+47,'3'./ &'()%BB%7''8%?'2+)'Y4+(C%/(,()-+'3,%202+(N2 &'()%BB%7'.8%G-)/(Q25-*(%:0C)363P() T!& &'()%BB%7.8%M4,'5'6-*%23*'C%P-2+( 0#+-%11%!###$%9+(),5:*#5+%(),).+(+,62%!"##$% 1,6+.-)6+5%4';#,+5%4')3%.)*#<#4)6#',% 6+47,'3'./ &'()%BB%7'8%T-2+(%53-* 0#+-%11%!###$%9+(),5:*#5+%(),).+(+,6 M$T ?$1I%RD T9 TO &'()%BB%7.'8%D(,()-+'3,%3@%(*(5+)'5'+0%4+'*'X',/% T?$ Y0Q6)3C45+2%3@%+:(%64*6',/%6)35(22%-,C% P33C%N-,4@-5+4)',/%6)35(22 &'()%BB%7.'8%D(,()-+'3,%3@%(*(5+)'5'+0%4+'*'X',/% T?$I%EGI%EGJ Y0Q6)3C45+2%3@%+:(%64*6',/%6)35(22%-,C% P33C%N-,4@-5+4)',/%6)35(22 103 Table 5.7: Alternative energy sources for Pennsylvania in PJM-EIS Generation Attribute Tracking System (GATS) !""#$%&'(&)* 452 48 48; <=2 <== A+B 5< DE+ 852 =+G I42 I2 +I8 G!3 G< GA+ GD GJA +),#-$./0$1&*$0/&*/2!3+ 46'.(/1,#*'-$/7'. 46'-9/6&:,)# 46'-9/6&:,)# <)'6/>&*$/>$(?'*$@7'. <)'6/>&*$/>$(?'*$ A$>'*0/.&0$/#$.C)*.$ 5,$6/-$66 DF0#)/C,>C/.()#'7$ 8'*01&66/7'. =,*&-&C'6/.)6&0/H'.($ I(?$#/"&)>'../7'. I(?$#/7'. +)6'#/C?)()%)6('&<)*%$*(&)*'6/?F0#) G'.($/-)'6/@/)(?$#/-)'6 G))0@H))0/H'.($/.)6&0. G'.($/?$'( G&*0 104 Chapter 6 Concluding Remarks Subsurface leakage pathways, abandoned oil and gas wells and faults, and their role in GHG emission estimates and mitigation are considered from a modeling, measurement, and policy perspective in this dissertation. The GHGs considered are CO2 and methane and the mitigation strategies considered are CO2 capture and storage in geologic formations and methane emissions reduction from abandoned oil and gas wells. The modeling and measurement work (Chapters 2, 3 and 4) focuses on understanding and determination of leakage potentials. Policies considered focus on methane emissions from abandoned oil and gas wells, specifically for Pennsylvania (Chapter 5). In terms of modeling, analytical solutions that capture vertical flow effects and consider properties of leakage pathways are found to be valuable for determining leakage rates. These analytical solutions can be used in a multi-scale framework to represent the local-scale effects of a leakage feature in a basin-scale model, as described in Chapter 3 for leaky faults. The solutions in Chapter 2 are developed for faults but much of the theory can be applied to abandoned wells. These models provide the capability to understand and quantify the risks of leakage through subsurface pathways. 105 Determining the leakage potential of subsurface leakage pathways requires field measurements and an understanding of the properties of the leakage feature. For abandoned wells, especially those that lack available records, direct measurement of methane fluxes from these wells at the land surface provides information at a low cost. Based on the orders of magnitude variation in the measured methane fluxes from these wells, the properties of these wells are likely to be highly variable. Further, the composition of the emitted gas, which includes ethane, propane, and butane isomers, is highly variable. Therefore, more measurements and field investigations of these wells are needed to determine the controlling processes and to develop useful modeling tools. The measured fluxes can also be used to determine the methane emissions potential of these wells. Methane is a GHG whose mitigation can have significant climate benefits in the short term. However, methane emissions and the sources of these emissions are highly uncertain. The measurements in Chapter 4 combined with the number of abandoned oil and gas wells estimated in Chapter 5 show that methane emissions from abandoned oil and gas wells can be 4 to 13% of total currently estimated anthropogenic methane emissions for the state of Pennsylvania. Similar numbers might be extrapolated to other states with a significant history of oil and gas development and production. The ability to quantify methane emissions from abandoned wells are necessary for future inclusion of leaky wells as a methane emissions source in emission inventories. It is also necessary to develop appropriate mitigation strategies. For accurate upscaling of these fluxes to emissions at the state and national levels, additional measurements are needed to characterize their distribution and determine the attributes that control leakage rates. Improving methane emission estimates from abandoned oil and gas wells requires not only measurements but also policy actions. Policy actions are needed to promote reporting and monitoring of abandoned oil and gas wells and also to promote 106 mitigation that results in methane emissions reduction. More research is needed to determine if current regulations for well abandonment are sufficient for methane emissions reductions and if alternative strategies should be considered. The use of gases emitted from abandoned wells and the qualification of these wells as an alternative energy resource may incentivize mitigation as wells as monitoring and reporting. 6.1 6.1.1 Future work Model development The analytical model for leakage through faults presented in Chapter 2 is designed for implementation within a multi-scale framework described in Chapter 3. The next step would be to further develop, test, and validate the multi-scale model for leakage through faults considering a wider range of fault configurations. This model will be valuable for considering CO2 leakage risk in basins containing faults and being considered for geologic storage of CO2 such as the Illinois Basin. Development of models of abandoned oil and gas wells and determination of the properties that govern flow processes are important for both estimating emissions and developing strategies for mitigation. Models can be used to improve measurement plans that inform emission estimates and optimize allocation of resources. An improved understanding of the drivers behind fluid leakage, based on model development coupled with field data, is also useful for mitigating potential fluid leakage from abandoned oil and gas wells to groundwater resources and to the atmosphere. 6.1.2 Measurements and policies Available measurements of methane emissions from abandoned oil and gas wells are limited to the 19 wells in northwest Pennsylvania presented in Chapter 4. However, abandoned oil and gas wells can be found in other regions within Pennsylvania and 107 in other states and countries. First, it would be interesting to measure abandoned oil and gas wells in southwest Pennsylvania where many abandoned wells exist and active shale gas production is occurring. These measurements would be in contrast to the measurements in northwest Pennsylvania where shale gas production has been minimal but many abandoned wells exist. Given the current and projected growth in shale gas development in the U.S. and the interest in developing shale gas resources around the world (e.g. China), it is important to understand the interaction between abandoned oil and gas wells and shale gas development. There are many other states such as California, Colorado, and Texas with relatively long histories of oil and gas production that face issues similar to those in Pennsylvania. Additional measurements should be performed in these states and measurement plans should be designed to inform policies specific to the corresponding regions. Thus, policies relevant to oil and gas well abandonment, reporting and monitoring, groundwater protection, air quality control, and any regulations on carbon or climate in these states should be studied. 108 Bibliography [1] David T. Allen, Vincent M. Torres, James Thomas, David W. Sullivan, Matthew Harrison, Al Hendler, Scott C. Herndon, Charles E. Kolb, Matthew P. Fraser, A. Daniel Hill, Brian K. Lamb, Jennifer Miskimins, Robert F. Sawyer, and John H. Seinfeld. Measurements of methane emissions at natural gas production sites in the united states. Proceedings of the National Academy of Sciences, 2013. [2] Anne E Altor and William J Mitsch. Methane flux from created riparian marshes: relationship to intermittent versus continuous inundation and emergent macrophytes. Ecological Engineering, 28(3):224–234, 2006. [3] American Refining Group, Inc. Bradford oil field history. www.amref.com, 2014. [4] Erik I. Anderson. Analytical solutions for flow to a well through a fault. Advances in Water Resources, 29:1790–1803, 2006. [5] Erik I. Anderson and M. Bakker. Groundwater flow through anisotropic fault zones in multiaquifer systems. Water Resources Research, 44, 2008. [6] Ralph Arnold and William J. Kemnitzer. Petroleum in the United States and Possessions. Harper and Brothers, New York and London, 1931. [7] George H. Ashley and J. French Robinson. The Oil and Gas Fields of Pennsylvania, volume 1 of Pennsylvania Geological Survey Fourth Series. Pennsylvania Department of Internal Affairs, Bureau of Topographic and Geological Survey, 1922. [8] Stefan Bachu and William D. Gunter. Acid-gas injection in the alberta basin, canada: a co2-storage experience. Geological Society, London, Special Publications, 233(1):225–234, 2004. [9] Jacob Bear. Dynamics of fluids in porous media. Elsevier, New York, 1972. [10] Brian Berkowitz. Characterizing flow and transport in fractured geological media: A review. Advances in Water Resources, 25(8-12):861 – 884, 2002. [11] A. R. Brandt, G. A. Heath, E. A. Kort, F. O’Sullivan, G. P´etron, S. M. Jordaan, P. Tans, J. Wilcox, A. M. Gopstein, D. Arent, S. Wofsy, N. J. Brown, R. Bradley, G. D. Stucky, D. Eardley, and R. Harriss. Methane leaks from north american natural gas systems. Science, 343(6172):733–735, 2014. 109 [12] Michael A. Celia, Jan M. Nordbotten, Benjamin Court, Mark Dobossy, and Stefan Bachu. Field-scale application of a semi-analytical model for estimation of CO2 and brine leakage along old wells. International Journal of Greenhouse Gas Control, 5(2):257 – 269, 2011. [13] Y. Chen, Kevin. K. Lehmann, J. Kessler, B. Sherwood Lollar, G. Lacrampe Couloume, and T. C. Onstott. Measurement of the 13c/12c of atmospheric ch4 using near-infrared (nir) cavity ring-down spectroscopy. Analytical Chemistry, 85(23):11250–11257, 2013. [14] G.V. Chilingar and B. Endres. Environmental hazards posed by the los angeles basin urban oilfields: an historical perspective of lessons learned. Environmental Geology, 47(2):302–317, October 2005. [15] K. H. Coats, J. R. Dempsey, and J.H. Henderson. The use of vertical equilibrium in two-dimensional simulation of three-dimensional reservoir performance. SPE Journal, 11(1):63–71, March 1971. [16] Benjamin Court, Karl W. Bandilla, Michael A. Celia, Adam Janzen, Mark Dobossy, and Jan M. Nordbotten. Applicability of vertical-equilibrium and sharpinterface assumptions in CO2 sequestration modeling. International Journal of Greenhouse Gas Control, 10:134–147, 2012. [17] Cheryl L. Cozart and John A. Harper. Oil and gas development in pennsylvania. Progress Report 205, Department of Environmental Resources, Pennsylvania Geological Survey, 1993. [18] E.A Davidson, K Savage, L.V Verchot, and Rosa Navarro. Minimizing artifacts and biases in chamber-based measurements of soil respiration. Agricultural and Forest Meteorology, 113(1–4):21 – 37, 2002. ¡ce:title¿FLUXNET 2000 Synthesis¡/ce:title¿. [19] R. De Loubens and T.S. Ramakrishnan. Analysis and computation of gravityinduce migration in porous media. Journal of Fluid Mechanics, 675:60–86, 2011. [20] C. L. Farmer. Upscaling: a review. International Journal for Numerical Methods in Fluids, 40(1-2):63–78, 2002. [21] D.R. Faulkner, C.A.L. Jackson, R.J. Lunn, R.W. Schlische, Z.K. Shipton, C.A.J. Wibberley, and M.O. Withjack. A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones. Journal of Structural Geology, 32(11):1557 – 1575, 2010. [22] Chas R. Fettke. Oil and gas development in pennsylvania in 1952. Progress Report 143, Department of Environmental Resources, Pennsylvania Geological Survey, 1953. [23] Q.J. Fisher and S.J. Jolley. Treatment of faults in production simulation models. Geological Society, London, Special Publications, 292:219–233, 2007. 110 [24] Sarah E. Gasda, Jan M. Nordbotten, and Michael A. Celia. Vertical equilibrium with sub-scale analytical methods for geological CO2 sequestration. Computational Geosciences, 13(4):469–481, 2009. [25] SarahE. Gasda, Stefan Bachu, and MichaelA. Celia. Spatial characterization of the location of potentially leaky wells penetrating a deep saline aquifer in a mature sedimentary basin. Environmental Geology, 46(6-7):707–720, 2004. [26] Anthony W. Gorody. Factors affecting the variability of stray gas concentration and composition in groundwater. Environmental Geosciences, 19(1):17–31, 2012. [27] A.E. Gurevich, B.L. Endres, J.O. Robertson Jr., and G.V. Chilingar. Gas migration from oil and gas fields and associated hazards. Journal of Petroleum Science and Engineering, 9(3):223 – 238, 1993. [28] William C. Haneberg. Steady state groundwater flow across idealized faults. Water Resources Research, 31(7):1815–1820, July 1995. [29] E. C. D. Hooper. Fluid migration along growth faults in compacting sediments. Journal of Petroleum Geology, 14:161–180, 1991. [30] Robert W. Howarth, Anthony Ingraffea, and Terry Engelder. Should fracking stop? Nature, 477(7364):271–275, 09 2011. Natural gas: [31] Robert W. Howarth, Renne Santoro, and Anthony Ingraffea. Methane and the greenhouse-gas footprint of natural gas from shale formations. Climatic Change, 106:679–690, 2011. [32] Ying-Kuang Hsu, Tony VanCuren, Seong Park, Chris Jakober, Jorn Herner, Michael FitzGibbon, Donald R. Blake, and David D. Parrish. Methane emissions inventory verification in southern california. Atmospheric Environment, 44(1):1 – 7, 2010. [33] Herbert E. Huppert and Andrew W. Woods. Gravity-driven flows in porous layers. Journal of Fluid Mechanics, 292:55–69, 1995. [34] Intergovernmental Panel on Climate Change. Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Technical report, Intergovernmental Panel on Climate Change, Cambridge, United Kingdom and New York, NY, USA, 2007. [35] Intergovernmental Panel on Climate Change. Summary for policymakers. in: Climate change 2013: The physical science basis. contribution of working group i to the fifth assessment report of the intergovernmental panel on climate change. Technical report, Intergovernmental Panel on Climate Change, Cambridge, United Kingdom and New York, NY, USA., 2013. 111 [36] International Energy Agency. Technology roadmap: Carbon capture and storage. Technical report, International Energy Agency, 2013. [37] IOGCC. Orphaned and abandoned wells: innovative solutions. The Journal of the Interstate Oil & Gas Compact Commision, Groundwork, October 2009. [38] R.E. Jackson, A.W. Gorody, B. Mayer, J.W. Roy, M.C. Ryan, and D.R. Van Stempvoort. Groundwater protection and unconventional gas extraction: The critical need for field-based hydrogeological research. Groundwater, 51(4):488–510, 2013. [39] Robert B. Jackson, Avner Vengosh, Thomas H. Darrah, Nathaniel R. Warner, Adrian Down, Robert J. Poreda, Stephen G. Osborn, Kaiguang Zhao, and Jonathan D. Karr. Increased stray gas abundance in a subset of drinking water wells near marcellus shale gas extraction. Proceedings of the National Academy of Sciences, 110(28):11250—11255, July 2013. [40] Rachhpal S. Jassal, T. Andrew Black, Zoran Nesic, and David Gaumont-Guay. Using automated non-steady-state chamber systems for making continuous longterm measurements of soil {CO2} efflux in forest ecosystems. Agricultural and Forest Meteorology, 161(0):57 – 65, 2012. [41] PD Jenden, DJ Drazan, and IR Kaplan. Mixing of thermogenic natural gases in northern appalachian basin. AAPG Bulletin, 77(6):980–998, 1993. [42] V.T. Jones and R.J. Drozd. Prediction of oil and gas potential by near-surface geochemistry. AAPG Bulletin, 67(6):932–952, 1983. [43] George E. King and Daniel E. King. Environmental risk arising from well construction failure: Difference between barrier and well failure, and estimates of failure frequency across common well types, locations and well age. SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, USA, SPE 166142, September 2013. [44] Larry Lake. Enhanced Oil Recovery. Prentice Hall, 1 edition, June 1996. [45] Jean Le Mer and Pierre Roger. Production, oxidation, emission and consumption of methane by soils: A review. European Journal of Soil Biology, 37(1):25 – 50, 2001. [46] Jennifer L. Lewicki, Jens Birkholzer, and Chin-Fu Tsang. Natural and industrial analogues for leakage of CO2 from storage reservoirs: identification of features, events, and processes and lessons learned. Environmental Geology, 52(3):457– 467, 2007. [47] G.P. Livingston and G.L. Hutchinson. Enclosure-based measurement of trace gas exchange: applications and sources of error, chapter 2, pages 14–51. Methods in Ecology. Blackwell Science Ltd., 1995. 112 ¨ Mander, Evelyn Uuemaa, Ain Kull, Arno Kanal, Martin Maddison, Kaido [48] Ulo Soosaar, J¨ uri-Ott Salm, Merje Lesta, Raili Hansen, Reili Kuller, et al. Assessment of methane and nitrous oxide fluxes in rural landscapes. Landscape and Urban Planning, 98(3):172–181, 2010. [49] T. Manzocchi, A.E. Heath, J.J. Walsh, and C. Childs. The representation of two phase fault-rock properties in flow simulation models. Petroleum Geoscience, 8(2):119–132, 2002. [50] T. Manzocchi, J.J. Walsh, P. Nell, and G. Yielding. Fault transmissibility multipliers for flow simulation models. Petroleum Geoscience, 5:53–63, 1999. [51] Bert Metz, Ogunlade Davidson, Heleen de Coninck, Manuela Loos, and Leo Meyer. Carbon dioxide capture and storage. Technical report, Intergovernmental Panel on Climate Change (IPCC), 2005. [52] Cass T. Miller, Clint N. Dawson, Matthew W. Farthing, Thomas Y. Hou, Jingfang Huang, Christopher E. Kees, C.T. Kelley, and Hans Petter Langtangen. Numerical simulation of water resources problems: Models, methods, and trends. Advances in Water Resources, 51(0):405 – 437, 2013. ¡ce:title¿35th Year Anniversary Issue¡/ce:title¿. [53] Lisa J. Molofsky, John A. Connor, Albert S. Wylie, Tom Wagner, and Shahla K. Farhat. Evaluation of methane sources in groundwater in northeastern pennsylvania. Groundwater, 51(3):333–349, 2013. [54] A Mosier, D Schimel, D Valentine, K Bronson, and W Parton. Methane and nitrous oxide fluxes in native, fertilized and cultivated grasslands. Nature, 350(6316):330–332, 1991. [55] Juan P. Nogues, Benjamin Court, Mark Dobossy, Jan M. Nordbotten, and Michael A. Celia. A methodology to estimate maximum probable leakage along old wells in a geological sequestration operation. International Journal of Greenhouse Gas Control, 7(0):39 – 47, 2012. [56] Jan M. Nordbotten and Michael A. Celia. An improved analytical solution for interface upconing around a well. Water Resources Research, 42, 2006. [57] Jan M. Nordbotten and Michael A. Celia. Geological Storage of CO2 : Modeling Approaches for Large-Scale Simulation. Wiley, Hoboken, New Jersey, 1 edition, December 2012. [58] Jan M. Nordbotten, Michael A. Celia, and Stefan Bachu. Analytical solutions for leakage rates through abandoned wells. Water Resources Research, 40:W04204, 2004. [59] Jan M. Nordbotten and Helge K. Dahle. Impact of the capillary fringe in vertically integrated models for CO2 storage. Water Resources Research, 47(W02537):11, 2011. 113 [60] Jan M. Nordbotten, Dmitri Kavetski, Michael A. Celia, and Stefan Bachu. Model for CO2 Leakage Including Multiple Geological Layers and Multiple Leaky Wells. Environ. Sci. Technol, 43:743–749, 2009. [61] Oil and Gas Technical Advisory Board (TAB). Draft minutes oil and gas technical advisory board meeting. files.dep.state.pa.us, August 2012. [62] Stephen G. Osborn, Avner Vengosh, Nathaniel R. Warner, and Robert B. Jackson. Methane contamination of drinking water accompanying gas-well drilling and hydraulic fracturing. Proceedings of the National Academy of Sciences, 108(20):8172–8176, 2011. [63] PA DEP. Pennsylvania’s Plan For Addressing Problem Abandoned Wells and Orphaned Wells. Document Number 550-0800-001, April 2000. [64] PA DEP. Pennsylvania climate change action plan update. Technical report, PA DEP, December 2013. [65] PA DEP. Permits Issued-Wells Drilled Maps. www.portal.state.pa.us, April 2014. [66] S. Pacala and R. Socolow. Stabilization wedges: Solving the climate problem for the next 50 years with current technologies. Science, 305(5686):968–972, 2004. [67] Donald W. Peaceman. Interpretation of well-block pressures in numerical reservoir simulation. SPE Journal, (18):183–194, 1978. [68] Mark Person, Amlan Banerjee, John Rupp, Cristian Medina, Peter Lichtner, Carl Gable, Rajesh Pawar, Michael Celia, Jennifer McIntosh, and Victor Bense. Assessment of basin-scale hydrologic impacts of CO2 sequestration, Illinois basin. International Journal of Greenhouse Gas Control, 4:840–854, 2010. [69] Gabrielle P´etron, Gregory Frost, Benjamin R. Miller, Adam I. Hirsch, Stephen A. Montzka, Anna Karion, Michael Trainer, Colm Sweeney, Arlyn E. Andrews, Lloyd Miller, Jonathan Kofler, Amnon Bar-Ilan, Ed J. Dlugokencky, Laura Patrick, Charles T. Moore Jr., Thomas B. Ryerson, Carolina Siso, William Kolodzey, Patricia M. Lang, Thomas Conway, Paul Novelli, Kenneth Masarie, Bradley Hall, Douglas Guenther, Duane Kitzis, John Miller, David Welsh, Dan Wolfe, William Neff, and Pieter Tans. Hydrocarbon emissions characterization in the colorado front range: A pilot study. Journal of Geophysical Research, 117(D04304):19, 2012. [70] George F. Pinder and Michael A. Celia. Subsurface Hydrology. Wiley Interscience, 2006. [71] David Pritchard. Gravity currents over fractured substrates in a porous medium. Journal of Fluid Mechanics, 584:415–431, 2007. 114 [72] PUC Bureau of Conservation, Economics and Energy Planning and PA DEP. 2007 annual report: Alternative energy portfolio standards act of 2004. Technical report, Pennsylvania Public Utility Commission, P.O. Box 3265, Harrisburg, PA 17105-3265, May 2008. [73] PUC Bureau of Technical Utility Services and PA DEP. 2012 Annual Report: Alternative Energy Portfolio Standards Act of 2004. Technical report, Pennsylvania Public Utility Commission, P.O. Box 3265, Harrisburg, PA 17105-3265, October 2013. [74] M. B. Rayment. Closed chamber systems underestimate soil CO2 efflux. European Journal of Soil Science, 51(1):107–110, 2000. [75] M. C. Reid, R. Tripathee, K. V. R. Sch¨afer, and P. R. Jaff´e. Tidal marsh methane dynamics: Difference in seasonal lags in emissions driven by storage in vegetated versus unvegetated sediments. Journal of Geophysical Research: Biogeosciences, pages 1802–1813, 2013. [76] K. Savage, E. A. Davidson, and A. D. Richardson. A conceptual and practical approach to data quality and analysis procedures for high-frequency soil respiration measurements. Functional Ecology, 22(6):1000–1007, 2008. [77] Martin Schoell. Multiple origins of methane in the earth. Chemical Geology, 71(1–3):1 – 10, 1988. [78] Chao Shan, Iraj Javandel, and PaulA. Witherspoon. Characterization of leaky faults: Study of water flow in aquifer-fault-aquifer systems. Water Resources Research, 31(12):2897–2904, December 1995. [79] Drew Shindell, Johan C. I. Kuylenstierna, Elisabetta Vignati, Rita van Dingenen, Markus Amann, Zbigniew Klimont, Susan C. Anenberg, Nicholas Muller, Greet Janssens-Maenhout, Frank Raes, Joel Schwartz, Greg Faluvegi, Luca Pozzoli, Kaarle Kupiainen, Lena H¨oglund-Isaksson, Lisa Emberson, David Streets, V. Ramanathan, Kevin Hicks, N. T. Kim Oanh, George Milly, Martin Williams, Volodymyr Demkine, and David Fowler. Simultaneously mitigating near-term climate change and improving human health and food security. Science, 335(6065):183–189, 2012. [80] Qing Tao, David Alexander, and Steven L. Bryant. Development and field application of model for leakage of CO2 along a fault. Energy Procedia 00, 2013. [81] SW Taylor, B Sherwood Lollar, and I Wassenaar. Bacteriogenic ethane in nearsurface aquifers: Implications for leaking hydrocarbon well bores. Environmental science & technology, 34(22):4727–4732, 2000. [82] C-F. Tsang and J.A. Apps, editors. Deep Injection Disposal of Hazardous and Industrial Waste: Scientific and Engineering Aspects. Academic Press, San Diego, CA, 1996. 115 [83] Chin-Fu Tsang, Jens Birkholzer, and Jonny Rutqvist. A comparative review of hydrologic issues involved in geologic storage of CO2 and injection disposal of liquid waste. Environmental Geology, 54:1723–1737, 2008. 10.1007/s00254-0070949-6. [84] UNEP. Near-term climate protection and clean air benefits: Actions for controlling short-lived climate forcers. A unep synthesis report, United Nations Environment Programme (UNEP), Nairobi, Kenya, 2011. [85] U.S. Energy Information Administration. International Energy Outlook 2013 With Projections to 2040. Technical report, U.S. Energy Information Administration, Office of Energy Analysis, U.S. Department of Energy, Washington, DC 20585, July 2013. [86] Dominic Vella, Jerome A. Neufeld, Herbert E. Huppert, and John R. Lister. Leakage from gravity currents in a porous medium. part 2. a line sink. Journal of Fluid Mechanics, 666:414–427, January 2011. [87] Theresa L. Watson and Stefan Bachu. Evaluation of the potential for gas and CO2 leakage along wellbores. SPE Drilling and Completion, SPE 106817, March 2009. [88] Y. C. Yortsos. A theoretical analysis of vertical flow equilibrium. Transport in Porous Media, 18:107–129, 1995. [89] Mehdi Zeidouni. Analytical model of leakage through fault to overlying formations. Water Resources Research, 48, 2012. [90] Quanlin Zhou, Jens T. Birkholzer, Edward Mehnert, Yu-Feng Lin, and Keni Zhang. Modeling basin- and plume-scale processes of CO2 storage for full-scale deployment. Ground Water, 48(4):494–514, 2010. 116