University of Southern California Dissertations and Theses

Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure

(USC Thesis Other)

Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content UNIVERSITY OF SOUTHERN CALIFORNIA Hydraulic Fracturing and the Environment: Risk Assessment for Groundwater Contamination from Well Casing Failure By Nima Jabbari A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF ENGINEERING (ENVIRONMENTAL ENGINEERING) LOS ANGELES, CALIFORNIA May 2016 © NIMA JABBARI 2016 ii Abstract Hydraulic fracturing is a well stimulation method used within the oil and gas industry to increase production. The recent developments of horizontal drilling, in combination with fracturing, have made this technique very popular for natural gas production from low-permeability formations such as tight sands and shales. Nevertheless, this technology is relatively new and the impact on the environment has been called into question, as the long-term effects are not thoroughly understood. From a technical standpoint, fracturing tight formations is accomplished by introducing highly pressurized fluid (fracturing slurry) into the formation to create fissures and cracks, which in turn increase the permeability. Fracturing slurry is a mixture of water, propping agent (mostly sand), and a combination of various toxic and benign chemicals which are responsible for reducing friction, inhibiting corrosion, and preventing scale. One of the main points of contention in fracturing concerns the specific chemical compounds that are being used. Debates regarding the risk versus reward ratio of fracturing are ongoing within the scientific community. Policy makers and the public have also expressed concerns about possible groundwater contamination during the hydraulic fracturing process. This study will take a close look at some of those concerns. Specifically, out of several risk pathways to the aquifer, leakage from a well casing during the injection periods has been addressed. Computer simulations of fluid flow in porous media and transport phenomena have been conducted and the sensitivity of results with respect to parameters of interest has been evaluated. In order to consider the risk to human health, Tetraethylenepentamine has been selected as the representative chemical with proven toxicity and with a pre-defined reference dose in literature, so that this agent can be used as a threshold to compare the results of various hypothetical scenarios. The uncertainty of the problem is addressed through a novel risk iii framework consisting of the numerical simulation of fluid flow and chemical transport in porous media and the Monte Carlo simulation. Hydrogeological and operational parameters are selected as random variables in the framework. This unique risk framework is developed to assess three environmental performance metrics: the solute concentration, the arrival times from source to receptor, and the ingestion hazard of the contaminated aquifer. The effect of parametric uncertainty in risk is also analyzed. The results show that risk strongly depends on water saturation, the anisotropy, and the heterogeneity of the permeability field. Initial water saturation of 0.3 and vertical to horizontal permeability ratio of 0.1 within the sand medium are associated with more adverse situation. Through the sensitivity analysis, it is indicated that leakage depth and leakage rate are more effective in controlling the measured risk value when compared to the hydrogeological properties. On the same line of thoughts, the critical depth for extending the surface casing below the depth of fresh water aquifer is introduced and calculated as an important result from the simulation sets. For the case of this study, the critical depth is 18 m below the base of freshwater. Results of the parametric uncertainty analysis for anisotropic homogenous scenarios show that on average it takes 276 years for the chemical to show concentrations above the threshold. Also, the hazard index is reported to be 0.45. However, introducing spatial uncertainty to the sand layer alleviates the situation by increasing the time and reducing the risk. Findings of this study can be applied by setting more stringent requirements for the well integrity to ensure that hydraulic fracturing is carried out in an environmentally safe and sound manner. In conclusion, this research will examine the likelihood of groundwater contamination caused by a specific failure in the well casing during a hydraulic fracturing job. However, in order to investigate the worst case scenario, the presented results are drawn from a specific situation in iv which the leakage is ongoing all over the injection phase of the operation. One implementation of the findings in this study is in enhancing oil and gas production standards, where it can be used to refine the surface casing requirements. In addition, the risk framework proposed by this study can be implemented in a real fracturing operation after enriching by real data. This practice will give site managers the power of making appropriate decisions in a timely manner in an event of an unwanted failure. It should be noted that the outlined concluding remarks are specifically related to the failure scenario, conceptual model, and geological settings discussed in this study. The results cannot and should not be generalized to other failure scenarios; however the approach taken here to construct the risk framework can be used in subsurface injection operations practiced in regions with more complex geological settings and operational parameters. Future research activities will comprise of increasing the uncertainty dimension of the failure scenario and the risk assessment by including more stochastic parameters (e.g. health related ones). Also, to make this study more comprehensive one might address additional failure scenarios. Direct leakage into the aquifer and spillage of chemicals and returned fluid from the surface can be proposed as extensions of the current study. The results of the other failure events can then be compared to the present work and hopefully lead to a better understanding of the magnitude and significance of each scenario, and possibly define a set of guidelines on how better to manage the risk during such episodes. v Acknowledgements I would like to express my deepest appreciation and thanks to my advisor Professor Fred Aminzadeh for his never-ending support, encouragement and advice throughout my studies at the University of Southern California. He continually and persuasively conveyed a spirit of adventure in regard to research and scholarship. I would especially like to thank my co-advisor Professor Felipe de Barros for his tremendous help in starting and completing this research work. This study would not have been possible without his continuous support. I would also like to thank Professor Roger Ghanem, Professor Najmedin Meshkati, Professor Kelly Sanders, and Professor Lucio Soibelman for their advice, recommendations, and support. In addition, I would like to express my appreciation to Professor Behnam Jafarpour and Professor Iraj Ershaghi for the constructive discussions during my studies. I would also like to thank Induced Seismicity Consortium at USC for all the support and help throughout this research work. Special and sincere thanks go to my friends and colleagues at Sonny Astani department of Civil and Environemntal Engineering and also the Mork Family Department of Chemical Engineering and Material Science, those who incented me to strive towards my goal. Especially, I would like to express my sincere gratitude to Dr. Farrokh Jazizadeh and Mr. Amirhossein Eftekharian for their friendship and supports. Last but not least, a special thanks to my family. Words cannot express how grateful I am to my mother, father, and brother for their unconditional love and support and for all of the sacrifices that they have made on my behalf. vi Table of Contents Abstract….. ..................................................................................................................................... ii Acknowledgements ......................................................................................................................... v Table of Contents ........................................................................................................................... vi List of Tables ................................................................................................................................. ix List of Figures ................................................................................................................................. x Introduction ................................................................................................................. 1 Chapter 1: 1.1 Background ................................................................................................................. 1 1.2 Motivation ................................................................................................................... 3 1.3 Research Objectives and Problem Statement .............................................................. 4 1.4 Organization of the Dissertation .................................................................................. 5 Literature Review ........................................................................................................ 6 Chapter 2: 2.1 Introduction ................................................................................................................. 6 2.2 Natural Gas Production ............................................................................................... 6 2.3 Introduction to Hydraulic Fracturing .......................................................................... 8 2.3.1 Operation Steps .............................................................................................................. 8 2.3.2 Fracturing Slurry Constituents ..................................................................................... 11 2.4 Hydraulic Fracturing Environmental Impact ............................................................ 13 2.4.1 Air ................................................................................................................................ 14 2.4.2 Induced Seismicity ....................................................................................................... 14 vii 2.4.3 Water and Wastewater ................................................................................................. 15 2.5 Risk Pathways of Fluid Flow to the Shallow Aquifers ............................................. 18 2.6 Application of Numerical Simulation of Fluid Movement in Hydraulic Fracturing 23 2.7 Health Risk Assessment from Contaminated Subsurface Resources ........................ 27 Well Casing Failure and Upward Flow Movement ................................................... 29 Chapter 3: 3.1 Introduction ............................................................................................................... 29 3.2 Problem Statement .................................................................................................... 29 3.2.1 Conceptual Hydrogeological Model ............................................................................ 31 3.2.2 Formulation of Risk Model .......................................................................................... 32 3.3 Modeling and Numerical Simulation ........................................................................ 35 3.3.1 Hydrogeological Parameters ........................................................................................ 35 3.3.2 Flow and Transport Model Formulation ...................................................................... 36 3.3.3 Numerical Implementation .......................................................................................... 43 3.4 Simulation Scenarios ................................................................................................. 50 3.5 Result and Discussion ............................................................................................... 52 3.5.1 Chemical Transport Mechanisms................................................................................. 53 3.5.2 Defining the Threshold and Health Parameters for the Chemical of Interest .............. 54 3.5.3 Main Case Result ......................................................................................................... 55 3.5.4 Effect of Changes in Leakage Rate (Scenario Set I) .................................................... 72 3.5.5 Effect of Changes in Leakage Depth (Scenario Set II) ................................................ 73 3.5.6 Effect of Changes in Absolute Permeability (Scenario Set III) ................................... 77 viii 3.5.7 Effect of Changes in Porosity (Scenario Set IV) ......................................................... 78 3.5.8 Deterministic Risk Quantification for the Most Critical Scenarios ............................. 79 3.5.9 Parametric Uncertainty Analysis and the Stochastic Characterization of Risk ........... 80 3.5.10 Effect of Spatial Uncertainty (Heterogeneity) of the Media on the Characterization of Risk… ................................................................................................................................ 88 Conclusion ................................................................................................................. 97 Chapter 4: Future Research Direction ....................................................................................... 102 Chapter 5: Appendix A: Defining the Normalized Threshold for the Chemical of Interest (TEPA) ........... 104 Appendix B: Calculation of Reference Dose for the Chemical of Interest (TEPA) ................... 106 Appendix C: Water Contamination Mapping Tool .................................................................... 107 Appendix D: List of Publications ............................................................................................... 119 Bibliography ............................................................................................................................... 121 ix List of Tables Table 2-1: Well-known gas shale formations in lower 48 states .................................................... 7 Table 2-2: Example constituents of a fracturing slurry ................................................................ 11 Table 2-3: List of the hazardous chemicals used in three fracturing jobs in Germany .............. 16 Table 2-4: Comparison of the numerical modellings in hydraulic fracturing .............................. 26 Table 3-1: ECLIPSE100 equations used....................................................................................... 41 Table 3-2: Different grid cells used in the simulation .................................................................. 44 Table 3-3: Schedule of the fracturing operation ........................................................................... 50 Table 3-4: Parameters used in the main case model ..................................................................... 51 Table 3-5: Parameters used for different deterministic scenarios ................................................. 52 Table 3-6: Health-related parameters used in the model .............................................................. 55 Table 3-7: Statistical distributions for uncertain parameters ........................................................ 83 Table 3-8: Steps of the sequential Gaussian simulation algorithm ............................................... 89 Table A-1: Tetraethylenepentamine concentrations in two fracturing fluids ............................. 105 Table A-2: Calculation of threshold concentration ..................................................................... 105 x List of Figures Figure 1-1: Hydraulic fracturing operation ..................................................................................... 2 Figure 2-1: Major shale plays in lower 48 states ............................................................................ 6 Figure 2-2: The first hydraulic fracturing project. .......................................................................... 8 Figure 2-3: Schematic drawing of well casing and cementing around it. ....................................... 9 Figure 2-4: Sample plot of pressure, injection rate, and proppant concentration. ........................ 10 Figure 2-5: Upward flow migration scenarios proposed by the USEPA. ..................................... 20 Figure 3-1: Schematic representation of the failure event. ........................................................... 32 Figure 3-2: Numerical model layers. ............................................................................................ 45 Figure 3-3: Simulation grid. .......................................................................................................... 49 Figure 3-4: Analytical aquifers on sides of the model. ................................................................. 50 Figure 3-5: Cross-section of the main case model in various time-steps ..................................... 58 Figure 3-6: Concentration plume expansion ................................................................................. 59 Figure 3-7: Temporal evolution of the concentration at the locations of interest ......................... 59 Figure 3-8: Contour plots of normalized concentration ................................................................ 60 Figure 3-9: Concentration-time curve general trend for the case with 3 . 0   w S and 1 . 0   ......... 61 Figure 3-10: Effect of permeability ratio on concentration magnitude ........................................ 63 Figure 3-11: Effect of permeability ratio on concentration plume geometry ............................... 66 Figure 3-12: Effect of permeability ratio on concentration plume geometry (plan view). ........... 67 Figure 3-13: Effect of initial water saturation on concentration magnitude ................................. 69 Figure 3-14: Effect of initial water saturation on concentration plume geometry ( 1 . 0   ) ......... 70 Figure 3-15: Effect of impervious layer on diffusion phenomenon ............................................. 71 Figure 3-16: Impact of  and  w S on the Hazard Index (HI) ....................................................... 71 xi Figure 3-17: Effect of the leakage rate (Q) on concentration values ............................................ 73 Figure 3-18: Effect of the leakage depth (H) on concentration values ......................................... 74 Figure 3-19: Breakthrough curves for leakage depths 5 (a) and 10 meters (b) ............................ 76 Figure 3-20: Breakthrough curves for various leakage depth scenarios ....................................... 77 Figure 3-21: Effect of absolute horizontal permeability (k) on concentration values .................. 78 Figure 3-22: Effect of the sand layer porosity (φ) on concentration values ................................. 79 Figure 3-23: Changes of hazard index (HI) with leakage depth (H) ............................................ 80 Figure 3-24: Schematic representation of sources of uncertainty in the simulation model .......... 81 Figure 3-25: Workflow of the random simulations - Coupling MATLAB and ECLIPSE ........... 82 Figure 3-26: Temporal evolution of the probability of exceedance  ......................................... 84 Figure 3-27: Histogram of concentration at the monitoring well ................................................. 85 Figure 3-28: CDF of concentration at the monitoring well .......................................................... 86 Figure 3-29: CDF of the arrival time for the threshold concentration .......................................... 87 Figure 3-30: Hazard index (HI) CDF............................................................................................ 88 Figure 3-31: Sample realizations for the sand layer permeability field ........................................ 90 Figure 3-32: Concentration CDF for sand layer variogram sill of 0.4.......................................... 91 Figure 3-33: Source-to-receptor arrival time variogram sill of 0.4 ............................................. 91 Figure 3-34: Hazard Index CDF for sand layer variogram sill of 0.4........................................... 92 Figure 3-35: Sample realizations for the aquifer permeability field ............................................. 93 Figure 3-36: Concentration CDFs of two extreme values of aquifer variogram sill .................... 94 Figure 3-37: Arrival time of two extreme values of aquifer variogram sill.................................. 94 Figure 3-38: Mean values of source-to-receptor arrival time ....................................................... 95 Figure 3-39: Hazard Index of two extreme values of aquifer variogram sill................................ 95 xii Figure 3-40: Mean values of Hazard Index for different cases of variogram sill ......................... 96 Figure A-1: Hazardous Material Identification Rating for Tetraethylenepentamine .................. 104 Figure C-1: Attribute categorization for water contamination mapping tool (WCMT) ............. 109 Figure C-2: WCMT interface ...................................................................................................... 110 Figure C-3: Selecting district of interest and the name of county .............................................. 110 Figure C-4: Selecting year of interest ......................................................................................... 111 Figure C-5: Selecting type of data to be displayed ..................................................................... 111 Figure C-6: Selecting data features ............................................................................................. 112 Figure C-7: Software output – volume and spatial distribution of injected water ...................... 112 Figure C-8: Software output – Tabulated results for injected water volume in Kern County .... 113 Figure C-9: Displaying TDS and hydraulic fracturing data at the same time ............................ 113 Figure C-10: Showing numerical values of data on the map ...................................................... 114 Figure C-11: TDS and injected volume data in Kern County in year 2013 ............................... 114 Figure C-12: TDS and injected volume data in Los Angeles County in year 2012 ................... 115 Figure C-13: TDS and injected volume data in Ventura County in year 2012........................... 116 Figure C-14: TDS and injected volume data in Glenn County in year 2012 .............................. 116 Figure C-15: TDS and injected volume data in Colusa County in year 2012 ............................ 117 Figure C-16: TDS and injected volume data in Sutter County in year 2012 .............................. 117 Page | 1 Introduction Chapter 1: 1.1 Background Population growth and the vital need for energy are among the important challenges to be tackled by nations across the globe. Fossil fuels have long been relied upon to meet the energy demand, despite their limited quantity and causing environmental damage. Greenhouse gas emission is broadly recognized as an example of a worldwide issue which is encountered when heavy hydrocarbon fuels are burnt. To move towards a more sustainable world, one might think of renewable sources of energy as the best alternative to conventional fuels. But renewables are not yet ready to go online due to their highly intermittent nature as well as some availability limitations. More extensive research and investigation is needed before renewables will be fully able to quench the world’s thirst for energy. Particularly in the United States the need to lessen dependency on imported conventional energy sources has always been a point of discussion. According to the United States Energy Information Administration (USEIA) the country is trending downward in the use of imported oil, after the peak year of 2005 (USEIA 2011). It should be noted that shifting from conventional sources to the more sustainable renewables requires substantial time and effort, primarily to ensure that the energy production from renewables is economically viable (Spellman 2012). Natural gas, as a candidate for short and even mid-term energy solutions, has garnered the attention of both energy companies and energy demanding countries. As a light hydrocarbon fuel, natural gas is easily accessible and is dependable, which means that unlike some renewable energy sources, natural gas’s performance is not a function of environmental conditions such as sunlight or wind, nor is it intermittent. When compared to the heavier hydrocarbon fuels, natural Page | 2 gas emits less harmful by-products when burnt. Natural gas produces 50% and 30% lower CO 2 emissions respectively when compared to coal and fuel oil (Spellman 2012; Howarth et al. 2011). A well-known technique for extracting natural gas is hydraulic fracturing (also known as hydrofracking). Hydraulic fracturing is an old reservoir stimulation method that has been acknowledged in the petroleum industry for the past fifty years (AAPG 2014). The technology aims to increase oil and gas production through a sequence of processes, which include drilling vertical or/and horizontal wells, injection of a highly pressurized fluid, cracking the matrix formation, and creating fissures which in turn increase permeability values of the reservoir (King 2012). Figure 1-1 pictures a typical hydraulic fracturing operation. Along with the horizontal drilling technology, hydraulic fracturing has recently become very popular among the oil and gas industry companies, primarily because of the applications of both technologies in producing natural gas from shale and other tight formations (Arthur et al. 2009; USEPA 2012). Figure 1-1: Hydraulic fracturing operation (Probublica 2014) Page | 3 As stated by Daneshy, hydraulic fracturing changes the flow pattern regime toward the well from radial to linear (Daneshy 1990), at the same time horizontal wellbore helps to obtain much larger reservoir to well contact and production is heightened accordingly. Although hydraulic fracturing is a major contributor to increasing natural gas production and maintaining a relatively low price, it is suspected to be responsible for environmental downsides ranging from consuming and polluting natural resources (i.e. water and air) to manipulating the regional faults beneath the ground and inducing earthquake tremors. Upward fluid migration (e.g. reservoir brine, fracturing fluid, or methane) and contamination of shallower aquifers are potential issues attributed to hydraulic fracturing (USEPA 2012; USEPA 2015), particularly in fracturing shallower reservoirs. The fracturing fluid consists of approximately 90 percent water, 9.5 percent propping agent and 0.5 percent of chemical additives (Holloway & Rudd 2013) including but not limited to, polymers, cross- linkers, biocides, surfactants, buffers, fluid loss agents, stabilizers, friction reducers, and corrosion and scaling inhibitors (USEPA 2004; Gregory et al. 2011). 1.2 Motivation Among the various possibilities for a negative environmental impact from hydraulic fracturing, groundwater contamination is the topic of this study because of a high risk to human health once the contaminant is found in drinking water. It is imperative to ascertain which combinations of operational and hydrogeological parameters could be responsible for polluting a fresh shallow aquifer. If, in a hypothetical situation, an aquifer becomes polluted, what would be the measure of risk to human health for a particular chemical compound or a mixture of chemicals that had Page | 4 been used in the fracturing process? Knowing the possibility of groundwater contamination and the extent of pollution will be helpful in defining remediation scenarios. Despite the existence of several works that reference groundwater contamination, some novel risk frameworks need to be developed which are tailored to the hydraulic fracturing process, while capturing the specific features of the subsurface (e.g. coupled saturated- unsaturated flow) on the risk response. It is essential that contaminant site managers quantify the uncertainty of various parameters, such as the location and intensity of the contamination source, anisotropy and heterogeneity of the subsurface permeability fields, and the water content of the shale and overburden formations. The impetus for this research is the need to develop a risk- based framework which links the hydrogeological and hydraulic fracturing operational parameters to the groundwater contamination and risk metrics. The ultimate goal is to ensure that, while maintaining the efficiency of the operation, water resources and human health are protected. 1.3 Research Objectives and Problem Statement This research is focused on implementing a systematic and scientific approach to track chemical species transportation inside a porous medium towards the aquifer, and to define associated arrival times and concentrations. One specific risk pathway (leakage through the casing while fracturing) has been addressed, leaving the following questions to be answered: 1. Under what situations is fracturing fluid able to reach the aquifer? What are the concentrations and arrival times associated with such an event? 2. What factors are most important in the simulation, and how sensitive are the results to those different parameters? Page | 5 3. What are the most important issues when characterizing the risk in a casing failure scenario? 4. What is the critical depth for the leakage point along the casing for a particular chemical species used in the fracturing slurry? 5. What is the contamination probability at different times after the leakage took place? How likely is it to exceed the thresholds set by drinking water standards and guidelines? 6. Which concern is of overriding importance for the regulators to reduce the human health risk in such complex systems? Hypothetical scenarios are designed and simulated correspondingly, with each one of them covering a set of strictures borrowed from real cases. Answers to the problems are broadly based on the range of hydrogeological and operational parameters, and the assumptions presumed in numerical models. Nevertheless, parameters leading to extreme case scenarios are in sharper focus in this study. Events with catastrophic consequences result in great risk, even if the likelihood of such incidents is low. Assessing the probability of aquifer contamination and crossing the threshold under different stochastic scenarios is the next goal of this study. 1.4 Organization of the Dissertation The remainder of this manuscript is organized as follows. Chapter two includes a brief literature review on hydraulic fracturing steps and the environmental impacts with the focus on groundwater contamination. Chemical possible groundwater contamination, and results of sensitivity analysis on parameters domain and the environmental performance metrics (EPMs) are discussed in Chapter three. Chapter three further includes effects of heterogeneity of the permeability field within the aquifer and sand medium on the risk measures. Chapter four concludes the study and finally Chapter five draws future research directions. Page | 6 Literature Review Chapter 2: 2.1 Introduction This chapter starts with natural gas production, challenges faced by the industry, and major shale plays in the United States. It then discusses the fundamentals of hydraulic fracturing as a stimulation method widely used in tight gas production. Next, the environmental impact of fracturing operations are addressed, as well as the risk pathways to the groundwater aquifer. The chapter is completed by introducing the application of numerical methods in subsurface fluid movement. 2.2 Natural Gas Production U.S. natural gas production is growing quickly and it is predicted that it will continue its growth with a rate of 1.3 percent per year until 2040, so that by 2020 the country will become an exporter of natural gas. According to the EIA, sources of natural gas production are: shale gas, tight gas, coal-bed methane, onshore conventional and offshore -- among which shale gas is the main contributor to the production. Major oil and gas shale plays in the lower 48 states are pictured in Figure 2-1. Figure 2-1: Major shale plays in lower 48 states (USEIA 2015) Page | 7 Shale gas formed 34% of the total natural gas production of the U.S. in 2011 and the number is expected to rise to 50% in 2040 (USEIA 2013). Table 2-1 lists the average properties of major gas shale plays in the United States as reported by the Energy Information Administration (USEIA 2011). Table 2-1: Well-known gas shale formations in lower 48 states separated by regions (USEIA 2011) Average Properties of the Shale Play Region Gas Shale Depth (m) Thickness (m) Porosity (%) TOC * (% wt) Northeast Marcellus 2,057 38 8 12 Devonian Big Sandy 1,158 53 10 3.75 Devonian Low Thermal Maturity 914 113 7 N.A. Devonian Greater Siltstone 887 190 5.8 N.A. New Albany 838 61 12 13 Antrim 427 29 9 11 Gulf Coast Haynesville 3,658 76 8.5 2.25 Eagle Ford 2,134 61 9 4.25 Floyd-Neal & Conasauga 2,438 40 1.6 1.8 Mid- Continent Fayetteville 1,219 34 5 6.9 Central Woodford 1,524 76 6 4 Western Woodford 2,896 46 7 6.5 Cana Woodford 4,115 61 7 6 Southwest Barnett 2,286 91 5 N.A. Barnett Woodford 3,109 122 N.A. 5.5 Avalon/Bone Springs 2,667 396 N.A. N.A. Rocky Mountain Hilliard-Baxter-Mancos 4,496 937 4.25 1.75 Lewis 1,372 76 3.5 N.A. Mancos 4,648 914 3.5 14 Shale is a low permeable sedimentary rock which consists of clay layers and serves as the source and reservoir rock for the dry gas (Arthur et al. 2009). Due to the impervious characteristic of the rock, gas does not readily flow through the tight shale medium. Advanced technologies such as hydraulic fracturing are then implemented to facilitate the production from the formation (King 2012) and make the gas production economically feasible (USEIA 2011). * Total Organic Content in weight percent Page | 8 2.3 Introduction to Hydraulic Fracturing Hydraulic fracturing has long been utilized in the oil and gas industry (Economides & Martin 2007; Adams & Rowe 2013), however this techniques is not limited to oil and gas industry and has other applications such as enhancing water well production, mining, determination of rock stress, geothermal operations, carbon sequestrations, and mitigation of ruck burst (Adams & Rowe 2013). “Fracking” and “hydrofracking” are other terms that are frequently used to describe the procedure (King 2012). The first hydraulic fracturing operation was carried out on Klepper gas unit 1 of the Hugoton gas field in 1947 (Figure 2-2). The field was located in western Kansas (Miskimins et al. 2006). Figure 2-2: The first hydraulic fracturing project in the Klepper gas unit number 1 in 1947 (AAPG 2014). 2.3.1 Operation Steps Hydraulically fracturing rocks is a complementary step in drilling, casing, cementing, and injection. Spellman describes the general steps of the operation as follows (Spellman 2012): - Drilling the wellbore to a certain depth. - Placing the casing in the center of borehole allowing an annular empty space separate the casing and surrounding formation. In most recent cases and as suggested by the Page | 9 United States Environmental Protection Agency (USEPA), an intermediate casing is placed in the well which is extended from surface to depths below fresh aquifer. - Completing the well by pumping cement through the casing and up the annular space so that the casing would not be in direct contact with the formation. The main reason is to save casing from corrosive environment in adjacent formation and also to protect groundwater table. - Repeating the steps above until reaching the depth of interest. - Inserting a pipe inside the casing (termed as tubing) and injecting high pressure fluid through the tubing to fracture the formation. The fluid is a mixture of water, tiny particles known as proppants, and chemical additives. Figure 2-3: Schematic drawing of well casing and cementing around it (Spellman 2012). Viscosity of the fracturing fluid is such that it requires high pressure pumps to inject the fluid down the well. A pump used in a fracturing operation usually has 700 to 2700 hydraulic horsepower (Economides & Martin 2007). In addition, injection rate of the fracturing operation is a parameter that must be chosen according to the formation of target and other geological and Page | 10 physical properties of the area. In general the value can be anywhere between 15 to 100 bbl/min or 144 to 960 m 3 /hr (Spellman 2012). An injection process typically consists of two stages. First is the pad stage which includes pumping fluid into the formation, followed by stage two, which is pumping the fluid and proppants mixture (aka slurry or proppant-laden fluid) into the formation. High pressure is needed to crack the target formation, so the fracturing interval (as long as 1500 m in some cases) is divided into smaller sections. Each section is then fractured separately in a stage which may take from 45 minutes to two hours. After the fracturing job is complete the fluid is pumped out of the well. The role of proppant (e.g. natural sand) is to withstand the closer pressure of the fracture – which is approximately equal to the minimum principal stress of the matrix formation- and to prop the created fracture open once the fluid has been taken out. A typical job plot of a fracturing process is shown in Figure 2-4 where STP refers to the surface treating pressure or the well head pressure, and BHTP means bottom hole treating pressure or the pressure needed to crack the formation (Economides & Martin 2007). Figure 2-4: Sample plot of pressure, injection rate, and proppant concentration of a frack job (Economides & Martin 2007). Page | 11 2.3.2 Fracturing Slurry Constituents A suitable fracturing fluid is the one which is inexpensive, does not need large volumes of water to be added, flows well with minimal friction, creates large fissures, is able to carry proppants effectively and in high concentrations and then obtain low viscosity when returning the fluid (Wiseman 2009). Fracturing slurry is a mixture of water (90%), propping agent(8~9.5%), and chemical additives (0.5~2%). Organic and inorganic chemicals are used in order to improve the efficiency of the operation. Table 2-2 shows common constituents of a fracturing slurry (Holloway & Rudd 2013). Table 2-2: Example constituents of a fracturing slurry (Holloway & Rudd 2013) Name Example Amount frequently used Fluid Water 98-99% of total vol. Acid Hydrochloric Acid 5-25 % in solution Corrosion Inhibitor Ammonium Bisulfate 0.2-0.5% of acid vol. Biocides Sodium Hypochlorite 0.005%-0.05% of total vol. Friction Reducers Polyacrylamide based compounds 0.025% of total vol. Gelling Agents Guar Gum Not often used Crosslinking Agent Boric Acide – Zirconium - Breaker Solution - - Oxygen Scavenger Ammonium Bisulfate - Surfactant - 0.5 to 2 gal/1000gal of total vol. Scale Inhibitor Ethylen Glycol Seldom used Proppants Sand 1-1.9 % of total vol. Fracturing fluid systems are divided into three groups -- slickwater, linear gel (gel used as a cross-linker), and a combination of both, with the addition of the proppants (Palisch et al. 2008). The difference between the two systems concerns the viscosity; in slickwater the addition of a friction reducer makes the fluid smoother and enables higher pumping rates when compared to linear gel (Palisch et al. 2008). As stated by Frackfocus, the percentage of chemical additives in a slickwater fluid ranges from 0.5 to 2% (FrackFocus, 2014). Page | 12 2.3.3 Practical Considerations: Conventional vs. Unconventional When designing a fracturing job, one should take a series of factors into consideration. Type and tightness of formation, for instance, are two of the most important aspects as the choice of fracturing technology is greatly influenced by the geological conditions of the oil and gas resources (Wiseman 2009). As mentioned earlier, hydraulic fracturing has been typically used as a well stimulation method in conventional production when the permeability number of the reservoir has declined. In this case, fracturing is perceived as a technology to remove near well- bore damages and increase the permeability value, whereas in the case of unconventional production fracturing is the key technology and without its presence the production will never start. Some of the key differences between fracturing applications in conventional versus unconventional operations are outlined below: - In conventional, a small fraction of wells usually undergo vertical fracturing operation, whereas most of the drilled unconventional wells are accompanied by horizontal drilling and fracturing as the final completion step prior to production inception (USEPA 2015). - Conventional wells are commonly shallow and vertical (i.e. 1,200~1,600 m), however in unconventional fracturing operators deal with deep target formations and long horizontal intervals (4~5 km) (USEPA 2015). Deeper formations and long intervals require higher injection pressures to be provided on the surface which in turn raise the risk of fracturing operation as a whole. - Unconventional jobs are quite water-intensive as opposed to conventional ones. - Proppants and proppant-laden fluids are used more extensively in fracturing tight formations. Use of proppants requires addition of new chemicals to give the fluid proppant carriage ability. Friction reducers and breakers are good examples of these Page | 13 fluids. In conventional fracturing treatment, however, proppants are used less frequesntly. In acid treatment of carbonate reservoirs, as an example, fractures are created by dissolving the formation matrix. - Production regime in conventional operations is a bell-shaped curve and it takes a considerable time for the production to drop. In unconventional operations, however, production peaks and drops quite rapidly which necessitates refracturing the formation in a year or two from the first fracturing job. As a result, the frequency of the fracturing operations in unconventional productions will be even higher. Given discussed key differences, hydraulic fracturing in unconventional resources over the last decade has gained a lot of attention from scientists and also public, as it potentially increases the load on environment and natural resources. 2.4 Hydraulic Fracturing Environmental Impact Despite the benefits of hydraulic fracturing in enhancing gas production and stimulating the economy, environmentalists continue to question its use (Kargbo et al. 2010; USEPA 2012; USEPA 2015). The questions remain unanswered because of the relative youth of the technology in gas production, and the unrevealed long-term effects on the environment (Gassiat et al. 2013). It is now critical to address these unique safety and environmental considerations (Heinecke et al. 2014; Jabbari, Ashayeri, et al. 2015). In broad strokes, some drawbacks can be mitigated by setting stricter regulations and improving the design. Potential environmental impacts can be classified into three major categories: air, induced seismicity, and water and wastewater. Page | 14 2.4.1 Air Air pollution in general is linked to natural gas development activities, including fracturing. For example, a human health risk assessment study by McKenzie et al. (2012) on a gas field in Garfield County, Colorado suggested that people living at a distance greater than half a mile from the site were in less danger than the ones living at a lesser distance. One chemical of concern that is used extensively in drilling and fracturing is diesel -- which is helpful in reducing the friction. Replacing diesel by natural gas in drilling can decline emission of VOC’s into the air by up to 85%. Methane, carbon dioxide and other VOC’s are causing major air pollution in areas with gas production activities in Wyoming, Colorado and Texas (Kargbo et al. 2010). According to a more recent study by Moore et al. (2014), particulate matters smaller than 2.5 μm in size (PM 2.5 ) and No x are found in diesel emissions during preproduction (i.e. drilling and hydraulic fracturing activities) of unconventional natural gas production. 2.4.2 Induced Seismicity Subsurface injection operations normally generate low magnitude (smaller than 2 in Richter scale) seismic events which are termed as micro-earthquakes or microseismic (Ellsworth 2013). In a few cases of hydrofracking jobs, seismic events have been felt by nearby residents and the operator has been forced to stop the injection. In 2011, two seismic events (2.3 and 1.5 magnitude in Richter scale) were recorded close to Blackpool, Lancashire in the U.K. -- the operator had to halt the fracturing operation in nearby Bowland Shale formation (Pater & Baisch 2011). Generally, the risk of generating serious earthquakes as a result of hydraulic fracturing is low when compared to the deep-well injection process, which shows higher probabilities of observed larger seismic events (Ellsworth 2013; Goebel et al. 2014; Aminzadeh et al. 2014). Historically, the oldest injection-induced seismic events are those of the Rocky Mountain Page | 15 Arsenal waste injection site in Denver, Colorado (Healy & Rubey 1968). The most recent events were generated in Youngstown, Ohio (ODNR 2012) and central Oklahoma (Holland 2011), both in 2011. The former is also known to be the largest injection-induced event with a 5.6 magnitude (Ellsworth 2013). 2.4.3 Water and Wastewater In a typical fracturing job for a single well, between 2 and 5 million gallons of fresh water, either from surface or ground sources, are being consumed (USEPA 2010), which is a considerable amount and can raise water conservation issues in areas with an arid climate. In addition, chemical additives to the fracturing slurry are various, and range from benign (e.g. guar gum) to toxic (e.g. Tetramethylammonium chloride) to very toxic materials (e.g. Kathon ® which is a biocide) (Gordalla et al. 2013). Injecting those chemicals underground is clearly increasing the risk of having ground or surface water bodies contaminated. Table 2-3 elaborates on hazardous chemical concentrations used in three different fracturing operations in Germany named as BT12, CZ3a, and D3. The allowable health reference concentration for drinking water is also tabulated. BT12 and CZ3a are operations carried out using gel-based fluid with chemical percents of 3 and 1.1 respectively of the total fluid without proppants. In contrast, D3 is using slickwater fluid and showing a chemical concentration of 0.16% of the fluid (Gordalla et al. 2013). Page | 16 Table 2-3: List of the hazardous chemicals used in three fracturing jobs in Germany (Gordalla et al. 2013) Chemical BT12 CZ3a D3 Reference Concentration (mg/l) Applied Concentration (mg/l) Amphoteric alkyl amines - 193 - 0.0003 2-Butoxyethanol 5,440 4,730 - 0.35 Citrus terpenes 184 - - 0.0003 Hydrotreated light petroleum distillates - - 208 0.0003 Inorganic borate 6 2.9 - 1 Kathon ® (1:3 combination of the followings) 5.5 4.4 3.6 0.0005 2-Methyl-2H-isothiazol-3-one 1.4 1.1 0.9 0.0001 5-Chloro-2-methyl-2H-isothiazol-3- one 4.1 3.3 2.7 0.0001 Methanol - 705 - 1.75 Propan-2-ol 184 141 - 8.4 Sodium bromate 348 218 - 0.01 Tetraethylenepentamine 736 173 - 0.0003 Tetramethylammonium chloride - 711 502 0.0003 Zirconium dichloride oxide 55 - - 0.0003 The produced wastewater not only includes the benign or toxic chemicals frequently used in the fracturing fluid, but it also bears natural salts, heavy metals (e.g. mercury and lead), naturally occurring radioactive materials (e.g. uranium and radium), gases (e.g. natural gas and hydrogen sulfide), volatile organic compounds (VOCs), bacteria, and pathogens from the deep under the ground (USEPA 2012; Gregory et al. 2011; Gordalla et al. 2013). In a study conducted by Hammer & VanBriesen (2012), storage tanks and ponds on a fracturing site, spatially extensive impacts of direct release of produced water into the environment, waste injection in deep wells, effluents to surface water, impact of land application such as dust and storm water runoff from the site, and management of residual wastes are introduced as potential impacts of shale gas wastewater to the water resources. Advanced treatment processes must be employed to treat the wastewater if it is meant to be discharged into natural streams because conventional wastewater treatment plants are simply incapable of removing all the compounds commonly found in the returned fluid. In addition, Page | 17 extreme caution must be used when storing the returned fluid on-site as well as during off-site transport. Another common method for handling wastewater, besides the treatment, is underground injection into deep formations through a class II well (USEPA 2010). Deep wastewater injection is an arguable technology and has its own requirements, considerations, and issues. As is briefly discussed in section 2.4.2 of this document, induced seismicity is one of the drawbacks which is thought to be a result of fluid injection down deep into the earth (Nicholson 1992). Surface water contamination is a concern in areas with gas operations (Hammer & VanBriesen 2012; VanBriesen et al. 2014; Glazer et al. 2014). If the returned fluid treatment is not handled properly (namely the wastewater treatment plant is not designed to remove specific chemicals), it may cause contamination in nearby streams. The Ohio River contamination by Strontium, Barium, and Bromides is an example of this case (Volz et al. 2011). Groundwater can also be a target of contamination from the chemicals found in returned fluid through a failure in appropriate onsite storage and management, leading to spills from the surface to the groundwater. Also gas migration from down deep, and contamination from additives of fracturing fluid during and after the high pressure injection process are among other potential risks to the groundwater (Jabbari, Aminzadeh, et al. 2015; Llewellyn et al. 2015; Osborn et al. 2011; Rozell & Reaven 2012; USEPA 2012; Vengosh et al. 2013). Moreover, there have been complaints about turbidity, color change, and the odor of water from people residing in the vicinity of oil and gas fracturing wells (in east coast cities and usually near unconventional gas production sites) (Holloway & Rudd 2013). Pavillion, Wyoming is an example of a community in close proximity to a hydraulic fracturing site where a broad range of organic Page | 18 compounds (e.g. trimethylbenzenes, phenols, 2-butoxyethanol, di-and tri-ethylene glycol, etc.) as well as elevated levels of potassium, chloride, and alkalinity were reported from deep monitoring wells (DiGiulio et al. 2011). From a human health point of view, addressing groundwater pollution is crucial, especially for regions with water scarcity and a high dependence on the groundwater tables to provide safe, potable water (Kargbo et al. 2010; USEPA 2012). In short, the possible risk of drinking water contamination as a result of hydraulic fracturing has been controversial enough to propel the residents of some American states to ask their local authorities to put moratoriums on hydraulic fracturing, which is reminiscent of the cases in Pavillion, Wyoming in 2009 and New York in 2010 (Rahm 2011). 2.5 Risk Pathways of Fluid Flow to the Shallow Aquifers In order to better analyze impacts of hydraulic fracturing on fresh groundwater basins, it is important to have a deep insight into the mechanisms the contaminants find their ways to the aquifers (risk pathways). Risk pathways to shallow aquifer can be divided into two major categories: (1) overground accidental spills from the pits used for temporary storage of the returned fluid, or from accidents while transporting fluids, or stemming from drilling activities and (2) underground failure and upward fluid migration (Rozell & Reaven 2012). Upward fluid migration to the aquifer can occur through different pathways. The USEPA suggests the following possibilities (Figure 2-5): (1) Damages in well integrity as a consequence of corrosion or an improper cementing job. Leakage of fracturing fluid into the formation is a possible outcome of this scenario. (2) Poor design of the fracturing job, causing fracture growth in overburden formations and providing hydraulic connections between the shale gas reservoir and aquifers above. In some Page | 19 cases and in the presence of hydrocarbon zones above the shale layer, the situation becomes more critical as the fluid is capable of carrying hydrocarbon while travelling to the shallower depths. (3) Fault reactivation resulting from the fracturing process and providing permeable pathways for the fluid. (4) Fracturing the overburden formation and creating pathways into wells in the vicinity of the fracturing well (abandoned wells in most cases). In this situation, fluid might find its way through the aquifer using the adjacent wells and through issues in well integrity (USEPA 2012). Well integrity issues described in Scenario 1 are further elaborated as follows (Holloway & Rudd 2013): - Behind the casing, movement of the fluid toward the well head in which fluid finds the space between cementing and casing as a possible pathway (annular flow) - Leakage from the well in a radial pattern so that the fluid can reach the nearby formation and the chemicals can undergo different transport mechanisms (leak flow). Page | 20 Figure 2-5: Upward flow migration scenarios proposed by the USEPA (2012). King (2012) has also named surface casing rupture as one of the potential leakage scenarios in a fracturing job. Real-time pressure monitoring helps to increase the discovery of a leak – in the same way that the pressure drop following a major leakage in a pipe enables the operators to discover the leakage in 5 minutes (Gordalla et al. 2013). However, when a smaller leak is taking place, there is a possibility that failure remains undetected during the hydraulic fracturing Scenario 1 Scenario 2 Scenario 3 Scenario 4 Page | 21 operation (Gordalla et al. 2013). It is, therefore, important to closely investigate the aftermath of all leakage events, particularly when the rupture is within or close to the aquifer. In a probability-bound analysis study from Rozell & Reaven (2012) on the risk of water pollution in Marcellus shale, it was concluded that among other risk pathways such as spills from on-site activities and wastewater disposal, the well casing failure was an important contributor to the contamination risk where it shows lower than 0.01 m 3 and 9 m 3 contamination volumes for best-case 50 th percentile and worst-case 50 th percentile respectively. A study by Browning & Smith (1993) on deep well injection in Louisiana, Michigan, Nebraska, and Pennsylvania suggests malfunctioning casings and tubings as two reasons for failing mechanical integrity tests. According to the same study, out of 10,000 wells that had undergone the integrity test, 50% of them fail to pass, and casing failure was the primary reason, as it was responsible for between 45 and 85% of the failures. Among the failed tests due to casing problems, 22% showed failure in the only protective layer, so that the injection waste was able to leak into the near formation. Although the data from this study might not be directly applicable to fracturing wells because of the difference between waste injection and fracturing fluid injection operations, it still gives us insight into possible failures of casings in a fracturing well. Reservoir fluid consists of liquid (e.g. brine or mixture of brine and the fracturing fluid) and gas (e.g. methane). Gas migration to the groundwater aquifer and surface has been acknowledged in petroleum literature (Leythaeuser et al. 1982; Schmitz & Cook 1996; Nunn & Meulbroek 2002). Increased risks of explosion and groundwater contamination are two probable negatives of upward gas migration. Despite the fact that, in general, transport arrival time to the Page | 22 fresh aquifer is a function of hydraulic, physical, and geometrical parameters of the region of interest, it is apparent that when compared to liquid, it is more likely for gas to reach the layers above the reservoir, due to its physical properties and the effects of buoyancy (Nunn & Meulbroek 2002). It also takes less time for gas migration to take place, as opposed to the upward movement of fluid, as is suggested by numerical simulation results (Kissinger et al. 2013; Gassiat et al. 2013). Upward fluid migration from deep formations is extremely uncertain. The absence of a hydraulic connectivity between aquifer and shale formation, due to geometrical features and physical barriers, makes it unlikely that reservoir fluid will migrate upward, as stated by several authors. To exemplify, the separation between shallow aquifer and the shale formation is too large to make the upward fluid migration into the aquifer happen in a short period of time (Arthur et al. 2009; Zoback et al. 2010; Howarth et al. 2011). Furthermore, shale layers are usually capped with tight formations acting as natural barriers (Arthur et al. 2008). Artificially created fractures may penetrate into the overburden formation; however, the distance up to which they can propagate upward is limited, and the fracture is very unlikely to pass through the overburden, reach the aquifer and connect it to the shale layer (Zoback et al. 2010; Fisher & Warpinski 2012; Flewelling et al. 2013). In addition, when production starts, the fluid flows toward the horizontal section of the well so that any upward migration would be less probable (Howarth et al. 2011; Kissinger et al. 2013). On the other hand, hydraulic fracturing operations near conductive faults might change the outcome. As declared by Gassiat et al. (2013), under some circumstances reservoir fluids could reach the aquifer through the permeable fault, although due to the slow rate of movement, that mechanism would take place over thousands of years. Page | 23 2.6 Application of Numerical Simulation of Fluid Movement in Porous Media in Hydraulic Fracturing Environmental Impact Assessment Numerical simulations have been used to model fluid movement under the ground during and after a hydraulic fracturing job. Myers (2012) addresses the problem in five MODFLOW simulation scenarios with and without fault presence, and concludes that upward migration is possible under certain circumstances while hydraulically stimulating Marcellus shale. However, because of the simplistic model and unrealistic assumptions, the results are dubious. First, the model proposed by Myers is a single-phase flow and is not capable of honoring the physical effects pertaining to gas-water interaction (e.g. imbibition phenomenon) which are common in un-saturated Marcellus shale and its overburden. Rather, the porous medium is perceived to be full of water. Assuming a fully water-saturated medium also influences water permeability throughout the formation, as in porous media with water saturation lower than unity, the concepts of relative permeability and connate water become prominent (Saiers & Barth 2012). As implied by widely-used mathematical functions, water relative permeability is directly proportional to the water content in unsaturated matrix (Brooks & Corey 1964; van Genuchten 1980). The relative permeability calculations and modeling is totally ignored in Myers work. Additionally, neither temperature gradient from shale to aquifer, nor brine density changes are incorporated in Myers model, which decreases the reliability of simulation results, in particular when extremely temperature sensitive solubility of solutes and gases are included in the model (Saiers & Barth 2012). When explaining transport of fluid from shale formation to surface aquifers, advection is named as the responsible mechanism for fluid displacement, but in fact, the model’s boundary conditions (i.e. no-flow boundaries on four sides and constant head boundaries on bottom and top of the model) force the fluid to flow upward (Saiers & Barth 2012). The Page | 24 spatial domain of the simulation model also adds to the problem where more than two third of the space is allocated to the injection zone, and fluid movement to the sides are limited (Cohen et al. 2013). Overestimation of hydraulic conductivity numbers, neglecting anisotropy of the layered model (equivalent conductivity calculation), mathematical errors, and failing to report MODFLOW’s flaw in considering rock compressibility, are among other scientific issues related to Myers research study (Saiers & Barth 2012; Cohen et al. 2013; Reagan et al. 2015). Using a multi-component and multi-phase mass balance equation and a numerical solver, Kissinger et al. (2013) have performed a set of simulations on three hypothetical scenarios. The first scenario models fracturing fluid and brine movement induced by the high pressure of injected fluid during a short period of a fracturing operation. The first situation shows how the fluid moves in a short time period into the formation above shale. Scenario two focuses on fluid movement from deep aquifer to the shallow one through a presumed fault zone in the region. The impact of faults on fluid movement is further investigated in scenario three, in which methane migration from within the fault and going up to the surface layers is modelled. Since methane migration can start from the beginning of the gas production and last until the end, the simulation work of the third scenario is completed over long time periods. Although this research is not a quantitative risk assessment, each scenario is simulated over a range of parameters so that the uncertainty is partially included. Scenario one results in a fluid movement with a maximum vertical distance of 50 meters on top of the shale layer. In the second scenario, assigning a large pressure gradient between the deep and shallow aquifers and ascribing high permeability values for the fault, lead to a condition under which the brine and fluid reach the shallow aquifer with diluted concentrations. And finally, the combination of a highly permeable Page | 25 fault zone and an overburden with small numbers of residual gas saturation force the methane to migrate up to the surface. In 2013 Gassiat et al released findings on the possible upward migration of fluid under certain circumstances defined by sensitivity analysis on parameters such as: overburden permeability and porosity, formation compressibility, fault geometrical and geological characteristics, and hydrofracking parameters such as: fracturing length and location in shale, and duration of the fracturing job. In a two dimensional simulation work of this study, shale formation is assumed to be an over pressurized layer which is hydraulically connected to more shallow aquifers system through a fault zone. Subsurface fluid movement and solute transport have been simulated using SUTRA-MS of US Geological Survey, which is a single-phase groundwater flow package. Results show that fracturing the upper portions of the over pressurized shale formation may lead to aquifer contamination in long time periods (thousands years) (Gassiat et al. 2013). Gassiat et al. (2013), however, did not consider multi-phase flow and transport and the associated capillary effects (Reagan et al. 2015). Table ‎ 2-4 summarizes the three studies discussed previously and the simulation work of the current study and compares approaches taken to model the flow and transport problem. As seen in the table, none of the three studies propose a human health risk assessment method for the case in which water resources are being contaminated in subsurface injection operations. Page | 26 Table 2-4: Comparison of the numerical modelling of upward flow migration in hydraulic fracturing Study Modelling Feature (Myers 2012) (Kissinger 2013) (Gassiat 2013) Current Study Spatial Dimension 3D 3D 2D 3D Simulator MODFLOW Dumu x SUTRA-MS ECLIPSE Flow Model Single-phase Multi-phase Multi-species Single-phase Multi-species Multi-phase Single-species Density Changes     Contaminants • Fracturing Fluid • Brine • Fracturing Fluid • Brine • Methane • Fracturing Fluid • Brine • Chemicals in place (e.g. NORMS) • Fracturing Fluid Tracer Used Passive Passive Passive Passive Scenario (s) • Natural upward migration • Flow through fracture • Natural upward migration (fluid and gas) • Flow through fault • Flow through fault • Upward migration through underground failure Simulation time 100 years 12 hours, 30 years, 100 years 1000 years 500 years Conclusion(s) Fluid reaches the aquifer • Gas reaches the aquifer • Fluid reaches the aquifer in extreme case Fluid reaches the aquifer Chemical reaches the shallow aquifer under certain situations Health Risk Assessment     Risk Framework for the Breach Scenario     Uncertainty in Parameters     Page | 27 2.7 Health Risk Assessment from Consumption of Contaminated Subsurface Resources Assessing the risks to groundwater reserves is a challenging task for several reasons. First, there exists a lack of detailed characterization of the geological system (Rubin 2003). Also, there is variability and uncertainty in the hydrogeological, operational, and health parameters that significantly affect the risk assessment (USEPA 2001; Ma 2002) . As a consequence, groundwater-driven health risk analysis should be treated within a probabilistic framework (R M Maxwell & Kastenberg 1999). Probabilistic human health risk assessment due to aquifer contamination has been the topic of intense research studies (Benekos et al. 2007; Ma 2002; López et al. 2008). It is shown that contaminant transport and exposure parameters are the main sources of uncertainty when assessing that risk (Maxwell et al. 1999). The relative contribution of uncertainty reduction in human physiological and hydrogeological parameters on the increased lifetime cancer risk assessment has also been addressed (de Barros & Rubin 2008; de Barros et al. 2009). In a more recent study, the concentration cumulative distribution function (CDF) is analytically derived with the goal of determining the probability of concentration exceedance (de Barros & Fiori 2014). The authors developed the concentration CDF model that incorporates the spatial uncertainty in hydraulic conductivity field, size of the concentration source, and the volume of sampled water at the monitoring well. The concentration probability density functions (PDF) of reactive and non-reactive species in aquifers have also been determined (Sanchez-Vila et al. 2009; Tartakovsky et al. 2009). In a review paper, different methods to probabilistically quantify groundwater risk analysis are discussed in detail (Tartakovsky 2013). Although several studies on groundwater-driven risk analysis are available in the literature, further work is needed to Page | 28 employ and develop these methodologies to the hydraulic fracturing process in order to improve our capacity to predict and control the associated risks. Page | 29 Well Casing Failure and Upward Flow Movement Chapter 3: 3.1 Introduction This chapter introduces the conceptual model, hypothetical scenarios, the simulation results, and risk characterization. Out of the discussed potential environmental impacts, shallow groundwater contamination as a result of failure in the injection system is the subject of emphasis herein. The goal is to evaluate the risks of such incidents to human health and determine how the hydrogeological and well setting control the magnitude of these risks. The conceptual model of this study is defined as a coupled geological system consisting of the groundwater table and the formation underneath. A numerical model is used as a core of the probabilistic risk-based framework to simulate fluid movement and chemical plume spread within the porous media and eventually translate the simulation results into risk measures. The results are reported in terms of three Environmental Performance Metrics (EPMs) at two environmentally sensitive locations – namely the aquifer boundary and the monitoring well. 3.2 Problem Statement A groundwater contamination problem resulted from a failure in the injection system is considered. The contamination is assumed to be released from a point source on the well vertical casing. Potential health problems from the groundwater exposure pathway are assessed through a human health risk framework. In this study, human health risk is termed as the hazard of being exposed to a non- carcinogen for a long time through the drinking water pathway. The adverse health effects are quantified by adapting USEPA’s chronic Hazard Index (HI) for Tetraethylenpentamine (TEPA). TEPA has applications in fracturing fluid as a stabilizer (Gordalla et al. 2013). This compound is Page | 30 known to be non-reactive, non-biodegradable, and stable in the environment (International Programme on Chemical Safety 2001). The main reason for choosing HI over other risk metrics (such as cancer risk) is the limited toxicology studies available for TEPA and many other contaminants of the fracturing fluid. The formulated risk for this study is a function of two different sets of field parameters – namely, operational and hydrogeological. Physical and hydraulic properties of the formation outside the well and near the leakage point define the hydrogeological parameters. The operational parameters, on the other hand, are selected as the leakage rate and the location of a leakage point with respect to the bottom of fresh water reserve. Effects of different parameters on the risk are first investigated through deterministic scenarios. These scenarios will help improve the understanding of the impact of model input parameters on the prediction of interest. Then, uncertainties in the parameters are acknowledged by conducting stochastic analyses through the Monte Carlo simulations. In all deterministic cases, a homogeneous and anisotropic permeability field is assumed. However, effects of spatial heterogeneities in the permeability field are further addressed through stochastic simulations with spatial uncertainty. Uncertainty in the results is quantified by the CDF in order to determine the probability of exceedance of three distinct EPMs that are relevant to risk analysis (de Barros et al. 2012): (1) the chemical concentration at an environmentally sensitive target location (2) the arrival time (from the source to the target location) and (3) the hazard index on the point of exposure. The hazard index is evaluated as a function of the hazard quotient (USEPA 2001). Page | 31 3.2.1 Conceptual Hydrogeological Model The approximate location of the failure point on the well casing, and the geological setting of the hypothetical scenario are shown schematically in Figure 3-1. The conceptual model consists of a series of different geological layers: Vadose zone, aquifer, impervious layer, overburden formation (sand), and the shale layer. The model is conceptualized as a 3D composite geological formation. The shallowest layer, i.e. the vadose zone (Layer I), is placed on top with the thickness of 20 m. A 20 m thick aquifer (Layer II) lies underneath the vadose zone. A thin shale layer (Layer III - 2 m thickness) separates the aquifer from the underneath sand medium (Layer IV) with the thickness of 800 m. This thin layer is denoted here as the impervious layer. The shale reservoir (Layer V) is assumed to be a shallow formation starting from a depth of 840 m and extending to 900 m. It should be noted that the thickness and depth values are based on the properties of shallow shale plays reported by the (USEIA 2011). A monitoring well is placed at a 100 m distance from the hydraulic fracturing well. It is necessary to be able to monitor the water quality at a specific distance from the hydraulic fracturing well. The sand formation over the shale is modeled as an under-saturated porous medium (Saiers & Barth 2012; Cohen et al. 2013; Reagan et al. 2015). As inferred from Figure 3-1, the vertical distance of the failure point to the aquifer bottom is an imperative parameter when studying chemical transport toward the aquifer since the situation becomes more critical if the leakage point is closer to the aquifer bottom. Page | 32 Figure 3-1: Schematic representation of the failure event and the hydrogeological conceptual model: Point source located on the vertical section of the well within the sand formation. The layers include Vadose zone (I), aquifer (II), impervious layer (III), overburden formation (sand layer-IV), and the shale layer (V). One could define the term Critical Depth as the depth below which the contamination plume does not reach the aquifer over a course of time, or it reaches the aquifer in insignificant amounts. The critical distance number is basically the interval during which more caution must be taken to ensure the pipe wall and cementing are, and will remain, in a healthy condition and the likelihood of leakage is reduced. 3.2.2 Formulation of Risk Model For a non-carcinogenic chemical, the hazard quotient (HQ) is defined as the exposure dose to a chemical over a period of time, divided by the daily exposure reference dose. The reference dose Page | 33 is that with a zero likelihood of occurrence for adverse health effects (USEPA 2001). The chronic exposure hazard quotient is given by: RfD CDI HQ  3-1 where CDI is the chronic daily intake [mg/(kg-d)] and RfD is the chronic reference dose [mg/(kg-day)] (i.e. calculated for the lifetime period). The chronic reference dose by definition is the estimation of a daily exposure level for the human population over a lifetime that is likely to be without an appreciable risk of deleterious effects (USEPA 1989). Protection from long-term (7 years to life-time) exposure to chemicals are addressed using chronic reference dose values (USEPA 1989). There are some uncertainty factors to applied in the chronic reference dose. For each noncarcinogenic chemical dose-responses studies will be conducted to determine the Lowest-Observed-Adverse-Effect-Level (LOAEL) and No-Observed-Adverse-Effect-Level (NOAEL). LOAEL is by definition the lowest exposure level where increases in severity or frequency of adverse effects are observed between the exposed population and the control group. These increases must be statistically or biologically significant to be considered. Similarly, the level at which no statistically or biologically significant increases are seen is called NOAEL (USEPA 1989). Uncertainty factors as well as modifying factors are being used to derive RfD from the NOAEL or LOAEL (USEPA 1989). The chronic reference dose can be determined from its Drinking Water Equivalent Level (DWEL) (USEPA 2009). DWEL is calculated by assuming an average body weight of 70 kg and water consumption of 2 L/d (USEPA 2009). For the purpose of this study, drinking water is assumed to be the only source of exposure. The CDI is given by: Page | 34    C CDI 3-2 where C is the chemical concentration (at the environmentally sensitive location) averaged over the exposure duration and  is a coefficient which incorporates the health parameters: AT EF ED BW CR     3-3 where CR is the contact rate of medium (water ingestion rate) [L/d], BW body weight [kg], ED exposure duration [yr.], EF exposure frequency [d/yr.], and AT is the averaging time [d] which is equivalent to 365×ED for non-carcinogens (USEPA 2001). It is assumed that C is calculated at a monitoring well location according to the moving averaged expressed below (R. M. Maxwell & Kastenberg 1999; de Barros & Rubin 2008):              dt t C ED C ED t t 0 0 , 1 max x x 3-4 with 0 t being the starting time of exposure and   t C , x being the chemical concentration at any point in space   x and time   t . It should be noted that defining the concentration as an average value is aligned with the USEPA guidelines, where it introduces the reasonable maximum exposure (RME) as the maximum exposure which is logically expected in a site (USEPA 1989). After calculating the hazard quotient for every single chemical, Hazard Index (HI) can be formulated as follows: Page | 35    n i i HQ HI 1 3-5 With n denoting the number of different chemicals in a compound. For this work, as there is only one chemical, HI is equal to HQ. According to USEPA (2001), HI =1 is referred to as the risk level of concern. It should be noted that risk in this study is regarded as the probability of harmful effects and adverse response induced in human’s body as a result of exposure to an environmental stressor (USEPA 2001). In other words, the risk to human health is defined as the probability of observing HI CDF value of greater than 1 (i.e.   1 Pr  HI ) (Siirila et al. 2012). In general, HI has an uncertainty which is routing from two different sources of risk pathway and health-related parameters such as the exposure duration, time and frequency, the body weight, the population, etc. This study, however, has only been focused on uncertainties from the risk pathway (namely uncertainties in hydrogeological parameters which translate into uncertainty in the concentration term of Equation 3-4). 3.3 Modeling and Numerical Simulation In this section, the ranges of geological and geometrical properties for different formations, governing mathematical equations, and the details of the simulation model are discussed. 3.3.1 Hydrogeological Parameters Shale is a tight formation with permeability values as low as 10 -5 millidarcies (10 -20 m 2 ) (Arthur et al. 2009) and not higher than 10 millidarcies (10 -14 m 2 ) (Freeze & Cherry 1979). The depths at which the formation exists within the United States are various and are dependent upon the geological features of different geographic regions. Barnett, Marcellus, Fayetteville, Haynesville, Woodford, Antrim, New Albany, and Lewis are among well-known U.S. gas shales (Table 2-1). Page | 36 These formations are from 4000 to 7400 feet (1220m to 2255 m) deep on average. New Albany is the shallowest shale formation in this set, with depth values of 500ft to 2000 ft. (150m to 610 m) followed by Antrim (600ft to 2200 ft. – 183m to 671 m) as the second shallowest (Arthur et al. 2009). As suggested by EIA’s study on U.S. gas and oil shales (Table 2-1), shale plays are found in five regions with porosity values ranging from 1% to 12% (USEIA 2011). Oil and gas production from shallow shale formation has long been of interest in the state of California. Formation overlaying the shale is normally more pervious. The geological units forming the overburden differ from one shale play to another but they usually consist of a set of sandstone, inter-bedded shale, siltstone, and mudstone as in the case of Marcellus formation (Saiers & Barth 2012). Generally, the reported permeability values for sandstone formations varies from 10 -17 to 10 -13 m 2 (Gleeson et al. 2011). In another classification by Bear (1988) fresh sandstone shows a permeability variation between 10 -15 and 10 -14 m 2 , whereas oil and gas bearing sandstones have reported values ranging from 10 -13 to as large as 10 -11 m 2 for naturally fractured reservoirs. Porosity of sandstone is in the range of 5% to 30% (Freeze & Cherry 1979). Also, the residual water content of the sandstones can be calculated using empirical relationships among permeability, porosity and residual water of sand formations discussed by Timur (1968). 3.3.2 Flow and Transport Model Formulation The solute transport mechanisms in this research include advection, dispersion, and diffusion. For the case analyzed here, the effects of advection and dispersion are more pronounced during the injection period because of the ongoing fluid movement from the point source into the medium. However, the transport mechanism is mainly diffusive, starting from the end of injection to the time at which the contamination plume diffuses across the impervious layer Page | 37 (Layer III-Figure 3-1) and enters the aquifer (Layer II-Figure 3-1). Chemical movement inside the aquifer is also advective-dispersive. Fundamental equations of flow and transport in porous media for a two-phase flow system are valid in the geological system addressed in this study, as discussed by (Saiers & Barth 2012). Three types of equations are concerned in such a system: conservation of mass (mass balance), conservation of momentum (Darcy’s law), and equation of state (Chen et al. 2006). Here, a 3D variably saturated flow field (within the sand medium) is considered with the domain of interest and the boundary denoted by  and  , respectively. The mass balance equation for the flow in such a system is as follows (Chen et al. 2006):            q t S       u , a w,   3-6 where  is the medium porosity,   is density of each phase per unit volume,  S is the saturation of each phase, t is time,  u is Darcy flux vector in 3 directions, q is sink or source (mass per unit volume per unit time), and D is the diffusion and dispersion tensor (Equation 3-17). The concentration varies in space   z y x , ,  x and time t . Darcy’s law for a two-phase flow system is written follows:   z g p k r             k u 1 , a w,   3-7 withk , , ,   p k r   denoting absolute permeability, relative permeability, pressure, and viscosity for the phase  . The relative permeability ( r k ) for each phase, is a scalar that relates the effective permeability tensor of that phase to the absolute permeability tensor. The physical description of the relative permeability is a tendency of the phase of interest in making the Page | 38 porous medium wet. Brooks-Corey relation is a commonly used empirical relation for acknowledging relative permeability (Brooks & Corey 1964): Wetting phase:     3 2   e rw S k 3-8 Non-wetting phase:               2 2 1 1 e e rnw S S k 3-9 where e S is the effective saturation and  is the pore-size distribution index. The existence of surface tension between two fluids causes the pressure in the wetting fluid to be less than the pressure of the non-wetting fluid. The pressure difference (capillary pressure) is then calculated as (Chen et al. 2006): w a c p p p   3-10 Capillary pressure can also be expressed as a function of w S . One of the most well-known empirical equations for capillary function is the Brooks-Corey equation (Brooks & Corey 1964):    / 1   ew e c S p p 3-11 where c p is the capillary pressure and e p is the air-entry pressure. Equation of state for the system is defined in the form of fluid compressibility relation with volume and pressure changes. For a slightly compressible fluid (Chen et al. 2006):     p p c T T f f e p p V V c               1 1 3-12 where   is the density at reference pressure  p . Rock compressibility can also be considered in the formulations in the same fashion, using the following equation: Page | 39     R R R c dp d p p c dp d c                 1 1 3-13 To solve the system of partial differential equations, boundary conditions must be defined. Boundary condition in this study is a Neumann boundary condition in which the mass flux is prescribed on the boundary (  ): g   υ u  on  3-14 where υ is the unit normal outward to  , υ u   gives the projection of the flux vector on unit normal vector of the boundary, and g is the known mass flux on the boundary. As it will be further explained is section 3.3.3, g will be 0 for all the boundaries of the overburden formation and also for the North and South boundaries of the aquifer layer. West and East sides of the aquifer, however, are ascribed a constant mass flux (refer to Table 3-4). The initial condition is given as:     x x 0 0 , p p  ,   x 3-15 where  denotes the entire domain of concern and   x 0 p is the hydrostatic pressure at point x in space. Finally, solving the system of equations with initial and boundary conditions will give out four unknowns (2 saturations and 2 pressures). The equation used for simulation of chemical transportation the medium is written as follows:     q c c t c         D u     3-16 In which c is the volumetric fraction of chemical in the fluid phase,  is the host phase density. The unit of c is the same as volume since   c has the mass dimension. D is the three dimensional diffusion and dispersion tensor defined as follows: Page | 40         u E u E u I D     t l m d d d  3-17 where m d is the molecular diffusion coefficient, l d and t d are the longitudinal (parallel to the flow) and transverse (perpendicular to the flow in two directions: horizontal and vertical) dispersion coefficients (dispersivities). l d in laboratory is between 0.1 to 10 millimeters and it ranges from 1 to 100 meters or larger in the field. Also, L d is 5 to 20 times larger than t d (Charbeneau 2006). A noteworthy point to be stated here is that in numerical simulations the values of dispersivities in each direction must be smaller than the cell dimension in the same direction. u is the Euclidian norm of u and has the value of 2 2 2 z y x u u u   .   u E is the orthogonal projection along velocity defined as:                2 2 2 2 1 z y z x z z y y x y z x y x x u u u u u u u u u u u u u u u u u E 3-18 and     u E I u E    3-19 In equation 3-16, u  c denotes the advective (convective) process and c  D  represents diffusion and dispersion. Solving the equations 3-7, 3-8 (for water) and 3-16 simultaneously gives the values of pressure and concentration. Numerical simulations of this study are conducted using ECLIPSE. As a well-known oil and gas reservoir simulator ECLIPSE is powerful in solving subsurface multi-phase and multi- component flow and transport equations through finite difference method (Schlumberger 2011). ECLIPSE has also applications in groundwater modelling (Mohammed et al. 2009) as well as Page | 41 environmental assessment and remediation (Zhou & Arthur 1994). This simulator uses Newton’s method, through a fully-implicit approach, to solve the non-linear material balance equations in a grid discretized by finite difference meshing. Using the fully-implicit approach helps in reaching stability when solving the equations over long time frames (Schlumberger 2011). Black oil is the famous model in petroleum engineering that suggests the flow and mass balance equations for a three-phase flow system (oil, water, and gas). Mass transfer between gas and oil phases is allowed in this model, but phase change is not included (Chen et al. 2006). Eclipse 100 is a black oil model simulator used to solve the system of equations for a two-phase flow in the sand medium above the shale formation (oil is not present in the system). A conservative tracer is introduced in the water phase to mimic the chemicals found in the fracturing fluid. The transport equation is then solved and the concentration of the tracer is reported in different time steps. Table 3-1 shows a list of equations solved by ECLIPSE in a discretized grid. Table 3-1: ECLIPSE100 equations used to solve the problem of chemical transport in a two-phase flow system Equation Description Q F dt dM R fl    0         dt t dt t t dt t fl Q F dt M M R Non-linear residual of phase (Mass conservation) for a grid block at each time step * .                               i n gni gi gn i n wni wi wn g g rg w w rw ni ni z z g p p z z g p p B k B k T     0 0 F Flow equation for neighboring cells n and i. B is the formation volume factor (the ratio of the fluid volume measured at pressure, p and temperature, T to the volume measured at standard condition).   w w rw bh iw wi wi i B k p H p T      Q Injection into block i.               dp dB B c 1 ,          dp dV V c R 1 Fluid and rock compressibility   0  X R System of equations to be solved in fully implicit approach (Newton’s method). Page | 42 Equation Description               0 0 g w R R R ,              g g w w B S B S PV M ,        g w S p X Vector of residuals, vector of Mass, and vector of unknowns. PV is the pore volume of the grid block. ij g g w g g w w w j i dS dR dp dR dS dR dp dR d d              X R The Jacobian matrix of the equation (i and j are neighboring cells).   QC l C DA S C gz p Tk B B C VS dt d i w r w w w w                         1 1 Tracer mass conservation in a grid block av D D D D m d m     Hydrodynamic dispersion (dispersion and diffusion combined). a is the dispersivity and v is the velocity. Note that u v   which is the Darcy flux. Parameters used in Table 3-1 are defined as follows: dM : Mass, per unit surface density, accumulated during the current time step, dt F : Net flow rate into neighboring grid blocks Q : Net flow rate into wells during the time step fl R : Non-linear residual for each fluid component at each time step and in each grid block ni F : The flow rate into cell i from neighboring cell n ni T : The transmissibility between cells n and i. Transmissibility between two cell blocks is defined as: l A k T ni   , where  k is the absolute permeability of phase  , A is the face area between blocks n and i, and l is the distance between centers of the two blocks.  r k : The relative permeability of phase   B : The formation volume factor of phase  * By setting 0   R equation will be solved Page | 43   : The viscosity of phase    i n ni i n ni z z g p p dp          : The potential difference for phase  between cells n and i, with p denoting the pressure,  denoting the fluid density, g denoting the acceleration due to gravity, and Z denoting the cell centers depth. iw H : The hydrostatic head correction bh P : The bottom hole pressure  S : Phase  saturation V : Block pore volume C : The flowing tracer concentration  : The porosity D: The hydrodynamic dispersion 3.3.3 Numerical Implementation Three dimensional simulation model of this research includes the aquifer and 60 meters of the sand layer underneath, with horizontal extents of 210 and 205 meters in X and Y directions respectively. It is assumed that the shale is hydraulically fractured somewhere 10 meters below its overburden sand formation, in other words at the depth of 850 m under the ground level. The shale layer, however, is not part of the numerical model. Fluid leakage is considered as a point source contamination on the vertical section of the well. Figure 3-2 illustrates the distribution of the cells in the numerical mesh. Three different sets of numerical block dimensions (cells) are used in this model: Primary cells (5 m × 5 m ×1 m), refined cells (1 m × 1 m × 1 m), and coarse cells (10 m × 10 m × 10 m). The reason for selecting a height value that is 5 folds smaller than lateral lengths is because of the importance of Z direction when reporting the concentration values. Refined or local cell is a smaller sized element used to capture the changes in the area surrounding the leakage point more precisely. The size of this cell is 1 m × 1 m × 1 m. The third cell type is a coarsened one with a size of 10 m × 10 m × 10 m, with application in those areas of Page | 44 the model that are distant to the injection well. Table 3-2 lists three different cell types discussed, along with their geometry and number of frequency in the entire model (Layers and skeleton of the model mesh are further demonstrated in Figure 3-2). Table 3-2: Different grid cells used in the simulation Cell Dimension Number Primary (global) 5 m × 5 m ×1 m 28040 Refined (local) 1 m × 1 m ×1 m 36750 Coarsened 10 m × 10 m ×10 m 1800 Figure 3-2 depicts three areas labelled as 1, 2, and 3. Areas 1 and 3 are those that are not of much interest (lower resolution needed) and they are originally meshed using global blocks and are then coarsened to form bigger blocks. The dimensions of these two cubes are 210 m × 70 m ×60 m and 210 m × 80 m × 60 m, respectively. Area 2, however, is the area of interest and is needed to be fine enough to give high resolution. This area is meshed using a group of global and local blocks to give a second order of refinement around the leakage point. Area 2 is 210 m × 55 m ×70 m in size with the local gird refinement of 35 m × 35 m × 30 m (bottom 10 meters of the aquifer and top 20 meters of the overburden). Figure 3-3 further elaborates on the grid’s specifications. Page | 45 Figure 3-2: Numerical model layers consisting of the groundwater and the geological settings underneath (a) and numerical mesh with injector and the monitoring well (b). An adaptive mesh was employed to accurately capture velocity and concentration gradients. Increasing the grid resolution near the injector is done using ECLIPSE’s Local Grid Refinement (LGR) option. This option enables the user to enhance grid definition near wells using models such as 2D radial, 3D radial or 3D Cartesian. Properties for the refined blocks can either be inherited from the global cells or be defined directly. When using this option, transmissibility from a local grid to the global block in positive X direction is calculated using following equations: J DIPC . TX CDARCY.TML TRANN l lg  3-20 Page | 46 where TRANN lg is the transmissibility from local (l) to global (g) block, CDARCY is the Darcy’s constant equal to 0.008527 in metric units, TMLTX l is the positive X multiplier for global cell g; and DIPC is the dip correction calculated as:         DVS DHS DHS DIPC 3-21 with 2 g l 2 DX DX DHS           3-22 and   2 g l DEPTH DEPTH DVS   3-23 J is defined as follows:           l l g g PERMX DX PERMX DX 2 1 J l A 3-24 where DX is the cell dimension and PERMX is the permeability in X direction. A l is the facing area given by: l l A NTG . DZ . DY l l  3-25 where DY and DZ denote the cell dimensions in Y and Z directions and NTG is the net to gross thickness ratio of the connecting grid block. As previously stated in section 3.3.2, lateral sides of the model are ascribed no-flow boundary conditions, with an exception to the sides of the aquifer perpendicular to the X axis. These boundaries are assigned open-flow conditions to simulate the actual horizontal movement Page | 47 of groundwater flow. Open-boundary is modelled by defining ECLIPSE’s analytical aquifer (constant flux) module on the sides. The type of grid block used for the aquifers is the global. Figure 3-4 shows the layers of the model, along with the place of analytical aquifers. It is important to state that this study deals with one possible boundary conditions setup, though uncertainty in boundary conditions can have remarkable effects on the results of simulation work. Larger grid domains lead to increased computational time and increased effect of the dispersion on the concentration breakthrough curves, so it is important to use appropriately- selected numerical grids for flow and transport simulations. The locally refined mesh of this work is selected after performing a grid refinement sensitivity analysis in which results for numerical grids with block dimensions of 5 m × 5 m, 2.5 m × 2.5 m, and 1.5 m × 1.5 m (all with z dimension of 1m) were evaluated. This analysis revealed that results for locally refined numerical mesh are not significantly different than those of the case with 1.5 m × 1.5 m × 1 m grid blocks. Also, it is verified that the dispersivity values are in accordance with the scale of the grid block to avoid further numerical errors (i.e. dispersivities are smaller than the finest grid block dimension of 1m) (Rubin et al. 1999; de Barros & Rubin 2011). A conservative tracer is used to model the fracturing fluid so that the chemical decay and the physical adsorption of the chemical on the rock surface are neglected. This assumption is in agreement with the study done by Kissinger et al. (2013) where a conservative modelling approach is being used. Considering the fracturing fluid as a persistent chemical goes hand in hand with the extreme hypothetical scenarios of concern in this study. Also, the host phase for the chemical is considered to be water and the viscosity used in the model is taken to be water Page | 48 viscosity, although in reality the fracturing fluid is a power-law fluid and its composition strongly depends on its constituents (Economides & Martin 2007). According to Saiers & Barth (2012), it is a not a valid assumption to consider the overburden formation as a fully saturated medium. In this research the overburden formation is modelled as an under-saturated porous medium in order to respect the multiphase flow system characteristics, and in particular to account for the effects of relative permeability changes due to the presence of air. As previously discussed, two constant flux aquifers are defined on the sides of the model to represent open boundary conditions of the aquifer layer. Simulation time is chosen to be 500 years so that the model can show the mid- and long- run response of the system. Table 3-3 shows the 6 stages and the 5 relaxation times separating each stage. As shown, leakage is presumed to happen during each fracturing stage (Table 3-3). Page | 49 Figure 3-3: Plan view (top), side view 1 (middle), side view 2 (bottom) of the simulation grid. Page | 50 Figure 3-4: Analytical aquifers on sides of the aquifer layer used to imply open-boundary condition across the aquifer in X direction (bottom). Table 3-3: Schedule of the fracturing operation Time Step 1 2 3 4 5 6 7 8 9 10 11 Duration (hr) 2.5 10 2.5 10 2.5 10 2.5 10 2.5 10 2.5 Accumulative Time (hr) 2.5 12.5 15 25 27.5 37.5 40 50 52.5 62.5 65 Accumulative Time (d) 0.10 0.52 0.63 1.04 1.15 1.56 1.67 2.08 2.19 2.60 2.71 Fracturing Stage 1 - 2 - 3 - 4 - 5 - 6 Relaxation - 1 - 2 - 3 - 4 - 5 - 3.4 Simulation Scenarios Contamination transport under the ground in the model is a function of different geological, flow and operational parameters. In order to better understand the system behavior under different situations, a main case scenario has been defined and sensitivity analysis on four different parameters has been conducted. Parameter values used in the main case are shown in Table 3-4. Sand layer is chosen to be a permeable formation which is topped by a thin impervious layer of shale with much lower permeability. The aquifer above is thought to be the most permeable layer of the system, with the permeability value twice that of the sand layer. It is assumed that the aquifer has a constant pressure gradient of 2%, or 2 meters of head in each 100 meters of length in positive X direction. Carrying out the related calculations for the aquifer, it gives out the Page | 51 required Darcy flux of 3.31×10 -3 m 3 /d/m 2 . When divided by the porosity of aquifer, the flow velocity of 1.1×10 -2 m/d or 1.1 cm/d would be achieved. The hydraulic fracturing operation is done in 6 stages, each 2.5 hours in duration, with 10 hours of relaxation between each stage. The injection rate is assumed to be 50 bbl/min (11,500 m 3 /day). Different sets of deterministic scenarios and the parameters range are tabulated in Table 3-5. For the main case scenario, the permeability ratio and initial water saturation of the sand formation are altered. Three values of 0.1, 0.2, and 0.3 have been chosen as water saturation ratios (based on the empirical relations proposed by Timur (1968) and the range of 0.1 to 1 (increase increment of 0.1) has been set for the vertical to horizontal permeability ratio. For all the other deterministic scenarios (sensitivity analysis) discussed in this section, the permeability ratio (Equation 3-26) is fixed at 0.1 and the water saturation is varied systematically. Two geological parameters of the sand layer, namely permeability (k) and porosity (φ), and two operational parameters: leakage point along the injector (H) and the rate of leakage (Q) are selected for sensitivity analysis of the model. The upper value of 20 m for the leakage depth is selected after performing a preliminary set of simulations in which for depths below 20 m very low concentration values for TEPA were observed in the monitoring well. Table 3-4: Parameters used in the main case model Parameter Value Aquifer horizontal permeability 2×10 -13 m 2 (200 mD) Aquifer porosity 0.30 Aquifer Darcy Flux 0.330765 cm 3 /d/cm 2 Aquifer Velocity 1.10255 cm/d Pressure Gradient in Aquifer 0.02 m/m Sand horizontal permeability 10 -13 m 2 (100 mD) Sand porosity 0.15 Tracer molecular diffusion coefficient 10 -5 cm 2 /s Tracer longitudinal dispersivity 0.5 m Tracer transverse dispersivity 0.05 m Tracer initial concentration 440 mg/l Page | 52 Parameter Value Rock compressibility 3×10 -5 1/bar Water compressibility 4.4×10 -5 1/bar Water density 1000 kg/m 3 Water viscosity 10 -3 Pa.s Gas density 1 kg/m 3 Gas viscosity Variable with pressure Brooks-Corey pore-size distribution index (sand) 1.124 Brooks-Corey air entry pressure (sand) 440 pa Number of fracturing stages 6 Duration of each stage 2.5 hr Relaxation time between stages 10 hr Injection rate 11,500 m 3 /d (50 bbl/min) Leakage rate 10 m 3 /d Leakage rate/injection rate 0.09% Leakage point depth to aquifer bottom 10 m Table 3-5: Parameters used for different deterministic scenarios Scenario Sets Parameter Notation Main I II III IV Leakage rate (m 3 /d) Q 5 0.1,3,7 5 5 5 Leakage depth * (m) H 10 10 5, 15, 20 10 10 Sand horizontal permeability (mD) k 100 100 100 50, 150,200 100 Sand porosity φ 0.15 0.15 0.15 0.15 0.05,0.1, 0.20, 0.3 3.5 Result and Discussion In this section the results of different simulations are presented. First, the main scenario and the general contamination plume propagation within the media are discussed. Then, the results of sensitivity analysis and effect of different parameters are shown. For all deterministic scenarios discussed in this section, the permeability ratio (Equation 3-26) is fixed at 0.1 and the water saturation is varied systematically. Finally, with the help of Monte Carlo analysis, the intent is to honor randomness of the parameters. * Below the aquifer Page | 53 3.5.1 Chemical Transport Mechanisms In the current study, transport is a function of different hydrogeological and operational parameters. A main case scenario is defined to assess the effect of initial water saturation (  w S ) and vertical to horizontal permeability ratio (  ) defined as: horizontal vertical k k /   3-26 where vertical k and horizontal k are the vertical and horizontal permeability values, respectively. During fracturing stages, the chemical is leaking into the formation and the transport mechanisms include advection, dispersion, and diffusion. The effects of advection and dispersion are more pronounced when the leak is ongoing, as there is fluid movement in the medium resulting from a high pressure difference between the well and surrounding formation. The effectiveness of the advection is in direct proportionality with presence of water content above the initial water saturation in the medium. In other words, as the medium is unsaturated, relative permeability is a parameter which plays an important role in defining the occurrence and facilitation of water movement. This fact can be deduced from Equation 3-7, where relative permeability and Darcy flux show direct proportionality. When the water content is equal to its initial value, relative permeability of the water phase takes on a zero value -- the initial water saturation is equal to the connate (irreducible or residual) water saturation of the medium. As the leak begins, more water is injected into the medium, water saturation is increased and so is the relative permeability of water. Equation 3-16 clearly shows that the advective term of the chemical mass transport formula is a function of Darcy flux. Therefore, increased Darcy flux values from higher relative permeabilities encourage advection as the governing mechanism. After the injection stops, diffusion becomes the most prominent Page | 54 mechanism for movement and expansion of the contamination plume in the overburden medium. Beside diffusion, a gravitational pressure gradient exists permanently in the medium, which results in downward advective flow in the areas where water saturation is larger than the initial water saturation of the medium – namely, around the leakage point. When chemical plume diffuses into the aquifer, it undergoes a horizontal pressure gradient and advection and dispersion become effective mechanisms once more. Molecular diffusion is also taking place. Different transport mechanisms are better understood by illustrating the simulation results in the following sections. 3.5.2 Defining the Threshold and Health Parameters for the Chemical of Interest In order to define whether or not a situation is critical from the vantage point of health, a chemical threshold (also known as contaminant’s reference level) must be defined and the simulated concentrations must be evaluated against that definition. The chemical of concern in this study is Tetraethylenepentamine (TEPA), which is used in fracturing fluid as a stabilizer even though it is known to be a toxic compound. This agent is stable and remains in the environment for a considerable length of time. There is no complete human toxicology database available, but according to the German Federal Environment Agency, for chemicals with some toxicological information, the reference dose of concentration is the same as the health related indication value (Gordalla et al. 2013). More information on this chemical and the calculations of threshold are discussed in Appendix A. The chemical threshold value is used to compare the results of the sensitivity analysis as well as the stochastic scenarios. The sensitivity analysis basically entails comparing the risk measures at the monitoring well and at the end of year 500. Page | 55 Health related parameters used in this study (see Equations 3-1 to 3-5) are shown in Table 3-6. Table 3-6: Health-related parameters used in the model Parameter Value Contact/ingestion Rate 2 L/d Exposure duration 30 y Exposure frequency 365 d/yr. Body weight 70 kg Averaging time 10,950 d Reference dose 8.57×10 -6 mg/kg.d 3.5.3 Main Case Result Selected two-dimensional snapshots of the main case model with initial water saturation (  w S ) of 0.3 and permeability ratio (  ) of 0.1 are presented in Figure 3-5. Pictures are showing the central X-Z cross-section of the model. Eight time-steps are shown, starting with the one at the end of the initial fracturing stage (hour 2.5) followed by the end of stage 6 (hour 65 or day 2.71), day 1000 (year 2.74), day 10,000 (year 27.4), day 36500 (year 100), day 109,500 (year 300), day 146,000 (year 400), day 182,500 (year 500). As is shown in the first two figures (top) the plume shape is almost the same in the first two time-steps. But the size of the plume in the right picture (hour 65) is larger. At year 2.74, the plume is expanded more within the medium from the apparent effect of gravity. Increase of water saturation locally around the leakage point makes the water relative permeability active and in turn results in the plume moving downward due to the gravity pressure gradient. In the next time-step (year 27.4) the plume front has passed the impervious layer and has found its way through the aquifer, where advection and dispersion start working, and give a perceptible asymmetry to the shape of the plume, as well as an apparent extension toward the flow direction in the aquifer (positive X). Page | 56 After 100 years, the aquifer shows chemical concentration ranging from 10 -5 to 10 4 mg/l. Some amounts of chemical diffuse back from the aquifer into the impervious layer, and some are transported of the model from the open boundary. The concentration plume located in the sand medium is still feeding the aquifer and as seen in the next snapshot, by the end of year 300, concentration gradient in the aquifer increases as a direct outcome. Although when compared to the previous time-step, it is clear that concentration plume in the sand is expanding even more with the warm colors turning into colder ones, especially at the plume center (i.e. where the leakage has originally taken place). Diffusion back into the sand from the aquifer bottom is still ongoing, which helps to mitigate the unfavorable condition by removing concentration from the aquifer. In the last snap-shots, concentration gradient of 10 -5 to 1 mg/l has filled considerable space of the aquifer. However, the concentration plume strength is diminishing and less chemical is leaking into the aquifer (the concentration value reaching the aquifer boundary is decreasing in this time-step as it is more elaborately shown in Figure 3-7). The trend of concentration value reaching the aquifer during the simulation time is increasing at first, up to a peak value, after which the concentration shows a decreasing trend as the plume in the sandy formation is weakened, and the concentration directly below the aquifer is lessening. In addition, some of the mass in the aquifer is taken out of the medium from the open boundary and some of it diffuses through the aquifer bottom back down. This phenomenon is further explained in the next section. Chemical movement within the aquifer is investigated more accurately in Figure 3-6, where a snapshot of the model at the end of year 300 is shown. Red arrows represent concentration movement directions from the plume front reaching the impervious layer beneath the aquifer, and the concentration movement within the aquifer. In addition to the arrows shown, Page | 57 there are concentration movements along Y axis (positive and negative directions) which are not depicted here, although they do exist and play an important role in reducing the concentration value. To elaborate more on the transport mechanisms and system behavior, Figure 3-7 and Figure 3-8 compare the concentration-time results for different permeability ratios and initial water saturations in the main case scenario. Concentration versus time figures are plotted for two different locations in the model: bottom of the aquifer (termed as aquifer boundary), and the monitoring well location. Additionally, concentration values are normalized with respect to the initial concentration of 22 kg/m3 (22,000 mg/l). In another graph representation, contour plots of the simulations are shown in Figure 3-8, where it is clear that the most critical situation (higher concentration in the aquifer) happens in the lowest permeability ratio and the highest initial water saturations of the sand. Moving towards fully homogenous medium ( 1   ) lowers concentration values as discussed later in this document. Page | 58 Figure 3-5: Cross-section of the main case model in various time-steps. The results are visualized for the case with 3 . 0   w S and 1 . 0   Concentration (mg/l) 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 2.5 hr 65 hr 2.74 yr. 27.4 yr. 100 yr. 300 yr. 400 yr. 500 yr. X Z Page | 59 Figure 3-6: Concentration plume expansion after 300 years. Movement directions of the chemical mass are shown by red arrows. Figure 3-7: Temporal evolution of the concentration at the locations of interest: aquifer boundary (a) and monitoring well (b). Results presented for different initial water saturations (  w S ) and different anisotropy ratios (  ). The orange vertical line depicts the end of injection process. Concentration (mg/l) 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 (a) (b) X Z Page | 60 Figure 3-8: Contour plots of normalized concentration at the sensitive locations as a function of time and  . The sensitive locations considered are: aquifer boundary (a) and monitoring well (b). Concentrations shown for different initial water saturations (  w S ) and anisotropy ratios (  ) over the simulation time. The concentration values are normalized by 0 C . The following sub-sections investigate the mechanisms of contaminant transport over the course of simulation, as well as the effect of changes in permeability ratio and initial water saturation of the overburden medium. (b) (a) Page | 61 3.5.3.1 Concentration-Time Curves: General Trend Figure 3-9 shows the general concentration-time plot for the aquifer boundary of the case with 1 . 0   and 3 . 0   w S . Concentration reaching the aquifer bottom remains zero for almost 40 years; later it increases up to a peak value. Thereafter, a decreasing trend is seen towards the end of the simulation. Figure 3-9: Concentration-time curve general trend for the case with 3 . 0   w S and 1 . 0   As discussed previously, and as seen in Figure 3-5, it takes time for the plume front to reach the lower boundary of the aquifer. An increase in concentration value in Figure 3-9 is a direct result of the diffusion from the plume front. As time evolves, the front covers a larger area of the impervious layer beneath the aquifer with elevated concentration values, so that the diffusion into the aquifer becomes more drastic. As inferred from the slope changes, the concentration increase rate (i.e. tangent to the curve) is low at first, then grows and roughly remains constant from year 90 up until year 240, then starts decreasing until year 320, where it approaches zero near the peak value, and finally goes negative from the peak point on. Declines Page | 62 in concentration values are mainly attributed to the mass displacement and advection through the boundary of the aquifer, and because of the diffusion of lower concentration values into the formation beneath the aquifer. The main transport mechanism of the concentration movement is diffusion, with an existing tendency for the plume to move laterally to the distant areas from the start of the simulation. However, as time elapses, more concentration is displaced and the value is decreased within the plume center, so that the plume is not as strong as it used to be and is not able to provide enough concentration into the aquifer to maintain the concentration value. As seen in the snap-shots from simulated model, the plume tends to fill up the empty rooms throughout the sand medium (Figure 3-6). The concentration-time plot of the monitoring well shows behavior similar to the aquifer boundary, but with a shift in time. Another difference is that unlike the case of the aquifer boundary, where the peak shows sometime in the middle of the simulation, the concentration in the monitoring well is close to its highest level by the end of the simulation. 3.5.3.2 Concentration-Time Curves: Effect of Permeability Ratio Each graph drawn for a specific permeability ratio shows the same general behavior. When comparing graphs for two different permeability ratios with fixed water saturation, however, one can see that at the start concentration values are directly proportional to the permeability ratio value. But as time elapses, the curve with lower  takes higher concentration values. Figure 3-10 demonstrates the concentration-time curves for initial water saturation of 0.3 and permeability ratios of 0.1 and 1.0. The spatial location of the reported concentrations is the aquifer boundary. Page | 63 Figure 3-10: Effect of permeability ratio on concentration reaching the aquifer boundary – comparing simulations with constant  w S of 0.3 and varied  of 0.1 and 1.0 In the beginning, both graphs show increases in concentration, but the values associated with the permeability ratio of 1 are larger due to the effect of permeability on facilitating water movement. In longer time-steps, however, diffusion becomes the governing mechanism and the concentration values for lower permeability ratio (0.1) become larger. This is primarily due to the shape of the concentration plume in two different cases. Figure 3-11 shows the plume geometry for two permeability ratios at different time-steps. There are main differences between the plume geometry in the two cases. When injection stops at hour 65, the plume is more spherical in 1   whereas the plume in 1 . 0   is ellipsoidal, clearly as a direct result of the permeability values. Hence, in time-steps below the crossing point, concentration reaching the aquifer grows faster and arrival-time is shorter (to the left of the crossing point in Figure 3-10). Furthermore, effect of gravity is more pronounced in 1   -- as it seen in time-step 2.74 years and beyond. Also, chemical reaches the open boundary of the aquifer sooner in the case of 1   and mass displacement across the open boundary starts sooner in this case. That is one of the reasons Page | 64 why the curve associated with 1 . 0   is steeper than the one of 1   after the crossing point, which further means the rate of increase in concentration is higher. At the crossing point the reported concentrations are very similar in both models, but the higher rate of increase observed in 1 . 0   is predominantly due to the stronger concentration plume feeding the aquifer from the layer beneath. As mentioned earlier, in the fully homogenous case a large part of the plume is pulled downward by gravity. As a result, the plume in the model with 1 . 0   does not expand as much as the one where 1   . The core remains near the aquifer bottom and keeps feeding the aquifer, whereas in the bottom model the plume is stretched downward. The effectiveness of the plume is, thus, weakened naturally in the fully homogenous case. In addition, mass transport out of the aquifer is more significant in the lower model, since the concentration reaches the open boundary sooner. Also, lower growth rate of concentration makes the curve for 1   to peak at a value less than the peak of 1 . 0   . Another noticeable event is the effect of concentration plume spread in the top of impervious layer in the case of 1 . 0   which results in a more drastic situation at the bottom of the aquifer. Figure 3-12 shows plan view of the top of the impervious layer and right beneath the aquifer for three time-steps of 2.74, 27.4, and 300 years. As shown in the top left figure, before the crossing point (2.74 years) plume is slightly less strong in the case of 1 . 0   . Later at the crossing point the same behavior still governs, indicating that the case of 1 . 0   is associated with lower concentrations in the impervious layer. After the crossing point the concentration is showing higher values in the case of 1 . 0   making the aquifer boundary concentrations larger. At 300 years, concentrations observed in the top of impervious layer originated from two different sources: 1- initial leakage point and upward migration of the chemical from the bottom, and 2- from the chemical which is already in the aquifer and diffuses back into the overburden formation. Finally, in the last two time-steps of Figure 3-11 one sees Page | 65 that the plume has reached its maximum expansion; meanwhile, diffusion, dispersion and mass displacement are all functioning as transport mechanisms for the chemical. Page | 66 1 . 0   1   Figure 3-11: Effect of permeability ratio on concentration plume geometry. Illustrations shown for the case with constant initial water saturation ( 3 . 0   w S ). Concentration (mg/l) 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 65 hr 2.74 yr. 27.4 yr. 100 yr. 300 yr. X Z 500 yr. Ellipsoidal shape Spherical shape Page | 67 1 . 0   1   Figure 3-12: Effect of permeability ratio on concentration plume geometry. Plan view of the aquifer bottom is shown for the case with constant initial water saturation ( 3 . 0   w S ). 3.5.3.3 Concentration-Time Curves: Effect of Initial Water Saturation Concentration versus time curves for two different cases are compared in Figure 3-13. Both cases share the same permeability ratio of 0.1, but in order to investigate the effects of water content the initial water saturations are different (0.1 and 0.3). For most of the simulation period, concentration values of 3 . 0   w S show higher levels when compared to the ones of 1 . 0   w S . In general, the curves are not substantially different in shape and the trend is the same as the Concentration (mg/l) 0 4.87×10 -9 9.75×10 -9 1.46×10 -8 1.95×10 -8 Concentration (mg/l) 0 2.95×10 -4 5.9×10 -4 8.85×10 -4 1.18×10 -3 Concentration (mg/l) 0 0.16 0.32 0.49 0.65 2.74 yr. 27.4 yr. 300 yr. X Y 1 2 1 2 Page | 68 customary trend discussed earlier: concentration starts going up, takes a peak value and then decreases. The reason pertains to the water permeability in a two-phase flow system which is defined as the product of relative permeability of the water phase and the absolute permeability of the formation ( abs w r w k k k   ). As deduced from the equation, when the absolute permeability is fixed changes in relative permeability must be very significant to result in substantial changes in water permeability. However, when calculating the relative permeability values from the Brooks-Corey equation, one notices that the case with 3 . 0   w S is associated with slightly larger numbers of r k when compared to the values of 1 . 0   w S . Almost the same numbers of relative permeability and the fixed value of absolute permeability explain why the curves behave similarly. In another event, as discussed in the last section, changing absolute permeability results in more clear differences in the two curves drawn for the same water saturation (Figure 3-10). Observing the change in the contaminant plume shape over the first 11 time-steps (toward the end of the fracturing job), clearly indicates that the plumes in both cases are almost identical. However, the core of the plume in the case with more water content shows lower concentrations as opposed to the case with less water saturation. That is because even though the same mass of chemical leaks into the formation in both cases, the presence of more water in the case with higher water saturation results in more diluted concentrations. Throughout the first 65 hours (end of the injection) diffusion is not yet too effective; however as time-steps increase, the plume spread driven by diffusion can be observed. On the other hand, the slight difference between the two curves is a direct result of the change in the water content ratios, which makes the diffusion across and within the impervious layer more effective (the porosity of the impervious layer is 0.1 and that of the sand is 0.15). Thus, even though increased water content Page | 69 can mitigate a harmful situation by diluting the chemical, it also provides larger and more effective pathways for the chemical, which in turn facilitates its transport. Figure 3-13: Effect of initial water saturation on concentration reaching the aquifer boundary – comparing simulations with constant  of 0.1 and varied  w S of 0.1 and 0.3 Figure 3-14 shows snapshots of the model for two cases of water saturation. In the end of first two time-steps (65 hours and 2.7 years) the geometry of the plume remains unchanged by shifting the water content from 0.1 to 0.3; however the concentrations are slightly lowered in the latter case. Next time-step (end of year 22) reveals the main difference between two cases, where by focusing at the top of the plume it is noted that the diffusion phenomena in two cases are not happening with the same strength. At that moment, the plume front has reached the impervious layer at the top of the overburden and the diffusion rate is decreased (Figure 3-15). For the rest of the simulation as seen in Figure 3-14, the plume continues moving upward and hits the aquifer. Then, it stretches to the right in accordance with the dominant flow movement direction in the aquifer. The difference in concentrations within the aquifer in both cases is indistinguishable Page | 70 from the snap-shots; however, the lesser values are still visible in the overburden formation in the case of 3 . 0   w S . Figure 3-14: Effect of initial water saturation on concentration plume geometry ( 1 . 0   ) Page | 71 1 . 0   w S 3 . 0   w S Figure 3-15: Effect of initial water saturation on concentration plume geometry ( 1 . 0   ) – impervious layer influence on diffusion phenomenon 3.5.3.4 Risk Quantification of the Main Case (Hazard Index Values) The simulated concentrations can now be translated into the risk measure (HI). The changes of HI with the permeability are illustrated in Figure 3-16. For all of the  values, HI shows numbers larger than 1, indicating a potential hazard. Similar to the case of concentration, increasing water saturation and the anisotropy ratio result in higher values of HI. Figure 3-16: Impact of anisotropy ratio (  ) and initial water saturation (  w S ) on the Hazard Index (HI). Values of HI remain above 1 for all the scenarios indicating a potential hazard in the consumption of drinking water. Concentration (mg/l) 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 22 yr. X Y Aquifer Sand Impervious Layer Page | 72 It is worthwhile to express that a metric like hazard index makes it possible to see the risk in terms of a binary indicator and determine whether or not a situation is deemed to be dangerous. In the case of Figure 3-16, irrespective of what the travel time or concentration values are, the hazard index is high (i.e. greater than 1) indicating a potential risk for the people exposed. 3.5.4 Effect of Changes in Leakage Rate (Scenario Set I) As seen in Figure 3-17, the scenario with a leakage rate of 0.1 m 3 /d concentration values at the aquifer boundary for all the water saturations pass the threshold. Concentrations of the monitoring well are slightly above the threshold. As there is a distance between the first point at which chemical enters the aquifer (aquifer boundary) and the point where the monitoring well is located, it takes longer for the concentration to reach the monitoring well – thus, the concentration value of the well is usually lower than that of the aquifer boundary. When comparing the results for different saturations, the green curve ( 3 . 0   w S ) shows slightly higher concentrations. The reason is the slightly larger relative permeability numbers of this case when compared to 1 . 0   w S and 2 . 0   w S . By increasing the leakage rate both concentrations show an increase in value as normally expected. At small leakage rates, the jump in the concentration values is high from one leakage rate to a larger one, whereas concentration values are less sensitive to changes in large leakage rates. Page | 73 Figure 3-17: Concentration values at the end of simulation (year 500) – Effect of the leakage rate (Q) 3.5.5 Effect of Changes in Leakage Depth (Scenario Set II) The deeper the leakage points the lower the concentration values reaching the aquifer. From depths between 5 and 15 meters below the aquifer a sharp decreasing trend is observed (Figure 3-18), but the rate significantly declines for depths below 18 meters so that the concentration values reaching the monitoring well from a leakage point of 20 meters below the aquifer are quite infinitesimal. When comparing the concentration values to the threshold line, it is clear that even for leakage depths around 20 meters below the aquifer, the concentration reaching the aquifer boundary is slightly higher than the specified threshold; however, in the case of the monitoring well, the concentration value is on the safe side. Main Case Page | 74 Figure 3-18: Concentration values at the end of simulation (year 500) – Effect of the leakage depth (H) As far as the water saturations are concerned, from the leakage depth of 10 m and beyond, a direct relation between the concentration and the water content value is seen, with 3 . 0   w S showing the largest values both for the aquifer boundary and the monitoring well. Whereas, at depths less than 10 m, a change in behavior takes place; at a depth of 5 meters below, 2 . 0   w S shows the highest concentrations followed by 3 . 0   w S and 1 . 0   w S . The underlying reason for this discrepancy is the ease and higher speed of concentration movement in 3 . 0   w S in which the concentration growth rate is the largest, as is the concentration decline rate (positive and negative slopes of the top graphic of Figure 3-19). Figure 3-19 also shows how the smaller concentration decline rate of case 2 . 0   w S ends up showing higher concentration values remaining at the conclusion of simulation. In another type of behavior, and as depicted in the bottom graphic of Figure 3-19, increasing the depth of leakage to 10 m dramatically drops the concentration values for all the saturations, while showing the same behavior for the plots. Main Case Page | 75 In one last representation, Figure 3-20 illustrates the semi-log time-concentration plots of all the leakage depth scenarios compared to the threshold concentration for the monitoring well. The water saturation is kept constant ( 3 . 0   w S ) for this set of plots in order to simply monitor the effect of leakage depths on the results. Of the simulated cases, a leakage depth of 20 m shows concentrations below the threshold value for the entire simulation time. However, the concentration still exists within the aquifer body and it crosses the limit at areas closer to the entrance point of the plume located at the aquifer bottom. On the other hand, the concentration in cases with leakage depths of 5, 10, and 15 m crosses the threshold line in 15, 120, and 325 years after the start of simulation. As expected, the smaller the leakage depth, the sooner the concentration exceeds the threshold. Page | 76 Figure 3-19: General trend of time-concentration plots (breakthrough curves) for leakage depths 5 (a) and 10 meters (b) (a) (b) Page | 77 Figure 3-20: Breakthrough curves for various leakage depth scenarios 3.5.6 Effect of Changes in Absolute Permeability (Scenario Set III) Increasing the horizontal absolute permeability has a positive effect on the result, as higher values for the permeability are linked to smaller concentration numbers reaching the aquifer and the monitoring well. In scenarios with higher horizontal permeability, the vertical permeability is also larger, but effects of the increase in horizontal permeability is greater as the permeability ratio is fixed ( 1 . 0   ) and the horizontal permeability is ten-fold that of the vertical. Thus, the order of growth for the horizontal permeability is 10 times as much, which makes it more effective. Larger horizontal permeability values guide the flow movement to the sides of the model more significantly and make the situation less critical as the chemical gets entrapped in the formation, and the upward migration is limited. However, when comparing the results to the threshold value it is clear that all the scenarios are beyond the acceptable range (Figure 3-21). One last point to be discussed here is that change in the concentration is nearly smooth and constant when the horizontal permeability is increased. Page | 78 Figure 3-21: Concentration values at the end of simulation (year 500) – Effect of the absolute horizontal permeability of the sand layer (k) 3.5.7 Effect of Changes in Porosity (Scenario Set IV) Keeping the total volume constant and increasing porosity result in increased pore volume. Larger pore volume of the sand medium provides more spaces for the concentration to spread. Also, for a fixed initial water saturation, when pore volume is larger, water volume in the medium is bigger and the concentration is lowered. According to Figure 3-22, both the aquifer boundary and the monitoring well show a slight decreasing trend by increasing porosity. Nevertheless, all the scenarios are beyond the acceptable level for the chemical of concern. Main Case Page | 79 Figure 3-22: Concentration values at the end of simulation (year 500) – Effect of the sand layer porosity (φ) 3.5.8 Deterministic Risk Quantification for the Most Critical Scenarios Of the four parameters discussed in sections ‎ 3.5.4 through ‎ 3.5.7, it was shown that the leakage depth has the strongest influence on the concentration prediction (see section ‎ 3.5.5). Here, the HI is calculated deterministically for different leakage depth scenarios. Figure ‎ 3-23 shows the impact of the leakage depth on the HI. Large numbers of HI are observed for a leakage depth of 5 m, whereas by moving to a depth of 10 m, values drop by nearly fivefold. HI equals one at the leakage depth around 18 m, so the depths below 18 m can be regarded as safe depths from the risk assessment standpoint; accordingly, the critical depth is set as 18 m below the aquifer. As demonstrated in Figure ‎ 3-23 for a leakage depth of 5 m the case of 3 . 0   w S shows the largest value, indicating that it is still the most critical water Main Case Page | 80 saturation, irrespective of the concentration value at the end of the simulation being higher for 2 . 0   w S . Figure 3-23: Changes of hazard index (HI) with leakage depth (H), shown for various water saturations. The solid gray line shows HI value of 1 3.5.9 Parametric Uncertainty Analysis and the Stochastic Characterization of Risk Uncertainty analysis is imperative for this study because of the highly uncertain parameters of interest used in the model. The goal of uncertainty analysis is to address most of the possible scenarios by conducting numerous simulations, each comprising a set of parameters. Figure 3-24 schematically illustrates the sources of uncertainty in the model. Uncertainty can be traced from two categories of parameters: hydrogeological and operational. Page | 81 Figure 3-24: Schematic representation of sources of uncertainty in the simulation model One of the most famous methods of incorporating uncertainty in groundwater contamination is Monte Carlo simulation (Meyer 1994; Ballio & Guadagnini 2004). In the Monte Carlo method, typically a significant number of realizations are simulated and the results are reported as probabilistic functions. Here, the first step toward a Monte Carlo run is to define the parameters of concern probabilistically; in other words, defining the parameters as random variables with specified probability density functions (PDF). Samples are then drawn from each PDF, a physical model is built based upon the selected parameters, and then the scenario is simulated by using ECLIPSE. Figure 3-25 shows the framework of the random analysis. Page | 82 Figure 3-25: Workflow of the random simulations - Coupling MATLAB and ECLIPSE Four parameters of interest are chosen for the stochastic analysis: leakage rate, leakage depth, sand layer absolute permeability, and porosity. From the results of the previous sections, it can be concluded that in cases with the lowest permeability ratio and highest saturation, there are higher concentration levels reaching the aquifer and, subsequently, the monitoring well. Accordingly, the focus of the stochastic simulation is selected to be on the cases with 1 . 0   and 3 . 0   w S . In order to choose appropriate PDFs for porosity and permeability of the sandstone, one should notice how these parameters are distributed in the field. Data from well-log samples show that for both sandstone and carbonate formations the permeability distribution follows a log-normal pattern, whereas the porosity tends to be normally distributed (Nelson 1994; Hohn 1999). Therefore, the permeability in the model is presumed to be log-normally distributed. Porosity is, however, sampled from a truncated normal distribution to ensure that negative values are assigned zero probability of selection. Leakage depth and rate are assumed to follow a uniform distribution (Rish 2005). To capture the effect of the aforementioned parameters and to make comparisons among different scenarios, hydrogeological features of the aquifer are set to be constant at the present time. The focus is on the sand layer surrounding the leakage point. The parameters are altered within a Monte Carlo based framework and CDF for three commonly used EPMs (de Barros et Samples from each distribution ECLIPSE run for each set of drawn samples Saving the results Page | 83 al. 2012) are evaluated at the environmentally sensitive location (i.e. the monitoring well). EPMs of interest in this study are the chemical concentration, the source-to-target arrival time, and the HI (see Equation 3-5). The analysis is performed with constant permeability ratio ( 1 . 0   ) and water saturation ( 3 . 0   w S ) assumptions. Statistical distributions of the uncertain parameters are shown in Table 3-7. The Monte Carlo simulation is carried out for 1000 realizations. Table 3-7: Statistical distributions for uncertain parameters Uncertain Parameters Leakage rate (m 3 /d) Uniform U(0.1,7) Leakage depth (m) Uniform U(10,30) Sand permeability (mD) Log-normal lnN(4.61,0.2) Sand porosity Truncated Normal N(0.125,0.025) φ (0.05,0.20) As part of the simulation results, Figure 3-27 shows the histograms of concentrations reaching the monitoring well at 8 different time-steps. Until year 27.4 all the realizations show zero concentration. In year 100, concentration values start acquiring non-zero numbers. As time elapses it can be observed that the concentrations start growing and showing higher values. CDF of the concentrations are also calculated to quantify the probability of crossing the threshold in different time-steps. Figure 3-28 shows the CDF, denoted by F, for eight distinct times at the monitoring well location. The threshold concentration for TEPA is depicted by the dashed red line and is chosen to be 3×10 -4 mg/l (Gordalla et al. 2013). The probability of exceedance (e.g. events in which the concentration is above the threshold) is defined as follows:       0 0 0 0 0 / 1 / / Pr 1 / / Pr C C F C C C C C C C C th th th         3-27 with th C denoting the threshold concentration, 0 C the initial concentration injected, and F probability value from the CDF. Page | 84 At the first four time-steps shown, the probability of passing the threshold value is zero. For the fifth time-step (year 200) the probability jumps to 0.15 and  further increases to 0.27, 0.38, and 0.46 for years 300, 400 and 500, respectively. The exceedance probability as a function of distinct times is plotted in Figure 3-26 . Figure 3-26: Temporal evolution of the probability of exceedance  (Equation 3-27).  measures the probability of observing concentrations equal or higher than the concentration threshold th C at different times. In summary, it is concluded that the Monte Carlo analyses for 1000 scenarios, with previously specified random parameters show a high probability of crossing the threshold value for large time-steps. Page | 85 Figure 3-27: Histogram of concentration reported for the monitoring well at time-steps of interest Page | 86 Figure 3-28: CDF of concentration reported for the monitoring well at time-steps of interest Page | 87 Next, the CDF of source-to-receptor arrival time is investigated (Figure 3-29). In this research, the arrival times are defined as the earliest period of time in which the simulated concentration in the monitoring well takes on a value equal to or greater than the chemical threshold. For the scenarios investigated, on average it takes 276 years for the contamination plume to reach the monitoring well, which is quite noticeable. Figure 3-29: CDF of the arrival time for the threshold concentration. Results illustrated for fixed water saturation and anisotropy ration (  w S =0.3,  =0.1). The CDF starts from 100 yrs as the first time in which concentration threshold ( th C ) is observed in the monitoring well. HI CDF (Equation 3-5) is plotted in for the last EPM of this study. HI values in Figure ‎ 3-30 are highlighted for the percentiles of interest in risk assessment (USEPA 2001). The probability of remaining in the safe zone, i.e.   1 Pr  HI , is 0.55 and the risk is 0.45, accordingly. By setting the critical HI at 1, one could use this CDF plot as a tool to answer questions such as: What are the safe margins and intervals for uncertain parameters? How do operational parameters in hydraulic fracturing affect the risk? Page | 88 Figure 3-30: Hazard index (HI) CDF plotted for fixed water saturation (  w S =0.3) and fixed anisotropy ratio (  =0.1). The risk level of concern (dashed red line) and percentiles of interest in risk assessment (solid and dashed gray line) are highlighted. 3.5.10 Effect of Spatial Uncertainty (Heterogeneity) of the Media on the Characterization of Risk This section will regard the impact of spatial uncertainty in the sand layer and aquifer on the EPMs. Here, one step further is taken towards a complete geological modeling of the system by honoring the heterogeneous characteristic of the porous media. In section 3.5.10.1, the permeability of sand formation is modelled using realizations constructed through geostatistical methods and the Monte Carlo analysis is performed similar to section 3.5.9 with random porosity, leakage rate, and leakage depth with PDFs defined in Table 3-7. Then the scenario with shortest arrival-time is selected from section 3.5.10.1 and associated parameters (i.e. sand layer heterogeneous permeability field, leakage rate, leakage depth, and porosity) are implemented in the model. Spatial uncertainty of the aquifer permeability field is, next, honored through 1000 realizations (section 3.5.10.2). For both sand and aquifer formation the permeability is deemed to follow a log-normal distribution (Nelson 1994; Hohn 1999). The mean permeability value for the Page | 89 sand layer and aquifer are assumed to be 150 mD and 200 mD, respectively. Sequential Gaussian Simulation (SGSIM) is used to randomly create 1000 realizations for each case. The sequential simulation formalism is used in SGSIM algorithm to generate a Gaussian random function. For a multivariate Gaussian random function (   u Y with u denoting the coordinates vector) and a given variogram model realizations will be generated using steps presented in Table 3-8. Table 3-8: Steps of the sequential Gaussian simulation algorithm SGSIM Algorithm I: Defining a random path which visits all the nodes of the grid II: for each node u on the random path do III: Get the conditioning data from previously simulated nodes IV: Local conditional CDF (CCDF) estimation as a Gaussian distribution (mean and variance: from krigging mean and variance) V: New value drawn from the CCDF and added to the data set VI: end for V: Repeat for another realization For the geostatistical analysis an omnidirectional exponential variogram is selected with constant range of 30 m for both cases. For the case of aquifer heterogeneity, the variogram sill (aka correlation variance) is altered to investigate the effects on the risk measure. Nugget effect is assumed to be 0. The mathematical formula for the variogram used is as follows:              a h h 3 exp 1  3-28 with h being the coordinates offset vector or lag vector and a representing the range parameter. 3.5.10.1 Heterogeneous Sand Layer, Random Parameters, and Homogeneous Anisotropic Aquifer Figure 3-31 shows four sample realizations of the sand layer permeability field. The stochastic analysis is performed for the variogram sill value of 0.4. The intention is to observe the effect of increasing heterogeneity on the risk measure. Figure 3-32 shows the concentration CDFs at Page | 90 various time-steps for two extreme values of variogram sills (i.e. 0.4 and 4). It is noted that changing the variogram sill of the permeability sill does not show a significant influence on the results. As an example, the probability of exceeding the threshold for cases of 0.4 and 4.0 are 0.06 and 0.08, respectively. Figure 3-31: Sample realizations for the sand layer permeability field. Page | 91 Figure 3-32: Comparing concentration CDFs for two extreme values of permeability field variogram sill in the sand layer: 0.4 (Left) and 4.0 (right). Values are slightly changing by as the sill is altered. In the following, the CDF of chemical arrival time from the source to the receptor for the two extreme cases of variogram sill is presented (Figure 3-33). On average it takes 291.7 years for the case with sill value of 0.4 to observe concentrations higher than accepted level, whereas this value slightly decreases to 287.1 years in the case of sill=4.0. Figure 3-33: Comparing source-to-receptor arrival time for two extreme values of variogram sill in the sand layer: 0.4 (Left) and 4.0 (right). Finally, comparing the results for the risk measure (HI) reveals that the probability of observing values higher than the critical HI of 1 is 0.23 in the scenarios with sill value of 0.4 and 0.24 in the case with sill value of 4.0. Page | 92 Figure 3-34: Comparing Hazard Index for two extreme values of variogram sill in the sand layer: 0.4 (Left) and 4.0 (right). HI value is almost deterministic for the lower value of sill. 3.5.10.2 Heterogeneous Sand Layer and aquifer In this section, it is strived to investigate the effect of uncertainty in the aquifer permeability field on the EPMs. Out of the simulations performed in section 3.5.10.1, the one with lowest source to receptor arrival-time (i.e. 170 years) is selected. This scenario is associated with parameters as follows: leakage rate of 2.84 m 3 /s, leakage depth of 23 m, and porosity of 0.15, and a specific heterogeneous permeability field. The focus is on acknowledging spatial uncertainty within the aquifer through 1000 realizations of permeability fields each with a log-normal distribution and a mean of 200 mD (Figure 3-35 shows four sample snapshots of the permeability fields). Page | 93 Figure 3-35: Sample realizations for the aquifer permeability field. The analysis in this section is carried out for four different correlation variances of 0.4, 1, 2 and 4 to capture the effect of heterogeneity in the aquifer. Figure 3-36 illustrates normalized concentration CDFs for two extreme values of correlation variances (sills) in the aquifer (i.e. 0.4 and 4.0). As observed, results for the case with lowest uncertainty (sill=0.4) are tend to be deterministic and the CDFs are essentially in the form of vertical lines. But for the case of sill=4.0, one can observe the effect of increased spatial uncertainty in curvature in CDFs. However, the CDFs are less uncertain as opposed to the scenarios investigated in section 3.5.10.1, where more effective parameters (namely operational ones) were selected to be random. Page | 94 Figure 3-36: Comparing concentration CDFs for two extreme values of aquifer variogram sill: 0.4 (Left) and 4.0 (right). Making a comparison between source-to-receptor arrival time CDFs for two extreme values of variogram sills reveals that the case with smaller uncertainty (i.e. sill=0.4) is almost deterministic with an average arrival time value of 107 years. However, increasing the spatial uncertainty results in a CDF with average value of 110.4 years (Figure 3-37). Figure 3-37: Comparing source-to-receptor arrival time for two extreme values of variogram sill in the aquifer: 0.4 (Left) and 4.0 (right). Figure 3-38 shows the mean values of source-to-receptor arrival times for different cases of variogram sills. Mean arrival time has a direct proportionality with the variogram sill. Page | 95 Figure 3-38: Mean values of source-to-receptor arrival time for different cases of variogram sill And as the last EPM of concern, comparing the results for the risk measure (HI) shows that for the scenarios investigated the HI values are above 1 indicating a potential hazard. This, however, is expected as the parameters chosen for this set of simulation in this section are in a way to give the shortest arrival time which in turn is associated with a more critical situation. HI value for the case of sill=0.4 is almost constant (9.8). But by increasing the uncertainty in the aquifer permeability field, HI CDF takes a shape of an S-curve within the range of 8 to 10 (Figure 3-39). Figure 3-39: Comparing Hazard Index for two extreme values of variogram sill in the aquifer: 0.4 (Left) and 4.0 (right). HI value is almost deterministic for the lower value of sill. Figure 3-40 demonstrates the mean values of Hazard Index for different cases of variogram sill. As seen, increasing the spatial uncertainty within the aquifer results manifests itself as slightly Page | 96 lowered values of Hazard Index. The reason pertains to the presence of pressure gradient and flow inside the aquifer which makes advective as the dominant transport mechanism. Increasing uncertainty in the permeability field results in less effective transport in the aquifer, lower observed concentrations in the monitoring well, and smaller numbers of Hazard Index. Figure 3-40: Mean values of Hazard Index for different cases of variogram sill Page | 97 Conclusion Chapter 4: This study was proposed to investigate one of the possible risk pathways of groundwater contamination during a hydraulic fracturing operation, and assess the long-term consequences from a water quality perspective. It is assumed that the contamination plume originates from a leakage incident during a six-stage fracturing job. The leakage point is located on the vertical section of the injector with variable distances from the aquifer bottom. Furthermore, different geological and operational parameters are altered in order to explore effects of each parameter on the simulated concentration and arrival-time. Following are the conclusion remarks based upon the analysis and simulations conducted: 1. In an underground failure scenario during high-pressure injection, there is a likelihood of shallow groundwater contamination incident under special circumstances. Hydrogeological and operational parameters are the key factors in increasing the risk of drinking water ingestion when keeping the health-related parameters constant. 2. Failure scenario becomes more critical when working with higher values of water saturation in the overburden formation. So that, in the case with 3 . 0   w S the highest concentrations are observed within the aquifer and in the monitoring well. 3. Homogeneity in permeability field has proven to have more adverse effects at the beginning of the simulation, but from a given time (near 82 years) towards the end of simulation the case with highest anisotropy ratio (i.e. horizontal vertical k k /   ) takes the lead and shows the highest concentration values. In other words, an increase in vertical permeability is alleviating the situation in the long-term. Page | 98 4. The results of the sensitivity analysis show that, in general, the concentrations are more sensitive to the operational parameters than geological features. 5. Higher leakage rates, expectedly, result in higher concentrations in the aquifer. But concentration changes are more significant at lower leakage rates. For example, moving from a rate of 0.1 to 3 m 3 /d shows a sharper change in concentration values as compared to the move from 5 to 7 m 3 /d. 6. It is shown that even at the deepest leakage depth considered (namely 20 meters below the aquifer) concentrations within the aquifer cross the allowable limit. Thus, it is recommended to extend the surface casing of oil and gas wells to depths more than 20 meters below the treatable shallow water (critical depth). It is worth mentioning that several states have mandatory surface casing standards with which the oil and gas operators must comply. Oklahoma, for example, requires the setting depth of the surface casing as the greater value of 90 feet below the surface or 50 feet below the base of treatable water (Oklahoma administrative Code). Pennsylvania regulations also indicate that casing must be cemented with ASTM-approved cement to a minimum of 50 feet deeper than the deepest fresh groundwater (The Pennsylvania Code). The calculated Hazard Index for the leakage depth scenarios is in accord with the knowledge that 20 m below the bottom of fresh aquifer can be regarded as the critical depth. 7. The worst case scenario takes place when the water saturation of 3 . 0   w S and 1 . 0   . Results of the Monte Carlo simulation for this scenario show that concentration of the monitoring well does not go beyond the threshold during early stages of the simulation. But the highest calculated probability is obtained at the end of simulation, Page | 99 with value of 0.46. Also, on average it takes 276 years for the chemical concentration in the monitoring well to cross the allowable limit. The probability of Hazard Index passing the critical value of 1 is reported as 0.55 and the risk is measured as 0.45. 8. Increasing the spatial uncertainty (i.e. permeability) in the sand layer slightly affects the results. The main reason is due to the fact that diffusion is the governing mechanism for the most of the simulation period (except for the injection period where advection and dispersion are the main mechanisms) and diffusion is dependent on the molecular diffusion constant as well as the pore volume of the formation and not on the permeability. Spatially varying porosity should result in more sensitive scenarios. 9. Fixing the hydrogeological and operational parameters and altering the heterogeneous permeability field within the aquifer results in slight changes in the HI and concentration and a remarkable change in the arrival time for the set of parameters investigated. The effect of increasing the heterogeneity of the permeability field of the aquifer is highlighted through CDFs of the EPMs, where in the case of lowest heterogeneity (i.e. variogram sill of 0.4) CDF will be in the form a vertical line. In general, the extent of stratification of the aquifer is an important factor when performing studies similar to the current work. Small changes in the Hazard Index are yet distinguishable as the correlation variance is being altered. Increasing spatial uncertainty helps in mitigating the adverse situation by increasing the arrival-time at one hand and lowering the Hazard Index on the other hand. In conclusion, it is beneficial to discuss some of the limitations of the numerical simulations and models used for this study. First, the model used is isothermal and does not take into account the Page | 100 property changes typically associated with variations in temperature (e.g. density). Second, the chemical transport model does not consider the chemical adsorption on the rock surface of the porous media, nor does it include the chemical reaction; however, the chemical reaction might be neglected when working with stable chemicals. In general, excluding reaction and adsorption features is in agreement with the assumption of critical and conservative scenarios. Fracturing fluid in this research has been modelled as a Newtonian fluid (i.e. viscosity is constant with changes in shear rate and it changes only with temperature), however in reality proppant-laden slurries are Power Law fluids and show changes in viscosity when the shear rate alters (Economides & Martin 2007). The viscosity of the fracturing slurry is high in the beginning to be able to carry the proppants but as the proppants are placed in the fractures, the viscosity drops to increase the flowability of the returned waste (Economides & Martin 2007). Using real viscosity number instead of the one used in this research will result in lower concentrations and reduced associated risk values. Also, this study has focused on one chemical of concern and does not cover the effects of other chemical agents typically found in a fracturing slurry. This assumption together with introducing water density in the model affects the results. In real settings, the density of fracturing slurry is higher than water which makes the upward movement of the fluid phase harder. When fluid phase density is replaced with larger numbers the concentrations and risk values associated will most likely decline. As the final point, the simulation work assumes that the overall geology of the region remains constant during the course of simulation (500 years) and the plume is propagated within the constant stratification; in reality there might be significant geological processes happening within 500 years. As previously discussed, geological formations are of highly complex nature so that in a more accurate practice, one should incorporate the effects of heterogeneity and spatial Page | 101 uncertainty in geological formations more precisely (i.e. presences of fault, fissures, channeling, etc.) to draw conclusive results. In that case, the geological parameters will likely play more remarkable roles and the results will be more sensitive towards the geological features as opposed to the operational ones. In the end it is noteworthy to state that the output of this research is a framework/methodology which allows the site managers to make critical decisions efficiently and in a timely fashion in an event of a failure accident similar to the one discussed herein. While the methodology is applicable for different sites with various geological conditions, the findings of this research shall, by no means, be generalized. Meanwhile, the calculated arrival- times for the simulation sets clearly show that an underground breach can potentially have aftermaths which are not immediately observable on the surface. Page | 102 Future Research Direction Chapter 5: Long-term effects of a failure scenario during a hydraulic fracturing operation have been the subject of this study. The goal was to demonstrate some of the important geological and operational parameters contributing to the problem, and to define the sensitivity of the results for each parameter. The uncertainty of parameters was captured through application of a Monte Carlo simulation. As for prospective directions for this type of research, one could think of two major topics of -- upward flow migration from either risk pathways discussed previously, or to focus on waste handling and management specifically, given that there are critical issues ongoing regarding application of unlined pits to temporarily store returned fluid in oil and gas operations within the state of California and elsewhere. The scope of future work is hereby presented based on the preliminary findings of this study:  This study can be further extended by increasing the uncertainty dimensions. The uncertainty in health-related parameters (e.g. ED, EF, BW, etc.) must be given the utmost attention when building risk assessment frameworks as such. Hydrogeological features of the aquifer (other than permeability which was in part addressed in this study) shall be considered as well.  As a complementary to the current scenario, heterogeneity of the aquifer and the overburden formation must be modelled by honoring real data. Statistical methods must be applied to create random fields for parameters such as permeability while respecting the bivariate statistics of the field. In other words interdependencies between variables must be acknowledged by working with higher order moments and the correlations between each pair of the spatial points. On top, in order to wind up to meaningful results, a specific region must be chosen and real geological data must be used for this set of Page | 103 analyses. Appendix C, introduces a water contamination mapping software as an example of a research work which honors the real data. This software will help in pinpointing areas where groundwater quality has been impacted by underground injection operations.  Direct leakage into the aquifer is another risk pathway that must be investigated. Simulations need to be conducted to define the contamination movement within the aquifer, and towards specific monitoring well targets. Similar to current work, the intensity of the concentrations should be outlined in order to ascertain under what situations it crosses the reference dose.  Other risk pathways of the contaminant are candidates to be the subject of a future study. Presence of faults in the system and possible upward migration of fluid under reservoir pressure should also be considered. Additionally, the spillage of chemicals from onsite storage or from temporary pits used to store returned fluid and waste are among other topics that are currently garnering much attention. Page | 104 Appendix A: Defining the Normalized Threshold for the Chemical of Interest (TEPA) Hazardous Chemicals list (Gordalla et al. 2013) is the main source of data for this calculation work. Assuming chemicals are two weight percent of the fracturing fluid, gives the initial concentration of 22 kg/m 3 (22,000 mg/l) for the tracer used in model. The hazardous chemical which is used as the chemical of concern in this study is termed Tetraethylenepentamine which is a stable compound. According to the Material Safety Data Sheet, absorption through skin, dermal contact, eye contact, inhalation, and ingestion are possible routs of entry to the body. It is reported to have acute oral toxicity to animals, and mutagenic effects on bacteria and yeast and chronic effects on human being. Also, this chemical is very hazardous when in contact with skin, and when it is ingested. In case of any eye contact or inhalation, it is also considered a hazardous material. The agent is categorized as level 3 in health according to Hazardous Material Identification System (HMIS) (Figure A-1). Figure A-1: Hazardous Material Identification Rating for Tetraethylenepentamine (www.sciencelab.com) Data from (Gordalla et al. 2013) compares the concentrations of this chemical in fracturing fluid of two sites (both using the gel-based fluid system). Table A-1 shows the details of data including the allowable limit of Tetraethylenepentamine in drinking water. Page | 105 Table A-1: Tetraethylenepentamine concentrations in two fracturing fluids and drinking water reference concentration Tetraethylenepentamine Concentration (mg/l) Total Concentration of Chemicals (mg/l) Weight Percent of Chemical in Fluid Tetraethylenepentamine Reference Concentration (mg/l) Site A 736 30,912 3.00% 0.0003 Site B 173 14,288 1.13% The weight percent of chemicals used in the fluid in this study is assumed to be 2. Hence, the approximate concentration of Tetraethylenepentamine can be calculated by interpolating between given values for site A and B. The concentration will be given as 435 mg/l which includes 1.98% of the entire chemical concentration. Next, the reference concentration is scaled up for the tracer concentration (sum of all chemicals) which is reported by the simulator. Finally, the threshold value is computed by normalizing the reference value with respect to the initial concentration of 22,000 mg/l. Table A-2 shows the result. Table A-2: Calculation of threshold concentration Hazardous Chemical Weight percent of the Chemicals Reference Concentration (mg/l) Modified Reference Concentration (mg/l) Normalized Reference Concentration (Threshold) (C ref /C init ) Tetraethylenepentamine 1.98% 0.0003 0.015 6.9×10 -7 Page | 106 Appendix B: Calculation of Reference Dose for the Chemical of Interest (TEPA) As stated in section 0, for a non-carcinogenic chemical, hazard quotient (HQ) is defined as the exposure dose to a chemical over a period of time, divided by the daily exposure reference dose. The reference dose is the dose with zero likelihood of occurrence for adverse health effects (USEPA 2001). The chronic exposure hazard quotient is given by: RfD CDI HQ  (B-1) where CDI is the chronic daily intake [mg/(kg-d)] and RfD is the chronic reference dose [mg/(kg-day)] (i.e. calculated for the life-time period). The reference dose can be determined from its Drinking Water Equivalent Level (DWEL) which is calculated by assuming an average body weight of 70 kg and water consumption of 2 L/d (USEPA 2009). Page | 107 Appendix C: Water Contamination Mapping Tool Water Contamination Mapping Tool (WCMT) is a MATLAB-based software developed to help visualize the injection and water quality data in the first phase. Data are borrowed from three sources: water quality data from Groundwater Ambient Monitoring Assessment (GAMA) program run by the California State Water Resources Control Board, and FrackFocus/SkyTruth and Division of Oil, Gas and Geothermal Resources (DOGGR) as online databases for hydraulic fracturing operations. The ultimate goal is to find correlations between injectant concentrations and the constituents found in groundwater in the State of California. Three public data sets have been collected. Conventional indicators that are often measured to establish a general picture of the aquatic environment include pH, Total Dissolved Solids (TDS) and Suspended sediment. TDS has been selected as the water quality metric for this work. The TDS information is matched to some hydraulic fracturing jobs in different regions and in different years. The aim of this graphic user interface (GUI) is to find a way to assist resource assessment and field operations planning in the future. Continuously updating data in the software helps in meeting the demand and concerns in the industry related to monitoring of water contamination and correlate the cause of those events with field operations in the vicinity of each county. Different data sets used in this work are as follows: Fracfocus/SkyTruth: The FracFocus system is an online database developed to document disclosures for wells fractured after January 1, 2011 in pdf format. SkyTruth is another online resource with FracFocus data tabulated in spreadsheets. Well name, job date, injected volume, well coordinates and depth are data features collected from SkyTruth. Page | 108 Division of Oil, Gas and Geothermal Resources (DOGGR): DOGGR is the repository for oil, gas, and geothermal well information and publishes statistics on drilling, production, and injection in California. Hydraulic fractured well information from 2001 to 2011 and the county boundary coordinates were extracted from this public website. GeoTracker GAMA: This database has been used to obtain TDS data from 2001 to present. Vast volumes of data contain redundant information, and thus only the pertinent information required by the user with simple queries are selected. The information stored in the software includes date, sample coordinates, and magnitude. Data are indexed by 15 county names, and grouped by 6 oil and gas districts in California. Each set of data covers Total Dissolved Solids (TDS) and hydraulic fracturing job information from respective counties. It should be noted only those counties that both contain TDS and fracturing jobs information are chosen for analysis. Database also contains some other geographic features such as the location of all major and minor faults in districts, which may help to improve the understanding the correlation between water indicators and field operations. The software provides a platform to map counties in California and the detailed hydraulic fracturing wells in the field, simultaneously. It also allows the user to correlate data with the TDS observed with their magnitudes in that region. In order to retrieve pertinent data efficiently and clearly display attributes, the GUI provides options to tabulate and display corresponding data on the map. Figure C-1 shows the attribute categorization for this software: Page | 109 Figure C-1: Attribute categorization for water contamination mapping tool (WCMT) Figure C-2 shows the software interface. User can perform the following steps to generate the output with detailed information about field operations and TDS: 1- Select district of interest and then name of county (Figure C-3). 2- Select year of events (Figure C-4). 3- Select the data type (Figure C-5). 4- Select the features of the specific data type (Figure C-6). 5- Press Report tab to generate a map (Figure C-7) or press Table tab to generate tabulated output (Figure C-8). Page | 110 Figure C-2: WCMT interface Figure C-3: Selecting district of interest and the name of county Page | 111 Figure C-4: Selecting year of interest Figure C-5: Selecting type of data to be displayed Page | 112 Figure C-6: Selecting data features Figure C-7: Software output – volume and spatial distribution of injected water in Kern County Page | 113 Figure C-8: Software output – Tabulated results for injected water volume in Kern County User may also add TDS data and feature (Figure C-9) with numerical values shown on the map (Figure C-10). Figure C-9: Displaying TDS and hydraulic fracturing data at the same time Page | 114 Figure C-10: Showing numerical values of data on the map Figure C-11 shows a sample map for Kern county. Figure C-11: TDS and injected volume data in Kern County in year 2013 Page | 115 Preliminary analysis and result In the rest of this section, a few outputs from the software are selected to qualitatively investigate possible correlations between TDS and hydraulic fracturing jobs in field. In the selected maps, TDS magnitude is represented by circles. Color and size of the circles are both indicative of the magnitude of TDS. The TDS data is color-coded with 3000 mg/l chosen to be the threshold number. Aquifers with TDS values of 3000 mg/l and below are considered as protected aquifers and injection in those resources are not allowed (CCST & LBNL 2015). Diamond represents the hydraulic fracturing jobs, and its volume injected and depth are also presented by color and size. Figure C-11 to Figure C-16 demonstrate this relationship as shown in these figures. Figure C-12: TDS and injected volume data in Los Angeles County in year 2012 Page | 116 Figure C-13: TDS and injected volume data in Ventura County in year 2012 Figure C-14: TDS and injected volume data in Glenn County in year 2012 Page | 117 Figure C-15: TDS and injected volume data in Colusa County in year 2012 Figure C-16: TDS and injected volume data in Sutter County in year 2012 North and Northwest of Kern County show correlations between injected volume of water and the heightened TDS levels. Kern County have experienced degraded water quality as a result of mistakenly injecting in protected aquifers and also from other field activities such as discharge to unlined pits – Saline water, formation fluids, HF fluid (CCST & LBNL 2015). In another case, TDS samples collected in the Los Angeles County (Figure C-12) show correlations with hydraulic fracturing jobs in the South. There are regions towards the West of the county which show high levels of TDS without any injection works in vicinity. In the Ventura County Page | 118 (Figure C-13), potential correlation is found close to the Western and Eastern boundaries of the county. In the Glenn County (Figure C-14), Colusa County (Figure C-15) and Sutter County (Figure C-16) there seems to be no significant correlations between underground injection works and the TDS data. Finding correlation is the first step and more detailed monitoring plans and analysis must be followed to back-track the traces of contamination to ascertain the existence of a correlation. Focusing on one or a few contaminants that are solely used in injection operations could be a great practice. Many factors are still not considered herein. First, leakage of chemicals can happen through existing wells (old reservoirs with high well density). Second, the effectiveness of the fracturing job design has a direct impact on water quality. Shallow hydraulic fracturing presents a higher risk of groundwater contamination, which groundwater monitoring may not detect. Third, some geologic events may also change the spatial distribution of observed chemicals. Finally, historical observation (Before and after an event) must be given enough attention in analyses as such. Page | 119 Appendix D: List of Publications Journal Publications Jabbari, N., Aminzadeh, F., & de Barros, F. P. J. (Under review). Hydraulic Fracturing and the Environment: Risk Assessment for Groundwater Contamination from Well Casing Failure. Stochastic Environmental Research and Risk Assessment Jabbari, N., Aminzadeh, F., & de Barros, F. P. J., 2015. Assessing the Groundwater Contamination Potential from a Well in a Hydraulic Fracturing Operation. Journal of Sustainable Energy Engineering, 3(1), 66-79. Heinecke, J., Jabbari, N., & Meshkati, N., 2014. The Role of Human Factors Considerations and Safety Culture in the Safety of Hydraulic Fracturing (Fracking). Journal of Sustainable Energy Engineering, 2(2), 130-151. Bardet, J. P., Jesmani, M., & Jabbari, N., 2014. Permeability and compressibility of wax-coated sands. Géotechnique, 64(5), 341-350. Jesmani, M., Bardet, J. P., Jabbari, N., & Kamalzare, M., 2013. Hydraulic and mechanical properties of wax-coated sands. Journal of Central South University, 20(12), 3667-3675. Bardet, J.P., Jesmani, M., & Jabbari, N., 2011. Effects of Compaction on Shear Strength of Wax-Coated Sandy Soils. Electronic Journal of Geotechnical Engineering, 16, 451-461. Conference Publications and Book Chapters Jabbari, N., Ashayeri, C., & Meshkati, N., 2015. Leading Safety, Health, and Environmental Indicators in Hydraulic Fracturing. SPE Western Regional Meeting, Garden Grove, CA, 27-30 April 2015. Page | 120 Thimmisetty, C., Khodabakhshnejad, A., Jabbari, N., Aminzadeh, F., Ghanem, R., Rose, K., ... & Disenhof, C., 2015. Multiscale Stochastic Representation in High-Dimensional Data Using Gaussian Processes with Implicit Diffusion Metrics. In Dynamic Data-Driven Environmental Systems Science (pp. 157-166). Springer International Publishing. Rose, K., Aminzadeh, F., Ghanem, R., Thimmisetty, C., Jabbari, N., Khodabakhshnejad, A., Sim, L., Bauer, J., Mark-Moser, M., and Disenhof, C., 2014. Risks and impact assessment for deepwater and ultra-deepwater Gulf of Mexico resources. Offshore Technology Conference, Houston, TX, 5-8 May Presentations and Posters Jabbari, N., Aminzadeh, F., de Barros, F. P. J., & Jafarpour, B., 2014. Hydraulic fracturing and potential groundwater contamination risk. Presented at the 248th American Chemical Society National Meeting & Exposition, Paper Number: 315, San Francisco, CA, August 10-14 2014. Jabbari, N., Aminzadeh, F., Jafarpour, B., & de Barros, F. P. J., 2013. Hydraulic Fracturing and the Environment. American Geophysical Union Annual Meeting, Poster Number: H53B- 1413, San Francisco, CA, 9-13 December 2013. Jabbari, N., Aminzadeh, F., & de Barros, F. P. J., 2013. Hydraulic Fracturing: Points and Counter Points. Presented at Urban Water Institute, 20 th Annual Water Policy Conference, San Diego, CA, August 2013. Page | 121 Bibliography Adams, J. & Rowe, C., 2013. Differentiating applications of hydraulic fracturing. ISRM International Conference for Effective and …. Aminzadeh, F., Tiwari, A. & Walker, R., 2014. Correlation between Induced Seismic Events and Hydraulic Fracturing activities in California. In American Geophysical Union Fall Meeting. San Francisco. Arthur, J., Bohm, B. & Layne, M., 2008. Hydraulic Fracturing Considerations for Natural Gas Wells of the Marcellus Shale. The Ground Water Protection Council Annual Forum, Cincinnati, Ohio. Arthur, J.D. et al., 2009. Evaluating Implications of Hydraulic Fracturing in Shale Gas Reservoirs. SPE Americas E&P Environmental and Safety Conference, (March), pp.23–25. Ballio, F. & Guadagnini, A., 2004. Convergence assessment of numerical Monte Carlo simulations in groundwater hydrology. Water Resources Research, 40(4). de Barros, F.P.J., Ezzedine, S. & Rubin, Y., 2012. Impact of hydrogeological data on measures of uncertainty, site characterization and environmental performance metrics. Advances in Water Resources, 36, pp.51–63. de Barros, F.P.J. & Fiori, A., 2014. First-order based cumulative distribution function for solute concentration in heterogeneous aquifers: Theoretical analysis and implications for human health risk. Water Resources Research, pp.1–20. de Barros, F.P.J. & Rubin, Y., 2008. A risk-driven approach for subsurface site characterization. Water Resources Research, 44(1). de Barros, F.P.J. & Rubin, Y., 2011. Modelling of block-scale macrodispersion as a random function. Journal of Fluid Mechanics, 676, pp.514–545. Available at: <Go to ISI>://000292093300021. de Barros, F.P.J., Rubin, Y. & Maxwell, R.M., 2009. The concept of comparative information yield curves and its application to risk-based site characterization. Water Resources Research, 45(6). Bear, J., 1988. Dynamics of Fluids in Porous Media, Courier Dover Publications. Benekos, I.D., Shoemaker, C.A. & Stedinger, J.R., 2007. Probabilistic risk and uncertainty analysis for bioremediation of four chlorinated ethenes in groundwater. Stochastic Environmental Research and Risk Assessment, 21, pp.375–390. Page | 122 Brooks, R. & Corey, A., 1964. Hydraulic properties of porous media. Hydrology Papers, Civil Engineering Department, Colorado State University, (3). Browning, L. & Smith, J., 1993. Analysis of the rate of and reasons for injection well mechanical integrity test failure. SPE/EPA Exploration and Production …, (Vic). CCST & LBNL, 2015. An Independent Scientific Assessment of Well Stimulation in California Volume II, Charbeneau, R.J., 2006. Groundwater Hydraulics And Pollutant Transport, Waveland Pr Inc. Chen, Z., Huan, G. & Ma, Y., 2006. Computational Methods for Multiphase Flows in Porous Media. SIAM, 2. Cohen, H.A., Parratt, T. & Andrews, C.B., 2013. Potential contaminant pathways from hydraulically fractured shale to aquifers. Ground water, 51(3), pp.317–9; discussion 319– 21. Daneshy, A., 1990. Hydraulic Fracturing To Improve Production. SPE Tech 101, pp.14–17. DiGiulio, D. et al., 2011. Investigation of ground water contamination near Pavillion, Wyoming, Economides, M. & Martin, T., 2007. Modern fracturing: Enhancing natural gas production, Energy Tribune Publishing Inc., Houston. Ellsworth, W.L., 2013. Injection-induced earthquakes. Science (New York, N.Y.), 341(6142), p.1225942. Fisher, M. & Warpinski, N., 2012. Hydraulic-fracture-height growth: Real data. SPE Production & Operations, (November). Flewelling, S. a., Tymchak, M.P. & Warpinski, N., 2013. Hydraulic fracture height limits and fault interactions in tight oil and gas formations. Geophysical Research Letters, 40(14), pp.3602–3606. Freeze, R.A. & Cherry, J.A., 1979. Groundwater, Prentice Hall. Gassiat, C. et al., 2013. Hydraulic fracturing in faulted sedimentary basins: Numerical simulation of potential contamination of shallow aquifers over long time scales. Water Resources Research, 49(2). Van Genuchten, M.T., 1980. A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Science Society of America Journal, 44(5), pp.892–898. Page | 123 Glazer, Y.R. et al., 2014. Potential for Using Energy from Flared Gas for On-Site Hydraulic Fracturing Wastewater Treatment in Texas. Environmental Science & Technology Letters, 1, pp.300–304. Gleeson, T. et al., 2011. Mapping permeability over the surface of the Earth. Geophysical Research Letters, 38(2). Goebel, T. et al., 2014. A probabilistic assessment of waste water injection induced seismicity in central California. In American Geophysical Union Fall Meeting. San Francisco. Gordalla, B.C., Ewers, U. & Frimmel, F.H., 2013. Hydraulic fracturing: a toxicological threat for groundwater and drinking-water? Environmental Earth Sciences, 70(8), pp.3875–3893. Gregory, K.B., Vidic, R.D. & Dzombak, D. a., 2011. Water Management Challenges Associated with the Production of Shale Gas by Hydraulic Fracturing. Elements, 7(3), pp.181–186. Hammer, R. & VanBriesen, J., 2012. In fracking’s wake: new rules are needed to protect our health and environment from contaminated wastewater. Natural Resources Defense Council, (may). Healy, J. & Rubey, W., 1968. The Denver earthquakes. Science (New York, N.Y.), 161(3848), pp.1301–1310. Heinecke, J., Jabbari, N. & Meshkati, N., 2014. The Role of Human Factors Considerations and Safety Culture in the Safety of Hydraulic Fracturing (Fracking). Journal of Sustainable Energy Engineering, 2(2), pp.130–151. Hohn, M., 1999. Geostatistics and Petroleum Geology, Springer Science & Business Media. Holland, A., 2011. Examination of possibly induced seismicity from hydraulic fracturing in the Eola Field, Garvin County, Oklahoma. Oklahoma Gological Survey Open-File Report OF1- 2011. Holloway, M.D. & Rudd, O., 2013. Fracking: The Operations and Environmental Consequences of Hydraulic Fracturing, John Wiley & Sons. Howarth, R.W., Ingraffea, A. & Engelder, T., 2011. Natural gas: Should fracking stop? Nature, 477(7364), pp.271–5. International Programme on Chemical Safety, 2001. TETRAETHYLENEPENTAMINE. Screening Information Dataset. Available at: http://www.inchem.org/documents/sids/sids/Tetraethylenepentamine.pdf [Accessed March 12, 2015]. Page | 124 Jabbari, N., Aminzadeh, F. & de Barros, F.P.J., 2015. Assessing the Groundwater Contamination Potential from a Well in a Hydraulic Fracturing Operation. Journal of Sustainable Energy Engineering, 3(1), pp.66–79. Jabbari, N., Ashayeri, C. & Meshkati, N., 2015. Leading Safety, Health, and Environmental Indicators in Hydraulic Fracturing. In SPE Western Regional Meeting. Garden Grove, CA. Kargbo, D.M., Wilhelm, R.G. & Campbell, D.J., 2010. Natural gas plays in the Marcellus Shale: challenges and potential opportunities. Environmental science & technology, 44(15), pp.5679–84. King, G., 2012. Hydraulic Fracturing 101: What Every Representative, Environmentalist, Regulator, Reporter, Investor, University Researcher, Neighbor and Engineer Should Know. SPE Hydraulic Fracturing Technology Conference, pp.1–80. Kissinger, A. et al., 2013. Hydraulic fracturing in unconventional gas reservoirs: risks in the geological system, part 2. Environmental Earth Sciences, 70(8), pp.3855–3873. Leythaeuser, D., Schaefer, R. & Yukler, A., 1982. Role of diffusion in primary migration of hydrocarbons. AAPG bulletin, 4(April), pp.408–429. Llewellyn, G.T. et al., 2015. Evaluating a groundwater supply contamination incident attributed to Marcellus Shale gas development. Proceedings of the National Academy of Sciences, 112(20), pp.6325–6330. López, E., Schuhmacher, M. & Domingo, J.L., 2008. Human health risks of petroleum- contaminated groundwater. Environmental science and pollution research international, 15(3), pp.278–288. Ma, H.W., 2002. Stochastic multimedia risk assessment for a site with contaminated groundwater. Stochastic Environmental Research and Risk Assessment, 16, pp.464–478. Maxwell, R.M. & Kastenberg, W.E., 1999. A model for assessing and managing the risks of environmental lead emissions. Stochastic Environmental Research and Risk Assessment, 13(4), pp.231–250. Available at: <Go to ISI>://000082342900001. Maxwell, R.M. & Kastenberg, W.E., 1999. Stochastic environmental risk analysis: an integrated methodology for predicting cancer risk from contaminated groundwater. Stochastic Environmental Research and Risk Assessment (SERRA), 13(1-2), pp.27–47. Maxwell, R.M., Kastenberg, W.E. & Rubin, Y., 1999. A methodology to integrate site characterization information into groundwater-driven health risk assessment. Water Resources Research, 35(9), pp.2841–2855. McKenzie, L.M. et al., 2012. Human health risk assessment of air emissions from development of unconventional natural gas resources. Science of The Total Environment, 424, pp.79–87. Page | 125 Meyer, P., 1994. Monitoring network design to provide initial detection of groundwater contamination. Water Resources Research, 30(9), pp.2647–2659. Miskimins, J. et al., 2006. Hydraulic Fracturing Technology and Case Studies in Tight Gas Sands and Shales. Mohammed, G.A. et al., 2009. Comparison of Two Mathematical Models for 3D Groundwater Flow: Block-Centered Heads and Edge-Based Stream Functions. Transport in Porous Media, 79(3), pp.469–485. Moore, C.W. et al., 2014. Air Impacts of Increased Natural Gas Acquisition, Processing, and Use: A Critical Review. Environmental science & technology. Myers, T., 2012. Potential contaminant pathways from hydraulically fractured shale to aquifers. Ground water, 50(6), pp.872–82. Nelson, P., 1994. Permeability-porosity relationships in sedimentary rocks. The Log Analyst, 35, pp.38–62. Nicholson, C., 1992. Earthquakes associated with deep well activities- Comments and case histories. , pp.1079–1086. Nunn, J.A. & Meulbroek, P., 2002. Kilometer-scale upward migration of hydrocarbons in geopressured sediments by buoyancy-driven propagation of methane-filled fractures. AAPG Bulletin, 86(5), pp.907–918. ODNR, 2012. Preliminary Report on the North Star 1 Class II Injection Well and the Seismic Events in the Youngstown, Ohio, Area, Osborn, S.G. et al., 2011. Methane contamination of drinking water accompanying gas-well drilling and hydraulic fracturing. Proceedings of the National Academy of Sciences of the United States of America, 108(20), pp.8172–6. Palisch, T., Vincent, M. & Handren, P., 2008. Slickwater fracturing: food for thought. SPE Annual Technical Conference and Exhibition, (August), pp.327–344. Pater, C. De & Baisch, S., 2011. Geomechanical study of Bowland Shale seismicity, Rahm, D., 2011. Regulating hydraulic fracturing in shale gas plays: The case of Texas. Energy Policy, 39(5), pp.2974–2981. Reagan, M.T. et al., 2015. Numerical simulation of the environmental impact of hydraulic fracturing of tight/shale gas reservoirs on near-surface groundwater: Background, base cases, shallow reservoirs, short-term gas, and water transport. Water Resources Research, p.n/a–n/a. Page | 126 Rish, W., 2005. A probabilistic risk assessment of Class I hazardous waste injection wells. Underground injection science and technology. Advances in water sciences, (52), pp.93– 138. Rozell, D.J. & Reaven, S.J., 2012. Water pollution risk associated with natural gas extraction from the Marcellus Shale. Risk analysis : an official publication of the Society for Risk Analysis, 32(8), pp.1382–93. Rubin, Y., 2003. Applied stochastic hydrogeology., Oxford University Press. Rubin, Y. et al., 1999. The concept of block-effective macrodispersivity and a unified approach for grid-scale- and plume-scale-dependent transport. Journal of Fluid Mechanics, 395, pp.161–180. Available at: http://www.journals.cambridge.org/abstract_S0022112099005868. Saiers, J.E. & Barth, E., 2012. Potential contaminant pathways from hydraulically fractured shale aquifers. Ground water, 50(6), pp.826–8; discussion 828–30. Sanchez-Vila, X., Guadagnini, a. & Fernàndez-Garcia, D., 2009. Conditional Probability Density Functions of Concentrations for Mixing-Controlled Reactive Transport in Heterogeneous Aquifers. Mathematical Geosciences, 41, pp.323–351. Schlumberger, 2011. ECLIPSE: reservoir simulation software, Technical Description, Schmitz, R. & Cook, T., 1996. A risk based management approach to the problem of gas migration. SPE Health Safety and Environment in Oil and Gas Exploration and Production Conference., pp.705–715. Siirila, E.R. et al., 2012. A quantitative methodology to assess the risks to human health from CO 2 leakage into groundwater. Advances in Water Resources, 36, pp.146–164. Available at: http://dx.doi.org/10.1016/j.advwatres.2010.11.005. Spellman, F., 2012. Environmental Impacts of Hydraulic Fracturing, CRC Press. Tartakovsky, D.M., 2013. Assessment and management of risk in subsurface hydrology: A review and perspective. Advances in Water Resources, 51, pp.247–260. Tartakovsky, D.M., Dentz, M. & Lichtner, P.C., 2009. Probability density functions for advective-reactive transport with uncertain reaction rates. Water Resources Research, 45(7), pp.1–8. Timur, A., 1968. An Investigation Of Permeability Porosity & Residual Water Saturation Relationships For Sandstone Reservoirs. The Log Analyst, 9(04). USEIA, 2013. Annual Energy Outlook, Page | 127 USEIA, 2011. Review of Emerging Resources: U.S. Shale Gas and Shale Oil Plays, USEPA, 2015. Assessment of the Potential Impacts of Hydraulic Fracturing for Oil and Gas on Drinking Water Resources, USEPA, 2009. DRINKING WATER STANDARDS AND HEALTH ADVISORIES TABLE, USEPA, 2004. Evaluation of Impacts to Underground Sources of Drinking Water by Hydraulic Fracturing of Coalbed Methane Reservoirs Study - Appendix A, USEPA, 2010. Hydraulic Fracturing Research Study, USEPA, 2001. Risk Assessment Guidance for Superfund : Volume III - Part A , Process for Conducting Probabilistic Risk Assessment, USEPA, 1989. Risk Assessment Guidance for Superfund. Volume I: Human Health Evaluation Manual (Part A), USEPA, 2012. Study of the Potential Impacts of Hydraulic Fracturing on Drinking Water Resources, VanBriesen, J., Wilson, J. & Wang, Y., 2014. Management of Produced Water in Pennsylvania: 2010-2012. In SHALE ENERGY ENGINEERING. pp. 107–113. Vengosh, A. et al., 2013. The Effects of Shale Gas Exploration and Hydraulic Fracturing on the Quality of Water Resources in the United States. Procedia Earth and Planetary Science, 7, pp.863–866. Volz, C.D. et al., 2011. Contaminant Characterization of Effluent from Pennsylvania Brine Treatment Inc ., Josephine Facility : Implications for Disposal of Oil and Gas Flowback Fluids from Brine Treatment Plants, Wiseman, H., 2009. Untested waters : The rise of hydraulic fracturing in oil and gas production and the need to revisit regulation. Fordham Environmental Law Review, 20, pp.115–170. Zhou, W. & Arthur, R., 1994. Environmental Applications of ECLIPSE-A Multiphase Flow Simulator Used in the Oil and Gas Industry. WASTE MANAGEMENT-TUCSON- AMERICAN NUCLEAR SOCIETY, pp.895–900. Zoback, M., Kitasei, S. & Copithorne, B., 2010. Addressing the environmental risks from shale gas development, Worldwatch Institute. Electronic Resources: American Association of Petroleum Geologists Website, Accessed on February 16 2014, http://www.aapg.org/explorer/2013/12dec/historical1213.cfm Page | 128 FrackFocus Website, Accessed on March 10 2014, http://fracfocus.org/water-protection/drilling- usage Material Safety Data Sheet for Tetraethylenepentamine, http://www.sciencelab.com/msds.php?msdsId=9925207 Oklahoma administrative Code 165:10-3-1, Available at: ftp://204.87.70.98/occrules/Ruleshtm/forweb04newrules.htm Probublica Website, Accessed on February 17 2014, http://www.propublica.org/special/hydraulic-fracturing-national The Pennsylvania Code, Chapter 78:Oil and Gas Wells, Available at: http://www.pacode.com/secure/data/025/chapter78/chap78toc.html The United States Energy Information Administration, Available at: https://www.eia.gov/oil_gas/rpd/shale_gas.jpg

Abstract (if available)

Abstract Hydraulic fracturing is a well stimulation method used within the oil and gas industry to increase production. The recent developments of horizontal drilling, in combination with fracturing, have made this technique very popular for natural gas production from low-permeability formations such as tight sands and shales. Nevertheless, this technology is relatively new and the impact on the environment has been called into question, as the long-term effects are not thoroughly understood. From a technical standpoint, fracturing tight formations is accomplished by introducing highly pressurized fluid (fracturing slurry) into the formation to create fissures and cracks, which in turn increase the permeability. Fracturing slurry is a mixture of water, propping agent (mostly sand), and a combination of various toxic and benign chemicals which are responsible for reducing friction, inhibiting corrosion, and preventing scale. One of the main points of contention in fracturing concerns the specific chemical compounds that are being used. Debates regarding the risk versus reward ratio of fracturing are ongoing within the scientific community. Policy makers and the public have also expressed concerns about possible groundwater contamination during the hydraulic fracturing process. ❧ This study will take a close look at some of those concerns. Specifically, out of several risk pathways to the aquifer, leakage from a well casing during the injection periods has been addressed. Computer simulations of fluid flow in porous media and transport phenomena have been conducted and the sensitivity of results with respect to parameters of interest has been evaluated. In order to consider the risk to human health, Tetraethylenepentamine has been selected as the representative chemical with proven toxicity and with a pre-defined reference dose in literature, so that this agent can be used as a threshold to compare the results of various hypothetical scenarios. The uncertainty of the problem is addressed through a novel risk framework consisting of the numerical simulation of fluid flow and chemical transport in porous media and the Monte Carlo simulation. Hydrogeological and operational parameters are selected as random variables in the framework. This unique risk framework is developed to assess three environmental performance metrics: the solute concentration, the arrival times from source to receptor, and the ingestion hazard of the contaminated aquifer. The effect of parametric uncertainty in risk is also analyzed. ❧ The results show that risk strongly depends on water saturation, the anisotropy, and the heterogeneity of the permeability field. Initial water saturation of 0.3 and vertical to horizontal permeability ratio of 0.1 within the sand medium are associated with more adverse situation. Through the sensitivity analysis, it is indicated that leakage depth and leakage rate are more effective in controlling the measured risk value when compared to the hydrogeological properties. On the same line of thoughts, the critical depth for extending the surface casing below the depth of fresh water aquifer is introduced and calculated as an important result from the simulation sets. For the case of this study, the critical depth is 18 m below the base of freshwater. Results of the parametric uncertainty analysis for anisotropic homogenous scenarios show that on average it takes 276 years for the chemical to show concentrations above the threshold. Also, the hazard index is reported to be 0.45. However, introducing spatial uncertainty to the sand layer alleviates the situation by increasing the time and reducing the risk. Findings of this study can be applied by setting more stringent requirements for the well integrity to ensure that hydraulic fracturing is carried out in an environmentally safe and sound manner. ❧ In conclusion, this research will examine the likelihood of groundwater contamination caused by a specific failure in the well casing during a hydraulic fracturing job. However, in order to investigate the worst case scenario, the presented results are drawn from a specific situation in which the leakage is ongoing all over the injection phase of the operation. One implementation of the findings in this study is in enhancing oil and gas production standards, where it can be used to refine the surface casing requirements. In addition, the risk framework proposed by this study can be implemented in a real fracturing operation after enriching by real data. This practice will give site managers the power of making appropriate decisions in a timely manner in an event of an unwanted failure. It should be noted that the outlined concluding remarks are specifically related to the failure scenario, conceptual model, and geological settings discussed in this study. The results cannot and should not be generalized to other failure scenarios

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Groundwater contaminant transport predictions in a sustainable water management scenario

$Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells$

PDF

Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells

$An extended finite element method based modeling of hydraulic fracturing$

PDF

An extended finite element method based modeling of hydraulic fracturing

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

PDF

Efficient stochastic simulations of hydrogeological systems: from model complexity to data assimilation

$History matching for fractured reservoirs with Hough-transform-based parameterization$

PDF

History matching for fractured reservoirs with Hough-transform-based parameterization

$Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model$

PDF

Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model

PDF

Hybrid physics-based and data-driven computational approaches for multi-scale hydrologic systems under heterogeneity

PDF

Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection

PDF

Reproducible and rapid computational approaches for assessing contamination in natural aquifers

PDF

Modeling and simulation of multicomponent mass transfer in tight dual-porosity systems (unconventional)

PDF

Evaluating the role of energy system decarbonization and land cover properties on urban air quality in southern California

PDF

Stochastic and multiscale models for urban and natural ecology

PDF

Inverse modeling and uncertainty quantification of nonlinear flow in porous media models

PDF

Integrated reservoir characterization for unconventional reservoirs using seismic, microseismic and well log data

PDF

Electrokinetic transport of Cr(VI) and integration with zero-valent iron nanoparticle and microbial fuel cell technologies for aquifer remediation

PDF

Investigating the role of urban emission sources on redox-active PM compounds and the chemical analysis of the standardized diesel exhaust particles

PDF

A study of diffusive mass transfer in tight dual-porosity systems (unconventional)

PDF

From tugboats to trees: investigating the coupled systems of urban air pollution and meteorology

PDF

Geothermal resource assessment and reservoir modeling with an application to Kavaklidere geothermal field, Turkey

Asset Metadata

Creator Jabbari, Nima (author)

Core Title Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Environmental Engineering

Publication Date 04/18/2016

Defense Date 01/07/2016

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag environmental impacts,groundwater,hydraulic fracturing,OAI-PMH Harvest,risk assessment

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Aminzadeh, Fred (committee chair), De Barros, Felipe (committee member), Ghanem, Roger (committee member), Sanders, Kelly (committee member)

Creator Email jabbari.nima@gmail.com,jabbari@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-231045

Unique identifier UC11278920

Identifier etd-JabbariNim-4273.pdf (filename),usctheses-c40-231045 (legacy record id)

Legacy Identifier etd-JabbariNim-4273.pdf

Dmrecord 231045

Document Type Dissertation

Format application/pdf (imt)

Rights Jabbari, Nima

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA