Modeling of methane migration from gas wellbores into shallow groundwater at basin scale

Methane contamination of drinking water resources is one of the major concerns associated with unconventional gas development. This study assesses the potential contamination of shallow groundwater via methane migration from a leaky natural gas well through overburden rocks, following hydraulic fracturing. A two-dimensional, two-phase, two-component numerical model is employed to simulate methane and brine upward migration toward shallow groundwater in a generic sedimentary basin. A sensitivity analysis is conducted to examine the influence of methane solubility, capillary pressure–saturation relationship parameters and residual water saturation of overburden rocks, gas leakage rate from the well, tilted formations, and low-permeability sediments (i.e., claystones) on the transport of fluids. Results show that the presence of lithological barriers is the most important factor controlling the temporal–spatial distribution of methane in the subsurface and the arrival time to shallow groundwater. A pulse of high leakage rate is required for early manifestation of methane in groundwater wells. Simulations reveal that the presence of tilted features could further explain fast-growing methane contamination and extensive lateral spreading reported in field studies.


Introduction
The global demand for energy is continuously increasing around the world and the International Energy Agency (IEA) anticipated more than 25% of growth by 2040 (IEA 2018). Natural gas production from unconventional resources is considered as a promising future source for energy supply as a bridge fuel toward a low-carbon energy system (Brown et al. 2009;McGlade et al. 2013). The surge for unconventional gas development has been projected to have impact on drinking water resources (Osborn et al. 2011;Sauter et al. 2012;Vengosh et al. 2014). The most frequent and severest threats are likely associated with (1) water acquisition in dry areas, (2) leaks and spills of fracturing fluids and wastewater (i.e., flowback and produced water) at the surface, (3) discharge of inadequately treated wastewater to the environment, and (4) migration of liquids (i.e., fracturing fluids and brine) and gas from deep hydrocarbon reservoirs into shallow aquifers (Melchers 2009;Rudolph et al. 2010;Sauter et al. 2012;Kissinger et al. 2013;U.S. EPA 2015;Rice et al. 2018a;Schout et al. 2020). Of the potential mechanisms identified, the focus of this study is on fugitive gas and brine migration from a leaky natural gas well through overburden rocks towards shallow groundwater.
Natural gas is mainly composed of methane with a small amount of carbon dioxide, oxygen, nitrogen, and hydrogen sulfide. Methane in stream flows and groundwater is either of biogenic and thermogenic origin. The origin of methane can be typically assessed by geological, hydrological, and geochemical analysis (Jackson et al. 2013;Darrah et al. 2014;Siegel et al. 2015). Biogenic methane is produced by microbiological processes from in situ fermentation or reduction of carbon dioxide. Thermogenic methane is derived from the upward migration of natural gas associated with deep hydrocarbon formations (Gorody 2012;Darrah et al. 2015). Moreover, thermogenic methane has been observed to reach 1 3 432 Page 2 of 16 shallow groundwater systems through anthropogenic activities such as hydraulic fracturing (Warner et al. 2012;Darrah et al. 2015). Methane does not react until dissolved in groundwater, in which case methane is rapidly being oxidized given respective redox conditions (Sherwood et al. 2016;Molofsky et al. 2016). Although methane is commonly regarded as a non-toxic substance, elevated levels of dissolved methane concentrations could change groundwater geochemistry, potentially leading to water-quality degradation (Duncan 2015;Hendry et al. 2017).
The presence of thermogenic methane in groundwater in the vicinity of oil and gas operations is most commonly associated with leakage from hydrocarbon wells (Jackson et al. 2013;Darrah et al. 2015;Rivard et al. 2019). Multiple studies have shown that leaked methane from the deep subsurface migrates upward by the buoyancy resulting from density contrasts between gas and brine (Rice et al. 2018a). Natural and anthropogenic permeable pathways such as leaky oil and gas abandoned wells, fault zones, and extensive fracture systems could facilitate the vertical movement of methane to shallow depths (Kissinger et al. 2013;Vengosh et al. 2014;Steelman et al. 2017;Yudhowijoyo et al. 2018;Tatomir et al. 2018;Riddick et al. 2019).
The quantitative analysis of the extent of methane leakage is complex even in simple geologic settings Cahill et al. 2017). Along with observational field studies, such as controlled methane-release experiments , numerical modeling of subsurface flow and transport of methane plays an essential role in advancing conceptual models to be able to better evaluate the relationship between groundwater quality and hydrocarbon development. Numerical modeling could be used to address unresolved questions associated with the observational studies on fluid flow and mass transport to groundwater. A major challenge associated with numerical modeling of methane migration between deep hydrocarbon reservoirs and shallow aquifers is due to the complex nature of multiphase, multicomponent flow, and transport phenomena involved in methane leakage through a variety of geological materials, wellbores, faults, fractures, and low-permeability rocks.
The major focus of previous numerical studies was on the rapid migration of methane to groundwater through leaky wellbores, faults, and fractures (e.g., Kissinger et al. 2013;Reagan et al. 2015;Schwartz 2015;Roy et al. 2016;Rice et al. 2018b;Moortgat et al. 2018;Soltanian et al. 2018;Nowamooz et al. 2018;D'Aniello et al. 2019), while less effort was directed at understanding the slow migration through overburden rocks (Nowamooz et al. 2015) as well as the influence of methane solubility on transport (Soltanian et al. 2018;Schout et al. 2020). In this study, we investigate the migration of free-and dissolved-phase methane from a leaky natural gas well to groundwater through overburden sediments at the basin scale. We specifically focus on the role of low-permeability sediments, tilted formations, and methane leakage rates and periods on the transport of methane to shallow subsurface.
We start with a generic conceptual model derived from an FEP knowledge database and the expert ranking of the key features, events, and processes that generate failure scenarios at basin scale (Tatomir et al. 2015. From the conceptual model, we build two-dimensional two-phase, twocomponent numerical models to study the flow and transport of methane and brine from a leaky gas well into groundwater for 25 years. The horizontal model domain has 1800 m depth and 10 km width. First, we evaluate the impact of methane solubility on fluid migration dynamics and the leakage rates into shallow subsurface. Then, we examine the influence of (1) the parameters describing capillarity pressure and relative permeability saturation of shallow and overburden rocks, (2) the low-permeability sediments, (3) the tilted formations (up to 3°), and (4) the methane leakage rate from the gas well. Finally, we compare the modeling results with the findings of observational studies.

Conceptual model for methane and brine migration at the basin scale
We build a generic basin-scale two-dimensional numerical model, Fig. 1, comprising three horizontal layers. The conceptual model is built based on permeability and stratigraphy data compilations of North America and Europe implemented by  and Potter van Loon (2016). We applied median values reported for the domain length, recharge, claystone thickness, and claystone depth. A thin (21 m) continuous claystone is assumed at shallow depth (171 m) overlain by coarse sediments. The model spans across a domain of 10 km length and 1.8 km depth.
A large section of the top boundary, i.e., 95% of the model length, is defined as a constant recharge boundary intersecting a discharge boundary (5% of the model length). The groundwater recharge rate amounts to 0.05 m/year. The discharge is modeled as a constant hydraulic head boundary . The lateral sides of the domain are set to no-flow boundaries.
It is common in the oil and gas industry to temporarily shut-in a well because of operational and economic reasons. The pressure in the gas reservoir does not decrease during the shut-in period. For this type of condition, should the well annulus intersect with a fracture/fracture system connecting the gas reservoir with the overburden, considerable leakage may occur (Reagan et al. 2015;Moortgat et al. 2018;Tatomir et al. 2018). The leakage rate depends on various factors such as reservoir pressure, reservoir gas saturation, in situ, and relative permeability (Reagan et al. 2015;Nowamooz et al. 2015). Kissinger et al. (2013) modeled constant leakage of methane from the reservoir into the overburden for 100 years. The authors emphasized that the constant escape of gas is regarded as "a conservative" assumption. Reagan et al. (2015) suggested that incidents of gas leakage from the reservoir are likely to be limited in time. Moortgat et al. (2018) considered a high methane inflow rate during the shut-in period, followed by a lower rate during the production.
The conceptual model considered here assumes fracture propagation into overlying formations during a hydraulic fracturing treatment (Kissinger et al. 2013;Reagan et al. 2015;Moortgat et al. 2018;Tatomir et al. 2018). The model assumes that during shut-in gas leaks into the overburden at a high flow rate, followed by a much lower rate once the production is resumed. We applied leakage rates published in Moortgat et al. (2018) as, similar to this study, they assumed methane leakage from a natural gas well into overlying formations. However, the geological setting, rock types, and tectonics are different from those considered here.
The leakage is assumed to occur at the top of the gas reservoir through the fracture system intersecting the overburden; thus, instead of explicitly considering the reservoir in the model, it is included as a boundary condition. The gas entering the domain is assumed to be pure methane and the presence of other thermogenic gases is neglected (Nowamooz et al. 2015;Roy et al. 2016;Rice et al. 2018a). The leakage source measures 100 m in length and it is situated at the bottom of the domain at 40% of the model length ( 3950 ≤ x ≤ 4050 ) for all subsequent simulations. The remaining of the lower boundary is set to no-flow condition. The methane leakage rate is estimated based on the production rate from the borehole. A typical production rate of 300 MSCFD (million standard ft 3 /day) (Moortgat et al. 2018) (35,680 m 3 /day at reservoir condition) is assumed for the base-case model. Methane leaks into the overburden at a high inflow rate of 8920 m 3 /day, i.e., 25% of the production during the 7 days of the shut-in period, followed by an inflow rate of 892 m 3 /day, i.e., 10% of the initial leakage rate, for the next 2 years. The upward migration of methane and brine is monitored for 25 years of simulation time.
The low-permeability unit, i.e., claystone, is a laterally extensive single layer overlying the entire basin which might contain discontinuities. Moreover, multi-layered systems where low-permeability units are interbedded within permeable sediments may exist. Each field configuration could lead to a different contamination distribution in the subsurface based on its unique stratigraphic, geophysical, and hydrodynamic characteristics (Sauter et al. 2012;Lange et al. 2013). Herein, we first assume a single continuous claystone across the entire model length to obtain the firstorder information on the influence of low-permeability sediments on methane migration. We then consider different configurations to explore the impact of claystone integrity and multi-layered clay systems.

Numerical modeling approach
The dominant physical and chemical processes for fluid migration are advection (viscous, buoyant, and capillary), dissolution, and diffusion (Nowamooz et al. 2015). Methane sorption is assumed to be negligible (Rice et al. 2018b). The flow and transport of fluids are simulated using both a miscible flow model and an immiscible flow model. The miscible flow model considers advection, diffusion, and dissolution processes, while the immiscible model considers migration via advection and diffusion.
All simulation runs are conducted with the open-source numerical simulator DuMu X (Flemisch et al. 2011). For spatial discretization, the box method, a vertex (or node)centered finite volume method based on a finite-element grid, is used. The grid resolution equals 500 cells in the x-direction and 180 cells in the y-direction, corresponding Generalized conceptual model of methane and brine migration from the gas reservoir into the overlying formations for the base-case model. The low-permeability layer, i.e., claystone, is indicated in gray and permeable sediments in light brown. The black dashed lines rep-resent the monitoring locations at distances Δ = 2000 m, 3000 m, and 4000 m from the leaky gas well. The leakage area is presented in yellow. Groundwater flow is from right to left, and the discharge area is modeled as a Dirichlet boundary condition to finite volume element length of 20 m in x-direction, and 10 m in y-direction. A fully implicit Euler time discretization scheme using the Newton method is applied to solve the linearized system of non-linear partial differential equations (Helmig 1997). The first time-step size is set to 0.01 s and the maximum time-step is 1e6 s. The time-step size is controlled by the number of iterations required by the Newton method to achieve convergence for the last time integration (Flemisch et al. 2011;Koch et al. 2020).
We use 2P2C (two-phase, two-component) and 2P (twophase) modules for simulating isothermal miscible and immiscible two-phase flow and transport of fluids in porous media, respectively. The binary diffusion and the diffusive mass fluxes in each phase are computed according to Fick's law. The mass balance equation for a two-phase two-component model can be formulated for each k as follows (Helmig 1997;Dietrich et al. 2005): where the index α represents the wetting and non-wetting phases. The index κ represents the components brine and methane. φ is the effective porosity and is the density of phase α. X k is the mass fraction of component κ in the phase α, S is the saturation, is the dynamic viscosity, and k r is the relative permeability. K is the intrinsic permeability tensor, P is the phase pressure, g is the gravitational acceleration, and D k Pm, is the effective diffusion coefficient of the porous medium.
Methane can exist as a free-phase gas with a limited amount of water vapor or dissolved in the aqueous phase. The amount of water in the gas phase based on Duan et al. (1992a). Methane solubility varies over a wide range of temperatures, pressures, and brine compositions (Duan et al. 1992a;Duan and Mao 2006). Methane dissolution in the aqueous phase is assumed to occur at thermodynamic equilibrium and gas-phase saturations were simultaneously updated at each time-step to consider the mass loss to brine. Assuming methane at equilibrium with respect to water, methane solubility (mg/L) can be computed following the Duan et al. (1992a)

approach as in:
where m is the molality of CH 4 or salts in the liquid phase and x is the composition in the vapor phase. P is the total pressure and R is the universal gas constant. l(0) standard chemical potential of CH 4 in the liquid phase and ∅ is the fugacity coefficient. The index a and c are anion and cation, respectively. CH 4 −ion and CH 4 −c−a are interaction parameters.
RT and the interaction parameters are calculated by Pitzer et al. (1984) method (Duan et al. 1992a).

Material properties
The variations in porosity and permeability are typically driven by changes in mineralogical composition, i.e., mineral grains and cement, depositional facies, burial depth, and post-depositional processes (Weibel et al. 2012). The mineral composition of the claystone comprises 50% kaolinite and 50% illite (Luijendijk and Gleeson 2015). Following Gleeson et al. (2011), we assumed that the shallow and overburden sediments contain 90% sand and 10% clay and obtained realistic permeability values. The permeability of sand is calculated using the Kozeny-Carman equation (Kozeny 1927;Carman 1937). The claystone permeability is determined using an empirical permeability-void ratio equation (Luijendijk and Gleeson 2015). The permeability of sediments is calculated following Luijendijk and Gleeson (2015). The porosity is calculated using lithology dependent exponential porosity-depth coefficients (Athy 1930;Sclater and Christie 1980). The estimated porosity and vertical permeability of claystone are 0.05 and 1.1 × 10 -20 m 2 , respectively (Fig. 2). Permeability and porosity of shallow and overburden sediments in the y-direction vary approximately from 6.8 × 10 -14 to 2.6 × 10 -15 m 2 and from 0.3 to 0.12, respectively, as shown in Fig. 2. The anisotropy ratio (K h /K v ) is set to 100 in the claystone and 10 in the surrounding sediments. The temperature at the top of the domain is 10 °C, with a vertical temperature gradient of 0.03 °C/m. The domain is initially water-saturated. The residual water saturation is set to 20% (Nowamooz et al. 2015;Rice et al. 2018b). In all simulations, we modeled a drainage-like process, because the gas saturation increases continuously with time. As gas cannot be trapped during the drainage, the residual gas saturation is assumed to be 0.01% (Rice et al. 2018b).
Methane is the predominant component in the gas mixture and small amounts of other components (e.g., N 2 , CO 2 , C 2 , and C 3 ) have a negligible influence on the density and viscosity of the gas phase (Nowamooz et al. 2015). For simplification, the model considers the leaked gas as pure methane. Brine is modeled as a pseudo-component and represents pure water with a given salinity. Brine density and viscosity are calculated following Batzle and Wang (1992) equations. Methane density and viscosity are calculated based on Duan et al. (1992b, c) and (Daubert et al. 1998), respectively. The diffusion coefficient of water in the gas phase is calculated based on Fuller et al. (1966). Since the empirical equations for estimating the diffusion coefficient in infinite solution presented in Reid et al. (1987) all show a linear dependency on temperature, we simply scale the experimentally obtained diffusion coefficient of the dissolved methane in the water phase given in Kobayashi (2004) by temperature. Kobayashi (2004) suggested the diffusion coefficient of 3.55 × 10 -9 m 2 /s. As methane leakage rates reported by Moortgat et al. (2018) are specified for surface conditions, a gas formation volume factor (B g ) of 0.0042 ft 3 /SCF is used to convert the values to reservoir conditions. The B g is calculated using the averaged values reported by Edwards and Celia (2018). The model property values are presented in Table 1. The capillary pressure and the relative permeability are calculated as a function of saturation using Brooks and Corey (1964) equations as follows: where P d is the entry pressure of the porous medium. The relative permeability functions, K r (S) , can be calculated by the following equations: The effective saturation S e is defined using the following equation. S wr and S gr are irreducible saturation of liquid and gas phases, respectively:

Simulation strategy
The base-case models explore miscible and immiscible flow and transport of fluids to shallow layers for a period of 25 years. The temporal and spatial behavior of the plume in the subsurface and the change in volumetric flow rates to shallow groundwater of methane are observed. Moreover, methane and brine flow rates at various distances from the leaky gas well (ΔL = 2000, 3000, and 4000 m) are monitored. These locations can be interpreted as intersections of the aquifer with conductive permeable pathways such as faults and fractures, i.e., preferential flow paths for the transport of fluids in the upper crust. This choice allowed us to track the methane plume at various distances from the gas wellbore in basins with or without vertical permeable pathways. A series of sensitivity analyses is designed to examine the roles of different geometric configurations and parameter combinations on the migration of fluids. Parameter selection is based on comprehensive features, events, and processes (FEPs) procedure highlighting the relevant importance of the key factors influencing the upward movement of contaminants . First, we explore the sensitivity of hydrodynamic properties ( P d , λ, and S wr ) of shallow and overburden sediments. The displacement or entry pressure (P d ) of the porous medium changes between 2 and 20 kPa and is set to 10 kPa in the base-case model. A large pore size distribution index (λ) indicates a comparatively uniform pore size distribution, while small values of λ describe a non-uniform distribution (Brooks and Corey 1964). λ varies We include simulations for sedimentary basins with dip angles of 0º, 1º, and 3º. We rotate the model coordinate system to tilt the geological geometric configuration. We keep the top of the domain at the same depth and increase the maximum depth by increasing the dip. Thereafter, the impact of methane leakage rate from the wellbore is explored. As an alternative to the high leakage with short duration, two low constant gas inflow rates of 60 and 120 m 3 /day for uninterrupted 25 years are considered.
The sensitivity analysis then considers a discontinuity, 10% of the model length, in the middle of the claystone to examine the role of low-permeability rocks on fluid flow and mass transport. Further simulations include the integration of a multi-layered claystone system with interbedded permeable rocks at shallow depth. We assume two different field settings by varying the geometry of the clay system, consisting of four clays with random horizontal position, depth, and integrity. All clays are embedded between 150 and 350 m depth and have the same thickness as the basecase value, 21 m. This sensitivity analysis is presented to examine how the migration of methane from a deep gas reservoir is impacted by sequences of low-permeability rocks.

Base-case model
The gas saturation and pressure profiles of methane leakage using the miscible flow model are shown in Figs. 3 and 4, respectively. In this model, methane exists as free and dissolved-phase gas in the aqueous phase. Methane migrates upward predominantly by buoyancy and reach the claystone after approximately 695 days. The strong vertical buoyancy controls the direction of methane migration; thus, the lateral spread is limited. The tube-shaped flow of gas is partially due to the consideration of a homogenous model domain. The diameter of the plume around the leakage source is approximately 320 m after 2 years and remains constant. The claystone overlying the overburden sediments constitutes an effective barrier to the upward movement of fluids. For this type of configuration, once methane reaches the claystone layer, the gas spreads laterally along the horizontally stratified layer due to the sustained pressure. As illustrated in Fig. 5, an important volume of gas phase accumulates at the claystone base after 25 years reaching eventually the aquifer. Permeable pathways, such as fractures, faults, and abandoned wells through barrier rocks, play an important role in fluid migration in the upper crust. In the context of hydraulic fracturing operations, the presence of these flow paths in the proximity of a leaky gas well could lead to the early occurrence of methane in groundwater (Darrah et al. 2014;Cahill et al. 2017;Moortgat et al. 2018).
For the time 25 years after the leakage, the dissolved methane (mg/L) in the aqueous phase is presented in Fig. 5b. The density of the gas phase and the composition of the aqueous phase strongly depend on pressure and temperature. Methane solubility in the aqueous phase is higher in deep formations due to the increase in pressure (Eq. 2), and thus, more methane is dissolved in brine before the appearance of the gas phase (Duan et al. 1992a;Duan and Mao 2006). Following the commencement of leakage, free-phase methane moves upward by buoyancy and progressively dissolves in brine. Dissolved methane is gradually transported by advection and diffusion. The high buoyancy force from density contrasts between methane (density > 110 kg/m 3 ) and brine (density > 1100 kg/m 3 ) controls the long-term upward migration of free-phase methane. As illustrated by the pressure profiles (Fig. 4), the hydrostatic pressure is lower at shallow depth. The reduction of the hydrostatic pressure leads to a significant decrease in the density of the gas phase, and the release of some of the dissolved methane into the gas phase. Accordingly, more highly buoyant free-phase methane appears at shallow depths with the consequence of increasing the upward flow rate. For such conditions, large amounts of free-phase methane migrate to the groundwater system within a relatively short time.
The gas-phase saturation and pressure profiles of the immiscible flow model are shown in Figs. 6 and 7, respectively. Methane reaches the base of the claystone after approximately 660 days. The methane accumulation at the claystone is followed by the extensive horizontal spreading. Similar to the methane migration in the miscible flow model set-up, the diameter of the plume around the leakage source measures approximately 320 m after 2 years and remains unchanged. In both miscible and immiscible flow models, the lateral spread of the plume before reaching the claystone is limited, partially due to the assumption of a homogenous domain.
The mass flow rates of methane and brine are obtained from the total non-wetting phase and wetting phase flow rates calculated from advection, diffusion, and dissolution processes. Figure 8a, b represents methane and brine volumetric flow rates at 2000, 3000, and 4000 m distances to the gas well, respectively. Methane reaches the monitoring locations during the extensive lateral spreading at the    claystone base. The methane flow rate increases rapidly to its maximum value and then shows a continuous decrease with time (Fig. 8a). The flow rates of fluids corresponding to the miscible model are slightly lower compared to the immiscible model, indicating the effect of gas dissolution on transport. Similar to the observations of Tatomir et al. (2018), the methane flow rate is higher where the conductive feature, e.g., fault, is closer to the leakage location. However, leakage continues to be active for much longer time periods when the conductive features are further located.
The pressure gradient created by the gas influx causes brine to slowly migrate in the subsurface (Fig. 8b). Brine flow rate at the measuring locations increases during the leakage and declines sharply after the drop-in pressure. The first peak value corresponds to the end of the shut-in period and the second peak corresponds to the end of leakage. Brine flow changes in the direction of regional groundwater flow after approximately 5 years due to the combined influence of regional hydraulic gradient, cease of gas leakage, and negative buoyancy. Although brine flows towards the monitoring locations (i.e., highly permeable pathways), the probability of deep brine reaching shallow groundwater is low (Osborn et al. 2011;Gassiat et al. 2013;Pfunt et al. 2016;Taherdangkoo et al. 2019) due to the absence of significant driving forces.
As illustrated by Figs. 3, 4, 5, 6, 7, and 8, using miscible and immiscible flow models results in slightly different methane plume sizes and travel times to groundwater. This is due to the low solubility of methane in the ambient water, particularly at shallow depths. The influence of methane solubility on gas transport becomes less significant once the plume spreads at shallow groundwater.

Sensitivity analysis
The sensitivity analysis is conducted employing the miscible flow model. First, the influence of key parameters including Brooks-Corey parameters, residual water saturation, tilted geometries, and methane leakage rates are examined. Then, the impacts of the claystone integrity and multi-layered claystone systems on the spatial distribution of the methane plume in the subsurface and the arrival time to shallow groundwater are discussed.

Geological and hydrogeological parameters
The impacts of hydrodynamic properties (i.e., P d , λ, and S wr ) of the shallow and overburden sediments on the transport of methane are illustrated in Fig. 9a-c. The entry pressure in the capillary pressure function is inversely correlated with methane flow rate at the monitoring locations. A higher entry pressure requires a higher non-wetting phase (i.e., gas phase) pressure to displace the ambient water resulting in the reduction of methane flow rates, as shown in Fig. 9a. The change of entry pressure has a substantial influence on the flow rate when P d is varied over a large range.
The pore size distribution index (λ) is a shape parameter in the capillary pressure-relative permeability-saturation relationships, and has influence on the phase saturation. An inverse relationship between λ and methane flow rate is observed as an increase in the λ value increases the degree of saturation of the water phase. Varying λ within a small range has negligible impacts on arrival times and flow rates of methane, while the influence of varying λ over a wider range is noticeable.  Residual water saturation (S wr ) is important with respect to its influence on the relative permeability. A positive correlation between S wr and methane flow rates is observed as an increase in the value of S wr decreasing the amount of gas phase trapped in the pore spaces. The effects of varying S wr on breakthrough time and volume of methane reaching monitoring locations are more substantial in our models in contrast to variations in P d and λ.
In a sedimentary basin with a larger dip angle (dip ≥ 3 • ), the lateral buoyancy component further contributes to the horizontal spreading of the plume. Methane moves up-dip along tilted layers and spreads over a large distance from the leaky wellbore. The lateral transport distances of methane are significantly higher even for only minor tilt (1°) as compared to a horizontal domain. In the formation with 3º tilt in the x-direction, methane reaches the monitoring location 4000 m away from the gas wellbore within 7.2 years, but it does not spread that far in the horizontal dimension (Fig. 9d). The results suggest that the presence of tilted formations/layers can encourage the lateral spread of the plume; thus, methane can be manifested at monitoring locations even at distances of several kilometers from a leaky gas well.

Methane leakage rate from the wellbore
The pulse of high leakage rates assumed in the base-case model was compared with the continuous slow leakage of methane. The lowest gas inflow rate, i.e., 60 m 3 /day, results in a delayed arrival of fluids at the monitoring locations in contrast to other cases, clearly showing less spatial spreading. The plume size is, as expected, larger with longer exposure to the leakage. Persistent long-term leakage leads to high methane concentrations in shallow layers, which induces higher contamination risks to groundwater aquifers. A pulse of high leakage rate is required to generate a larger plume and methane manifestation in shallow groundwater within short times. As illustrated in Fig. 9e, leakage rate and leakage period have significant impacts on the propagation of the plume in the overburden and accumulated methane volume in shallow groundwater.
Varying gas inflow rates in a tilted formation (~ 1°) (Fig. 9f) has significant influence on the temporal and spatial distribution of the plume in the overburden. The persistent leakage in a tilted formation leads to the extensive spreading in the direction of the tilt enhancing the risk to groundwater within a short-time scale, particularly with the presence of permeable pathways within this range. Under these conditions, free-and dissolved-phase methane could affect a larger area and millions of kilograms of methane are most likely to be transported to the aquifer. This scenario indicates that gas leakage in a tilted formation pose additional risk to shallow groundwater, partly due to the increase of lateral methane migration.

Claystone layer characteristics
We assumed a discontinuity in the center of the claystone at the location 4500 ≤ x ≤ 5500 m , to assess the influence of barrier layer integrity on the vertical movement of fluids. Methane migrates upward from the leakage source and accumulates at the claystone base. The trapped methane is compressed with the increase of pressure and tends to move upward through permeable pathways, e.g., fractures, and discontinuities. Once methane reaches the discontinuity, it preferentially migrates vertically along sediments with higher permeabilities (Fig. 10a). The low pressure and temperature at shallow depths have strong impacts on the density of the free-phase gas and its solubility in the aqueous phase; subsequently, fast vertical migration above the claystone formation slows the lateral spreading. Methane enters shallow groundwater after 2.4 years at a flow rate of 26.15 m 3 /day and peaks at 3.3 years at a flow rate of 26.15 m 3 /day (Fig. 11).
The impact of clay sequences is evaluated by considering two different field configurations (see the top panels in Fig. 10b, c). Each multi-layered claystone formation comprises four thin clay layers embedded in permeable sediments. The claystone layers have random horizontal position, integrity, and depth. In case I (Fig. 10b), and methane migrates upward until it is trapped by the claystone and spreads laterally in the direction of groundwater flow. The plume then flows upward from the left edge of the upper parts of the clay sequence because of the lack of an effective impermeable layer. Methane reaches the aquifer base after 4.8 years with a flow rate of 4.41 m 3 /day. The plume evolution in case II (Fig. 10c) is somewhat similar to case I. Methane migrates predominantly in the direction of groundwater flow, after being deviated by the claystone sequence. The methane plume arrives at the aquifer base after 3.7 years with a flow rate of 1.01 m 3 /day.
As illustrated in Fig. 11, continuous thin claystone layers constitute major safeguards against leakage of methane (free and dissolved phase) to shallow aquifers. Time to breakthrough and methane flow rate to the aquifer depend on the integrity of the claystone and the positions of discontinuities with respect to the leaky gas well. The sequence of claystone layers strongly controls the temporal and spatial behavior of the plume in the subsurface and cumulative volume of methane reaching the shallow groundwater.

Comparison with field studies
The occurrence of thermogenic methane in groundwater may indicate that the gas has migrated upward from a hydrocarbon reservoir into the shallow subsurface. The presence of connective permeable pathways, e.g., faults, interconnected bedding planes, and fractures, in the vicinity of the leaky gas well could result in an early manifestation of methane in groundwater wells (Gorody 2012;Darrah et al. 2014;Cahill et al. 2017). For the leaked methane encountering low-permeability layers, methane preferentially flows along higher permeability sediments, either up-dip (Fig. 9d) or in the direction of groundwater flow (Fig. 10b, c)  ). For such conditions, methane can reach distances of hundreds of meters from the leaky well, delaying the breakthrough to the shallow aquifer system. Elevated levels of methane concentrations have been observed in groundwater wells within approximately 1 km distance to hydrocarbon wells (Osborn et al. 2011;Jackson et al. 2013;Darrah et al. 2014;Brantley et al. 2014;Heilweil et al. 2015). Our results show that methane can be manifested even at distances of larger than 1 km from the leaky well because of the flow deviation along low-permeability rocks.
The temporal and spatial behavior of the methane plume in the subsurface highly depends on the geological and hydrogeological characteristics of the formations overlying the gas reservoir Rice et al. 2018b).
The results indicate that methane can migrate to large distances from the source of leakage and, thus, the concentration of the methane in groundwater wells (if at all) can vary considerably with both time and space. Some observational is being analyzed. The top panels show the geometry of the claystone layers at shallow depths, in black color studies (Molofsky et al. 2013;Li and Carlson 2014;Siegel et al. 2015;McMahon et al. 2017;Nicot et al. 2017;Botner et al. 2018) found no correlation between thermogenic methane in water samples from groundwater wells and just the distance to hydrocarbon wells corroborating our findings.

Conclusions
We built two-dimensional generic numerical models to assess hydrogeological conditions that could possibly lead to groundwater contamination from the upward migration of methane from leaky natural gas wells. We applied miscible and immiscible multiphase flow models to evaluate the influence of methane solubility on the transport and examine the behavior of methane plume by means of each flow model. Furthermore, we conducted sensitivity analyses to investigate the relative importance of (1) hydrodynamic parameters including entry pressure and pore size distribution index in the Brooks-Corey equation, and residual water saturation, (2) methane leakage rate from the wellbore, (3) tilted formations, and (4) integrity and sequences of claystone layers. The use of relatively simple conceptual models allowed us to provide detailed evaluations on the influence of each parameter studied. Our simulations provide the following insights: • This study demonstrates that a generic FEP database can be employed for the development of conceptual models representing a potential scenario that may lead to the contamination of shallow groundwater by hydraulic fracturing. The FEP analysis can be used to understand important mechanisms and identify key risks to groundwater resources ). • Miscible and immiscible flow models result in a slightly different temporal and spatial behavior of the plume in the overburden due to the relatively low solubility of methane in the aqueous phase. Methane solubility has a minimal impact on the transport, especially at shallow depths. • Gas leakage rate and leakage period are among the most important factors controlling the magnitude of methane migration to shallow depths. The persistent long-term leakage leads to the accumulation of large amounts of methane in the overburden, which could potentially be transported into the shallow aquifer system. A pulse of high methane flux is required to generate methane contamination in shallow groundwater for short times; for instance, methane arrival times at the measuring location at a horizontal distance of 2000 m from the leakage source were 11.3 and 5.7 years for the continuous inflow rate of 120 m 3 /day and the base-case model, respectively. Results indicate that long-term methane leakage poses higher risks to groundwater in comparison to the rapid migration along fractures (Moortgat et al. 2018) and faults (Reagan et al. 2015) or during drilling operations (Zhang and Soeder 2016). • The buoyancy component of methane migration can be significant in formations/layers with a sufficient dip angle ( ≥ 3 • ) in the horizontal and vertical direction. Under such conditions, the presence of low-permeability formations further contributes to the horizontal spreading of the plume and methane can be observed in monitoring wells kilometers away from the leaky well. For example, instances of methane occurrence in groundwater reported in some of observational studies (Molofsky et al. 2013;Li and Carlson 2014;Siegel et al. 2015;McMahon et al. 2017) were not spatially correlated with locations of oil and gas wells, which could further confirm the migration of methane to different distances and directions from the gas well. Our simulations corroborate the results of previous studies that methane observation in groundwater is not concentrically aligned around the leakage source (Jackson et al. 2013;Steelman et al. 2017;Cahill et al. 2017;Rice et al. 2018a;Moortgat et al. 2018;Botner et al. 2018). • The low-permeability layer, i.e., claystone, is assumed to be horizontal. In reality, such layers may contain domeshaped and wave-like patterns, similar to surface roughness at the pore-scale level. These structures may retain stratigraphically large portions of the leaked methane. Therefore, the overall contamination may be limited where the lateral extent of the low-permeability layer is larger than a few kilometers. Furthermore, mapping of such features is difficult and future models can be directed at investigating such effects. • In all modeled scenarios with a continuous claystone at shallow depth, the claystone constitutes an effective flow barrier. Simulations demonstrate that methane can be trapped at the base of a low-permeability unit and spread laterally until it reaches a discontinuity or a permeable vertically oriented pathway. Time to breakthrough and flow rates of methane to monitoring wells can vary strongly depending on the integrity, depth, and distribution of low-permeability sediments with respect to the gas well. For example, in the scenario with a discontinuity in the center of the claystone (Fig. 10a), methane is observed in groundwater within 1 km radius around the well, while in scenarios with multiple claystone layers (Fig. 10b, c), methane is observed at a distance of more than 2 km from the gas well. Our simulations indicate that the complex shape of the methane plume in the subsurface, the arrival time to groundwater (in case of occurrence), and distances from a leaky gas well vary significantly based on hydrogeological characteristics of formations intercalated between aquifer and gas reservoir (Osborn et al. 2011;Jackson et al. 2013;Darrah et al. 2014;Steelman et al. 2017;Cahill et al. 2017;Botner et al. 2018).