Impact of groundwater flow on methane gas migration and retention in unconsolidated aquifers.

Methane leaking at depth from hydrocarbon wells poses an environmental and safety hazard. However, determining the occurrence and magnitude of gas migration at ground surface is challenging, as part of the leaking gas is retained during upward migration. We investigated migration through unconsolidated sedimentary aquifers using a two-phase, two-component (water and methane) flow and transport model constructed in DuMux. A sensitivity analysis for migration through a 60 m thick sandy aquifer showed that retention by dissolution can be significant even with low groundwater Darcy velocities of 1 m.yr-1. Retention was negligible in the absence of groundwater flow. Besides groundwater velocity, both hydrogeological (permeability, entry pressure, pore-size distribution, and residual gas saturation) and leakage conditions (depth, magnitude and spatial dimensions) determined model outcomes. Additional simulations with interbedded finer grained sediments resulted in substantial lateral spreading of migrating gas. This delayed upward migration and enhanced retention in overlying sandy units where groundwater velocities are highest. Overall, the results of this study show that for unconsolidated aquifer systems and the most commonly observed leakage rates (0.1-10 m3.d-1), significant amounts of migrating methane can be retained due to dissolution into laterally flowing groundwater. Consequently, resulting atmospheric methane emissions above such leaks may be delayed with decades after the onset of leakage, significantly reduced, or prevented entirely.


Introduction
The reliance on fossil fuel resources for the majority of the world's energy supply has resulted in a large number of onshore hydrocarbon wells, conservatively estimated at over 4 million (Davies et al., 2014). In spite of efforts to maintain the vertically isolating function of geological formations that are penetrated when installing and operating such wells, failure of the wellbore system is a commonly observed problem (Davies et al., 2014) and can lead to leakage of hazardous liquids and gases. Furthermore, research has shown that this risk continues or may develop even after the active life-time of wells and their abandonment (Kang et al., 2014;Townsend-Small et al., 2016). Particularly, upward leakage of methane through anthropogenically opened, unintended connections between hydrocarbon reservoirs and the shallow subsurface has become a growing concern worldwide as it may contribute to greenhouse gas emissions (Kang et al., 2014), deteriorate water quality (Vengosh et al., 2014), and form an explosion hazard (Chilingar and Endres, 2005). On top of that, leaky wells could serve as pathways for the migration of other fluids when the downhole conditions are actively changed, for example when hydraulic fracturing is carried out (Brownlow et al., 2016) or CO 2 is stored (Gasda et al., 2004) in nearby wells. Thus, their presence may also hamper the safe implementation of future uses of the subsurface.
To be able to quantify and mitigate these risks, the accurate detection of gas leakage is vital. This is typically achieved by measurements of either sustained casing pressure (SCP) or surface casing vent flow (SCVF) at the wellhead (King and King, 2013). However, these measurements may only reflect a part of the total leakage flux, as a significant fraction of leaking gas can escape the wellbore system entirely and enter the surrounding geology (Forde et al., 2019), which is also referred to as gas migration. Lackey and Rajaram (2018) identified three main mechanisms that may lead to gas migration: (1) gas circumvention, when gas migrates through a cemented outer annulus and a section of low quality cement is overlain by a section of higher quality cement, causing the gas to migrate outwards, (2) groundwater crossflow, when gas leaks through an uncemented outer annulus and is transported into the surrounding formation by lateral groundwater flow, and (3) SCP-induced gas migration, when the gas pressure at the bottom of the surface casing exceeds the hydrostatic pressure leading to gas 'overflow'.
As gas migrating outside the wellbore may become trapped, dissolved or degraded (Cahill et al., 2018), surficial or shallow groundwater wells may only detect leakage after long periods of time, or possibly never. Measurements of soil gas migration are typically carried out using surface flux chambers (Erno and Schmitz, 1996). In the vadose zone, oxidation and dispersion of methane has indeed been shown to be capable of masking leaking wells from being detected at the surface entirely (McMahon et al., 2018;Schout et al., 2019). Furthermore, a recent field experiment where methane was released in a shallow aquifer showed that the combined effects of trapping and dissolution of leaking gas in the saturated zone significantly reduced the amount of gas that reached the surface (Cahill et al., 2018). Lastly, several field-based studies have shown that anaerobic methane oxidation coupled to either sulphate reduction (Van Stempvoort et al., 2005;Wolfe and Wilkin, 2017) and/or reduction of iron and manganese oxides (Schout et al., 2018;Woda et al., 2018) can lead to attenuation of a migrating methane plume. With increasing depth and longer migrating pathways leaking gas can become more severely affected by these attenuation and retention processes, and in turn surficial or shallow subsurface measurements become increasingly less reliable tools for detecting gas migration. Indeed, the fate of methane from leaks occurring at depths greater than 10 m was identified as a key knowledge gap in studying gas migration .
Given the complexity and related costs of measuring gas leakage in field or experimental settings, methane migration has also been assessed by means of multiphase numerical flow and transport simulations. Several studies have aimed to determine ranges of possible gas migration flow rates from the reservoir depth to shallow aquifers through high permeability pathways, such as faults or improperly cemented annuli (Kissinger et al., 2013;Nowamooz et al., 2015;Reagan et al., 2015;Schwartz, 2015). Tatomir et al. (2018) simulated migration of methane at great depth (~1500 m) into an inclined, regional aquifer. Breakthrough times of gaseous methane at various offset distances from the leak origin were calculated, that could for example correspond to conductive faults connecting the deep aquifer with an overlying freshwater aquifer. Rice et al. (2018) investigated SCP-induced gas migration through a shale formation present at the bottom end of the surface casing into an overlying shallow aquifer. They showed that for given source zone pressures, the permeability distribution and the parametrization of the capillary pressure saturation and relative permeability functions of the shale formation controlled the flow rate of methane into the overlying aquifer. Depending on these properties, flow rates at the base of the aquifer can be slow, which may allow methane contamination to go undetected.
Other studies assumed leakage to a shallow aquifer system occurred, and focused on the migration and attenuation of methane therein. Roy et al. (2016) coupled their multiphase model to a reactive transport simulator and showed that in confined aquifers, anaerobic methane oxidation coupled to sulphate reduction could attenuate a migrating dissolved methane plume significantly. This attenuating effect was even larger in unconfined aquifers, where aerobic methane oxidation was also possible. Moortgat et al. (2018) simulated gas phase methane migration through an aquifer system characterized by the presence of highly permeable pathways such as fractures and fluvial channels. They showed that rapid lateral gas migration over distances of several kilometers is possible through these features, but that migration is much less rapid in unfractured media. D'Aniello et al. (2019) simulated leakage of methane from a geothermal well into a 2 m thick surficial unconsolidated aquifer. Given the limited thickness and because mass transfer of methane to the aqueous phase was ignored, retention of methane in this aquifer was limited. Klazinga et al. (2019) carried out 2D simulations imitating the field experiments by Cahill et al. (2017) where methane gas was injected up to 10 m depth in a sandy aquifer. Lateral migration of the gas phase was shown to be significant even over such a relatively shallow interval, as a result of anisotropic sediments, and the presence of low permeability layers. Furthermore, wider plumes and larger groundwater flow velocities resulted in larger amounts of methane that were retained in the aquifer.
The effect of methane retention by dissolution into laterally flowing groundwater has not been considered in detail in previous studies of gas migration. In spite of the relatively low solubility of methane, dissolutive retention could be important in the shallow part of unconsolidated groundwater systems, where groundwater velocities are generally higher. Methane migration through unconsolidated aquifers is also less likely to be dominated by quick gas phase flow through preferential flow paths, which reduces the potential for dissolutive retention. Other factors also possibly play a role, such as the increase in methane aqueous solubility with increasing hydrostatic pressure (i.e. depth) and the lateral spreading of the methane plume as a result of anisotropy and low permeable layers. The interaction of these transport processes and their influences on gas migration were studied in a parameter sensitivity analysis based on 3D, two-phase, two-component numerical simulations across a range of realistic conditions. The impacts of horizontally layered unconsolidated aquifer systems on methane migration and retention were illustrated by simulation of a number of additional scenarios, based on the geology encountered at two recently identified leaking wellbore sites in the Netherlands. The overall aim of this study was to determine whether, and if so, to what extent subsurface methane migration through laterally flowing groundwater is impacted by methane dissolution.

Governing equations and constitutive relations
Numerical modelling was carried out using the open source multiphysics simulation package DuMu x (Ackermann et al., 2017;Flemisch et al., 2011). DuMu x is capable of calculating multiphase multicomponent flow and transport at the continuum scale (also referred to as miscible two-phase flow or compositional flow). For this study, two phases α (liquid and gas) and two components k (H 2 O and CH 4 ) were considered. The model accounts for mass transfer between the two phases as well as the compressibility of both phases. Phase velocities are calculated using Darcy's law. Furthermore, binary diffusion is assumed, and the diffusive fluxes in each phase are calculated according to Fick's law. The fully coupled mass balance equation is then formulated as follows: where Φ is the porosity, ρ is the density, X is the mass fraction, S is the saturation, κ is the intrinsic permeability tensor, k r is the relative permeability, μ is the dynamic viscosity, P is the pressure, g is the gravity vector, D eff is the effective diffusion coefficient and q is a source or sink term. A number of constitutive relationships are needed to close this system of equations. Firstly, the sum of the saturations and that of the mass fractions of the two components in each phase must equal 1: where the subscripts g and w represent the gaseous and liquid phase, respectively. The phase pressures are related by the capillary pressure P c as given by: The formulations by Brooks and Corey (1964) were used for the capillary pressure -saturation relationship and relative permeabilities: where P e is the capillary entry pressure, S e is the effective saturation and λ is the pore size distribution index. S e is defined as: where S lr and S gr are the residual (or 'irreducible') liquid and gaseous saturation, respectively. According to these formulations, the relative permeability of the gas phase is 0 as long as S g < S gr . Hence, the residual saturations in each cell have to reach at least up to the (threshold) residual gas saturation before a neighboring cell can be invaded by the gas phase (Fig. S1). A regularization of the Brooks-Corey P c (S g ) relationship was used for S e > 1 and S e < 0.01. To avoid numerical issues with having an infinite gradient of the P c (S g ) curve, the slope of the curve when approaching these limits was extended up to S l = 1 and S l = 0, respectively (Fig. S1). Both water and methane components can be present in either phase. Concentrations in each phase are assumed to be at thermodynamic equilibrium, meaning that mass transfer between the two phases occurs instantaneously at every time step in each control volume. Methane solubility is determined according to Henry's law and is temperature and pressure dependent, with the Henry's coefficient based on the IAPWS formulations (Fernández-Prini et al., 2003). A freshwater aquifer is assumed and therefore the effects of variable salinities on methane solubility are not considered. The solubility of methane also determines the phase state of the model, with a gas phase appearing once the mole fraction of methane (x l CH4 ) exceeds that of the equilibrium (solubility) mole fraction, and vice versa. Phase appearance and disappearance are accompanied by a primary variable switch: from P l and S g when both phases are present to P l and x l CH4 when only the liquid phase is present. Energy transport is not simulated, and the model is considered to be in thermal equilibrium. For the spatial discretization a vertex centered finite volume method ('box method') is used and a fully implicit, backward Euler method for the temporal discretization (Helmig, 1997). The density of the gaseous phase follows the ideal gas law. Viscosity is calculated as described in Reid et al. (1987). In the range of temperature and pressure relevant for this study, the density and viscosity of methane gas calculated in this manner were shown to be nearly equal to more complex equations of state (Kissinger et al., 2013). The density and viscosity of the liquid phase are both pressure and temperature dependent, according to the IAPWS definitions (Wagner and Pruss, 2002). The influence of the phase composition on the liquid phase properties is not taken into account, since even at the greatest depth considered in this study (480 m of water column) the maximum CH 4 mass concentration is negligible compared to the aqueous phase density (~0.1%). The diffusion coefficient of water in the gaseous phase is determined according to the method in Fuller et al. (1966) and the diffusion coefficient of methane in the aqueous phase according to Reid et al. (1987). Mechanical dispersion was not implemented in the model. In multiphase systems, the dispersion coefficient is saturation dependent, which would result in a much more complex system of equations (Helmig, 1997) and excessive runtimes. Computations were performed on a workstation and parallelized over 16 processors, resulting in runtimes of up to 2 days for the sensitivity analysis and 4 days for the layered scenarios.

Geological context and conceptual model
The range of hydrogeological conditions for which gas migration was studied are representative for unconsolidated sedimentary groundwater systems globally, but were inspired by the conditions that prevail in the subsurface of the Netherlands. The Netherlands is one of the major oil and gas producing nations in Europe with around 2500 onshore oil and gas wells (Schout et al., 2019). A survey of 986 gas wells by the State Supervision of Mines (SodM) revealed some form of well barrier failure at 227 (23%) of these wells (SodM, 2019). Observations of gas bubbles in flooded well cellars showed that leakage of thermogenic gas occurred at at least 13 wells (1.3%). Hydrogeologically, the country is characterized by the presence of shallow sandy aquifer units of Plio-Pleistocene age. These aquifers can be either phreatic or, in the northern and western part of the country, can be confined by overlying Holocene deposits (Fig. S2). A succession of thick, marine clays of Neogene and Paleogene age form the bedrock below the aquifer units in nearly the entire country (de Vries, 2007). Typically, the surface casing of oil and gas wells would at least extend down into these clays, the depths of which can be up to around 500 m below surface in the oil and gas producing areas (Kombrink et al., 2012). The modelled domain is a hypothetical rectangular block of a sandy aquifer within such a unconsolidated sedimentary sequence (Fig. 1a). Gas migration into the base of the model is assumed. Conceptually, it could be caused by a number of possible failure scenarios occurring below the simulated part of the aquifer, including SCP-induced gas migration or gas circumvention.

Domain discretization, boundary and initial conditions
The dimensions of the simulated aquifer were set to 100 × 60 × 60 m (XYZ), sufficiently large so that no boundary effects on gas migration occurred. Owing to symmetry in the y-axis, only half of this domain is simulated. Initial numerical testing showed that model outcomes stabilized when using a grid size smaller than 1x1x1 m. Hence, a grid size of 0.5 × 0.5 × 0.5 m was used, resulting in a total of 1,440,000 cells. The simulation time is 20 years, during which injection of CH 4 is continuous. For the CH 4 flux at the inlet (mol.m −2 .s −1 ) a Neumann type boundary was used over 8 cells from x = 40 to 42 m (2 × 1 m 2 ), or what would be 2 × 2 m 2 area when accounting for the half of the domain that is not modelled (Fig. 1b). The remainder of the bottom face of the model and the two lateral faces parallel to the direction of groundwater flow are treated as no-flow boundaries for both components. The two lateral faces perpendicular to groundwater flow are assigned Dirichlet conditions with a pressure equal to the hydrostatic pressure, albeit with a minor increase in pressure on one side to induce groundwater flow. Lastly, the top face of the model is assigned a noflow condition for H 2 O and an outflow condition for CH 4 . This allows methane to escape the model freely through the top boundary but keeps water flowing strictly horizontally. Conceptually, this corresponds to cases where the gas either escapes to the atmosphere or continues to migrate upwards to overlying layers, depending on the depth of the top of the domain. The model is initially fully water saturated (S g = 0) and no dissolved methane is present anywhere. A thermal gradient of 31.3°C.km −1 is imposed in all simulations, equivalent to the average thermal gradient found in the Netherlands (Verweij et al., 2018), on top of the yearly average ambient air temperature of 10°C. The density and viscosity of water inside the model and at its lateral boundaries are determined during an initialization phase that precedes the actual model run.

Parameter space sensitivity analysis
A total of 17 simulations were run, together encompassing a parameter space representative of the expected conditions for gas migration through unconsolidated sandy aquifers (Table 1). Groundwater flow velocities up to 100 m.yr −1 were considered, as observed for regional groundwater flow in the Netherlands (Bloemendal and Hartog, 2018). The thickness of the model is 60 m. For the base case a depth of 60 m was chosen for the bottom of the model domain, with depths of 240 and 480 m also considered. Preliminary simulations showed that model outcomes are insensitive to the temperature gradient however, as simulations ran with a constant temperature of 10°C yielded virtually equal outcomes. Hence, variations in the geothermal gradient were not considered.
Methane flow rates at the inlet were varied between 0.1 and 10 m 3 CH 4 atm .d −1 (i.e. the volumetric CH 4 flow rate at atmospheric pressure and 10°C) as the majority of reported SCVF rates fall within this range. For example, 68% of wells with SCVF in British Columbia, Canada, had flow rates below 1 m 3 .d −1 and 25% between 1 and 10 m 3 .d (Wisen et al., 2019). Similarly, the vast majority of reported SCVF rates in Alberta, Canada, are also below 10 m 3 .d −1 (Dusseault et al., 2014). It should be noted that very large SCVF rates exceeding 1000 m 3 .d −1 have also been reported (e.g. Nowamooz et al., 2015), but were not considered in this study. Imposed leakage rates were sustained for the full 20 year simulation period. The assumption of a constant leakage rate over long time periods was also made by Rice et al., 2018. Field measurements of leakage from abandoned wells in Pennsylvania that remained virtually constant over a 3 year time period support this assumption (Kang et al., 2016), as do measured SCPs that sustained over measurements periods up to nearly 10 years (Lackey and Rajaram, 2018). While properties of sandy aquifers are well known for single-phase problems, experimental studies where multiphase properties are determined are sparse. To reduce the number of possible scenarios, three sets of experimentally determined values for an air-water system were used for the porosity, permeability, entry pressure, and pore size distribution (Clayton, 1999). As a base case assumption, the properties of a coarse grained fluvial sand (D50 of 0.61 mm) were taken ( Table 1). The effect of migration through finer grained sands was considered by implementing the values determined for a 'hydraulic fill' (D50 of 0.16 mm) and a 'clayey alluvium' (D50 of 0.09 mm). Hereafter, we refer to these as a fine sand and very fine sand, respectively, according to the Wentworth scale of grain size classifications (Wentworth, 1922). Given that the transverse permeability is not known, an anisotropy factor (k h / k v ) of 5 was assumed in all simulations. Variations in anisotropy with a k h /k v of 1 and 10 were simulated in scenarios 14 and 15, respectively (Table 1). A residual gas saturation of 1% was assumed. The sensitivity of model outcomes to the residual gas saturation was investigated in scenarios 12 and 13, with values of 0% and 15%, respectively. The residual wetting phase saturation was set to 10% for all simulations. In scenarios 16 and 17 the inlet boundary size was changed (to 1 × 1 m 2 and 4 × 4 m 2 ), while maintaining the total gas flow rate over the boundary.

Setup of two layered case studies based on real sites
For the sensitivity analysis, a simplified, homogeneous, sandy aquifer was considered. However, in unconsolidated sedimentary basins these aquifers are typically alternated by layers of both coarser and finer grained sediments, ranging from gravel to clay. To investigate the effect of such stratification on gas migration, additional simulations Table 1 Summary of relevant input parameters used for the 17 simulated scenarios in the sensitivity analysis. were carried out based on the hydrogeology observed at two locations where methane leakage was recently shown to occur. At the first location, near the village of Sleen in the east of the Netherlands, gas migration resulted from a blowout that occurred in 1965 (Schout et al., 2018). At the second location, in a village called Monster in the west of the Netherlands, gas migration was detected above a fully decommissioned cut-and-buried gas well (Schout et al., 2019). This leak was more recently closed off by the responsible operator, in accordance with Dutch law. Hydrogeologically, these two locations are quite different as they are situated at opposite ends of the Netherlands (Fig. S2). Hence, they serve as two distinct case studies used to investigate the effects of horizontal layering on gas migration. However, it is important to note here that the goal is not to reproduce exactly the leakage conditions at these sites, as required information about the depth and magnitude of the leaks is not known. Lithological and permeability data were retrieved from the publicly available national hydrogeological model REGIS II (Vernes and Doorn, 2005). REGIS II divides the subsurface into sandy, complex, and clayey layers. Complex layers typically consist of successions of more and less permeable sediments that are taken together as one regional layer. As a result, they often have a large anisotropy with a much lower vertical than horizontal permeability. Where not given, anisotropy was assumed to be 5 for sandy layers and 10 for clay layers. Porosities were assumed to be 30% for sandy layers, 35% for complex layers, and 40% for clay layers. The Brooks-Corey pore size distribution index was assumed to be 1.5 for sandy layers, following the experimentally determined value for the sand used in the base case scenario. Accordingly, a smaller pore size distribution index was chosen for the clayey (0.75) and for complex (0.5) layers, reflecting the larger variation in pore sizes that may be expected for such lithologies. The Leverett J-function (Eq. (7)) was used to scale the entry pressure based on the ratio between porosity and permeability. Separate reference values were used for the complex, sandy and clayey units (Table S1).
The bedrock at the Sleen site can be considered the top of the Breda Formation at 126 m depth. The case study modelled after this site was therefore constructed from this interface up to the surface. At the Monster location, the thick clay unit below the base of the Maassluis formation at 231 m depth was considered as the bedrock. Furthermore, the surficial Holocene deposits (55 m thick) are not included in the model as REGIS II does not provide data on them. The modelled section is therefore 176 m thick. The resulting hydrogeological input data used for both case studies is shown in Fig. 2. For each case two simulations were carried out, one with a groundwater head gradient of 25 cm.km −1 and one with 100 cm.km −1 , on top of the hydrostatic pressure. A flow rate of 10 m 3 CH 4 atm .d − was assumed for all simulations. Otherwise, the boundary and initial conditions were kept equal to those used in the sensitivity analysis. The width of the model was extended to 80 m to avoid gas phase pools below interfaces of entry pressure or permeability reaching the boundaries of the model. Given the extended size of these models and the more complex hydrogeology, grid refinement was slightly reduced to 1 × 1 × 1 m to improve run times.

Results
First, an assessment of the potential for methane retention in the subsurface is briefly presented, and the relative magnitude of forces relevant to methane migration for the unconsolidated aquifers under study. Then, the results of the sensitivity analysis are presented, followed by the results of hydrogeologically layered case studies.

Dimensional analysis of methane retention and migration
Within the range of temperature and pressure considered in this study, that result from a maximum depth of 480 m below water table, methane aqueous solubility and gas phase density increase from 31 to 1112 mg.L −1 and 0.69 to 32.3 kg.m −3 , respectively (Fig. S3). While methane mass density as a gas phase is much greater than its solubility at equal depths, the mass stored in the aqueous phase is roughly equal at gas phase saturations of 4% (Fig. 3). For gas saturations exceeding 4%, retention in the gas phase starts to rapidly exceed aqueous retention. At a gas saturation of 4% and a depth of 60 m below the water table, the cumulative storage in both phases is~0.2 kg.m −3 (Fig. 3). In comparison, the flow rate in our base case scenario is 0.1 m 3 CH 4 atm .d −1 , equal to a methane mass flow rate of 0.07 kg.d −1 or flux of 0.02 kg.m −2 .d −1 (given the inlet area of 4 m 2 ). Therefore, at this depth, it would require 10 days to saturate 1 m 3 of aquifer with methane, assuming a gas saturation of 4%. Although this is just a first approximation, which notably does not take into account the replenishment of available groundwater for methane to dissolve in due to groundwater flow, it confirms that there is indeed potential for subsurface methane retention to significantly affect the monitorability of gas migration originating at depth.
Following the definitions in Kopp (2009) a dimensional analysis of the balance of forces in the system was carried out. Preliminary testing showed that for the gas migration velocities encountered in our sensitivity analysis, viscous forces were insignificant compared to both gravity and capillary forces. Hence, the system is characterized by the dimensionless Bond number: where p cr and l cr are the critical pressure and length, respectively, and g is the gravity constant (9.81 m.s −2 ). The critical pressure is defined as the capillary pressure drop over the saturation front length and is therefore roughly equal to entry pressure. Typically, for advectiondriven flow systems, the critical length is taken to be equal to the length of the saturation front width (Kopp, 2009). When assuming a critical length equal to the discretization length (0.5 m), capillary forces equal gravitational forces for entry pressure values of 5 kPa (Fig. S4). This shows that for the three sands considered in the sensitivity analysis, which have a maximum entry pressure of 3.9 kPa (Table 1), gravitational forces likely dominate. However, the complex and clayey layers in the hydrogeologically layered case studies have a maximum entry pressure of 7.8 kPa (Fig. 2). Depending on the actual width of the saturation front, capillary forces likely exceed gravitational forces for these clayey sediments.

Retention by CH 4 dissolution
In the base case scenario (Table 1) the upward, buoyancy-driven migration of the gas plume from the base of the simulated aquifer to its top (60 m interval) takes 1.65 years ( Fig. 4a and Table 2). Due to the coarse grain size and relatively low methane flux, migration is mostly vertical and the gas phase stays within the column overlying the inlet boundary. This observation is in line with the findings of the dimensional analysis, and confirms that gravitational forces indeed dominate in this case. Maximum gas phase saturations remain low at only 2.4%. As the gas phase migrates, the surrounding water column is saturated throughout, and methane saturated groundwater is advectively transported away from the gas phase plume, allowing more gaseous methane to be dissolved (Fig. 5). Migration time reduces to 1.14 years when the Darcy velocity in the aquifer is reduced to zero (scenario 2) and just 0.35 years if in addition the initial concentration of methane throughout the aquifer is raised from 0 to 95% of the solubility (scenario 3). Therefore, methane dissolution, from the bottom to the top of a 60 m thick sandy aquifer with a groundwater Darcy velocity of only 1 m.yr −1 , causes the migration of a 0.1 m 3 CH 4 atm .d −1 leak to occur G. Schout, et al.

Journal of Contaminant Hydrology xxx (xxxx) xxxx
1.3 years more slowly than it would have without dissolution. The difference in results of the first three scenarios is more pronounced when looking at the percentage of inflowing methane that exits the model through the top boundary (Q out /Q in ). For scenario 3, Q out /Q in reaches 82% in a year and goes to 100% after approximately 5 years (Fig. 4a). Ultimately, this fraction even slightly exceeds 100% as some of the initially dissolved methane also passes through the top boundary (Fig. 4a). For scenario 2, when both the initial methane concentration and groundwater velocity are zero, Q out /Q in increases gradually with time to 91% after 20 years. The remaining 9% is trapped in the aquifer by dissolution and subsequent diffusion away from the gas phase plume. With time, this fraction would steadily increase further as the aquifer saturates with methane and the rate of diffusion slows down. However, for the base case scenario these outcomes are entirely different, as Q out /Q in stabilizes after 4 years when just 10% of the imposed methane flow rate migrates on through the top boundary. This shows that dissolution in combination with advective transport exerts a much stronger control on gas migration than dissolution with diffusive transport, even for groundwater with a Darcy velocity of only 1 m.yr −1 .

Impact of methane flow rate
The relative impact of dissolutive retention is significantly reduced when increasing the methane flow rate through the inlet by a factor of 10, from 0.1 to 1 m 3 CH 4 atm .d −1 (scenario 4). While maximum gas phase saturations increased to 4.6%, migration is still primarily vertical and migration time reduced from 1.65 year in the base case scenario to just 0.16 year (Table 2). Therefore, a smaller fraction of migration methane can be dissolved and retained in the aquifer, given that the groundwater velocity was kept equal, and Q out /Q in stabilizes after just 2 years at 85% compared to 4 years and 10% in the base case scenario (Fig. 4a).

Impact of sand properties
Due to the smaller permeability, larger entry pressure and smaller pore size distributions of the finer grained sands simulated in scenarios 5 and 6 (Table 1), gas phase propagation patterns are changed considerably with respect to the base case (Fig. 5). Particularly for the case of a very fine sand (scenario 6), with a 10 times smaller pore size   3. Potential mass density of methane as a gas phase (green lines) or in the aqueous phase (blue lines) as a function of depth, for three different values of the gas phase saturation (S g ). Temperature is based on a thermal gradient of 31.3°C km −1 and a porosity of 38% was assumed, equal to that of the base case scenario. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) distribution index (0.15 versus 1.5 in the base case), the gaseous plume becomes more rounded (Fig. 5a). This also leads to some minor upstream migration of the gas phase, and is the result of an increase in capillary force, which acts in all directions equally, relative to the gravitational force, which acts only vertically upwards. In the fine sand case (scenario 5), both vertical and horizontal permeabilities are~7 times lower than in the base case. Due to this increased resistance to flow, the maximum gas phase saturation becomes higher than in the base case: 3.0% versus 2.4%. The gas phase can be seen to migrate slightly along with groundwater flow in both scenarios with finer average grain sizes (Fig. 5b). Ultimately, the change in shape and decreased speed of upward propagation allows more methane to be dissolved in both cases, given that Darcy velocities were kept equal. This causes the gas phase propagation to stagnate at 6.5 m and 18.5 m below the top of the model domain for the fine sand and the very fine sand cases, respectively (Table 2). Hence, Q out /Q in remains 0. The distance with which methane has been transported at the base of the model is 28 m in the base case scenario, 24 m in scenario 5 and 34 m in scenario 6, calculated from the inlet at x = 42 m to where the concentration equals half that of the solubility at 60 m depth (204 mg.L −1 ). The difference between these distances is caused by the differences in porosity, resulting in varying effective groundwater velocities. It should be noted that due to diffusion, the distance where methane concentrations are sufficiently large to be readily detected in groundwater samples (~0.1 mg.L −1 ) is almost twice as large (Fig. 5b).

Impact of groundwater flow velocity
Gas migration resulting from a 1 m 3 CH 4 atm .d −1 flow rate is not impacted greatly when raising groundwater velocity from 1 to 10 m.yr −1 (scenarios 4 and 7). Migration time increases from 0.16 to 0.25 year and Q out /Q in is reduced from 85% to 76% (Table 2). However, when groundwater velocity is further increased to 100 m.yr −1 (scenario 8), the gas plume movement is severely impacted and migration stagnates at 53.5 m below the top of the aquifer, having moved up just 6.5 m (Fig. 4b). The impact of a 100 m.yr −1 groundwater flow velocity was also assessed for an even higher methane flow rate of 10 m 3 CH 4 atm .d −1 (scenario 10). The higher flow rate results in larger maximum gas saturations (6.2%) and much more rapid gas phase  Table 1, with simulated variations with respect to the base case indicated in figure legends.
migration, as the top of the aquifer was reached in 0.06 years. However, Q out /Q in stabilized at just 2%. This shows that even large gaseous flow rates can be severely impacted by dissolutive retention, in spite of a very rapid upward migration of the gas phase.

Impact of aquifer depth
The depth of the aquifer was increased 4 and 8 times in scenarios 11 and 12, respectively, resulting in an increased depth of the base of the aquifer of 240 m and 480 m depth (aquifer thickness remains 60 m).
Simulations were run with a methane flow rate of 1 m 3 CH 4 atm .d −1 , and are thus compared to that of scenario 4 ( Table 1). The average methane solubility over the resulting depth intervals is 4.8 and 8.6 times higher than in the base case scenario, and the average gas phase density increases 5.3 and 10.9 times (Fig. S3). Both aqueous and gaseous   Not  (Table 2). By extension, Q out /Q in stabilized later and decreased from 85% to 45% at 240 m depth and to 16% at 480 m depth (Fig. 4).

Impact of the residual gas saturation
In scenario 12 the residual gas saturation was lowered to 0%. Migration time to the top of the aquifer decreased to 1.09 years, compared to 1.65 years in the base case scenario. Q out /Q in increased from 10 to 16%. This shows that even very small variations in residual gas saturation impact gas retention and migration. In scenario 13 the residual gas saturation was raised to 15%. As the residual gas saturation acts as a sort of threshold value, below which flow cannot start, this causes gas phase saturations to become much larger than in the base case scenario, and thus more gas is stored in the porous medium. Migration of the gas phase plume slows down considerably and had nearly stagnated after 20 years of simulation at a depth of 4 m below the top of the aquifer (Table 2).

Impact of anisotropy
Variations in anisotropy were modelled by raising the vertical permeability to equal the horizontal permeability in scenario 14 (k h / k v = 1) and lowering it to 1/10th of the horizontal permeability in scenario 15 (k h /k v = 10, Table 1). The expected higher and lower gas migration speed results in the gas reaching the top of the aquifer in only 0.80 years for the lower anisotropy case and in 3.19 years in the higher anisotropy case, compared to 1.65 in the base case scenario (Table 2). In turn, the change in gas migration velocity affects the amount of methane that is dissolved and Q out /Q in after 20 years of leakage changes to 24% and 2% for scenarios 14 and 15, respectively.

Impact of the area over which the influx occurs
The results of earlier simulations and the dimensional analysis demonstrate that for the flow rates considered in this study and parameter values typical of sandy aquifers, buoyancy driven, vertical gas migration dominates. As a consequence, the degree of dissolutive retention is not only related to the groundwater velocity, flow rate, and sediment properties, but also depends on the area of the inlet. This determines the leakage flux (L 3 .L −2 .d −1 ) and hence the amount of water that is available for methane to dissolve in. To investigate the importance of this effect, the inlet boundary size was reduced from 2 × 2 m 2 to 1 × 1 m 2 in scenario 16, while keeping Q in equal. Upward migration velocity increased considerably and gas reached to the top of the aquifer in 1.12 years while Q out /Q in after 20 years increased to 20%, compared 1.65 years and 10% in the base case. When increasing the inlet area to 4 × 4 m 2 (scenario 17), migration velocity was significantly reduced such that stagnation occurred at 17 m below the top of the aquifer ( Table 2).

Case study 1: Sleen site
The modelling domain of the first case study consists of a simplified hydrogeology with 7 main layers, cumulatively 126 m thick (Fig. 2). The lowest layer is characterized by a high anisotropy factor of 125 (k h / k v ) and low vertical permeability of 4.2•10 −14 m 2 . This complex layer is overlain by 5 sandy sections with slightly varying permeabilities. Notably, the third layer from the bottom is a coarse sand and has the highest horizontal permeability of 1.1•10 −10 m 2 . The sandy units are only interbedded by a 1 m thick complex layer with a low vertical permeability from 1 to 2 m depth. Two scenarios were carried out using this parametrization, one with a constant head gradient of 0.25 m.km −1 and the other with 1 m.km −1 . Given the constant head gradient, groundwater velocities in the coarse sandy layer are highest, with Darcy velocities of 8.6 and 34.2 m.yr −1 in the first and second scenario, respectively (Fig. 2).
The low k v in the highly anisotropic bottom layer causes significant lateral spreading of the methane plume to occur, with gas phase saturations greater than 6% close to the inlet. In the 0.25 m.km −1 head gradient scenario, this spreading occurs more or less symmetrically in each direction (Fig. 6). In the 1.00 m.km −1 scenario the plume is tilted more strongly in the direction of groundwater flow. Gas migration through this 49 m thick layer takes roughly 3.2 years in both scenarios. In comparison, migration through the remaining 77 m thick modelling domain occurs in only 1.9 years in the first scenario (Table S2). Dissolution of gaseous methane in the horizontally flowing groundwater is highly significant, as evidenced by the stagnation of the gas phase plume in the high permeability coarse sand layer in the 1.00 m.km −1 scenario (Fig. 6). The high groundwater velocities here, combined with the lateral spreading of the methane plume that occurs in the underlying layers, causes the 10 m 3 CH 4 atm .d −1 flow rate to be completely dissolved once the plume reaches a depth of around 52 m. In the 0.25 m.km −1 scenario, Q out /Q in stabilizes after 20 years at 41% (Table  S2). The radius of the gaseous plume at the top of the simulated section (which in this case represents the groundwater table) reaches a maximum value of around 20 m, and is slightly elongated in the direction of groundwater flow.

Case study 2: Monster site
Similar to the first case study, the lowest layer of the second case study consists of a complex unit with a low vertical permeability and high anisotropy. However, the overlying stratigraphy is notably different than in the first case study and is characterized by the presence of 5 low permeability, high entry pressure clay layers that alternate the  (Fig. 2), with darker shaded layers representing the clayey and complex layers that have a higher entry pressure and a lower permeability than the sandy layers. sandy aquifers (Fig. 2). Particularly the first and second clay layers from the bottom, at depths of 155-143 m and 119-110 m, have much lower vertical permeabilities (~4•10 −15 m 2 ). Accordingly, they are also assigned higher entry pressures (~7.8 kPa) than the sandy units, and the values considered in the sensitivity analysis. The sandy aquifers are also slightly less permeable than in the first case study, resulting in Darcy velocities of~0.7 m.yr −1 and~2.7 m.yr −1 for the low and high head gradient scenario's, respectively. The modelled section is 176 m thick, and the top of the model is at 55 m depth. Gas migrating through the top boundary would enter the overlying Holocene deposits, which were not considered in our simulations.
As in the first case study, the complex unit at the base of the model causes substantial lateral spreading (Fig. 7). After around 3 years, the migrating gas phase encounters the first clay layer which results in significant gas phase pooling. Invasion of the clay layer only occurs after gaseous saturations exceed 6%. At that saturation, the combined capillary and gravitational forces were large enough to overcome the entry pressure barrier. In the low head gradient scenario, this gaseous pool grows to saturations exceeding 10% and a diameter of roughly 70 m. In the high head gradient scenario, less pooling occurs due to the larger amount of methane that is already dissolved in the underlying units. However, the pool is more extensive in the direction of groundwater flow. Similar behavior is observed at the 2nd clay layer, but entry pressures assigned to the remaining three clay layers were not sufficient to cause significant pooling. In the low head gradient scenario, the gas phase plume reaches the top of the modelled domain after 16.2 years. Even at the end of the 50 year simulation time, Q out /Q in had not fully stabilized reaching a value of 46.3% (Table S2). Due to the relatively slow upward migration, pooling of gas below the clay layers, and low groundwater velocities in the clay layers, the dissolved methane plume takes on a pine tree-like shape (Fig. 7). In the high head gradient scenario, gas phase migration stagnates at 103 m depth, or 48 m below the top of the modelling domain. Notably, the center of the gaseous plume at that depth can be seen to have shifted by roughly 30 m in the direction of groundwater flow. This is caused by the imposed hydraulic gradient, which results in slightly lower water pressures downstream. Given that the entry pressure is defined as P e = P g -P l , the migrating gas will preferentially invade cells those cells and the plume shifts in the direction of groundwater flow.

Upward migration and retention of methane gas
For the simulated homogeneous sandy aquifers and leakage flow rates up to 1 m 3 CH 4 atm .d −1 , gas migration was shown to be primarily buoyancy driven and vertical (Fig. 5), as expected based on an analysis of the relative magnitude of gravitational and capillary forces (Fig. S4), and as observed by other researchers for weakly isotopic sandy porous media (Klazinga et al., 2019). A limited amount of up or down gradient flow of gas only occurred for two simulations with more fine grained sediment properties (Fig. 5), due to their lower vertical permeability and an increase in capillary forces. Gas phase saturations remained low and reached up to around 5%, depending on the imposed flow rate and permeability. At these saturations, the potential mass storage of methane in the aqueous phase is either larger than or equal to the storage in the gaseous phase (Fig. 3), even without considering the additional dissolution capacity with groundwater flow. Gaseous retention has been shown to be significant in simulations of methane migrating from gas reservoirs towards shallow aquifers, when assuming large residual gas saturations of up to 30% (Kissinger et al., 2013). However, there is little consensus on appropriate values to use for simulating gas migration through unconsolidated aquifers. While Klazinga et al. (2019) used a single value of 10%, Rice et al. (2018) chose a negligibly small value, arguing that gas cannot be trapped during processes with strictly increasing gaseous saturations (i.e. drainage), which is the case in our numerical simulations. However, gas migration through real-world, heterogeneous sediments is likely to be a more dynamic process. As a result, intermittent periods of imbibition may still occur, particularly when leakage at the leakage point is not entirely continuous. In spite of this uncertainty, model outcomes showed relatively little sensitivity to a change in residual gas saturation from 1% to 0% (scenario 12, Table 2), because gas migration is more strongly controlled by retention in the aqueous phase at those concentrations. On the contrary, an imposed residual saturation of 15% exerted a strong control on upward gas migration leading to a much more slowly developing gas plume (scenario 13, Table 2).

Retention of methane by dissolution
The sensitivity analysis showed that dissolutive retention of methane in unconsolidated sedimentary aquifers can play a major role in limiting upward gas migration. Even low groundwater flow velocities (1 m.yr −1 ) resulted in significant methane retention (Fig. 4a). On the contrary, retention was negligible in the absence of groundwater flow. This shows that it is primarily driven by advective transport of dissolved methane away from the gas phase plume, rather than diffusive transport. Hence, retention is proportional to the groundwater velocity  (Fig. 2), with darker shaded layers representing the clayey and complex layers that have a higher entry pressure and a lower permeability than the sandy layers. and methane solubility. The effect of retention on methane migration was therefore also much less significant for a simulation where the initial aqueous CH 4 concentration was raised to 95% solubility. In this scenario, migration time through the 60 m thick aquifer was just 0.35 years, compared to 1.65 years in the base case with fully unsaturated groundwater (Fig. 4a). This also highlights that upward methane gas migration can be more rapid in areas with pre-occurring concentrations of dissolved methane. However, natural methane concentrations near solubility are rarely found in shallow groundwater. As an example, the methane solubility at just 60 m depth below the water table (~200 mg.L −1 , Fig. S3) already exceeds the maximum concentration ever observed in Dutch shallow groundwater (~120 mg.L −1 , Cirkel et al., 2015). Hence, dissolutive retention can still be significant even for locations with relatively high methane concentrations.
Besides groundwater velocity and initial methane concentrations, model outcomes indicate a strong control of the sediment properties and the magnitude of the leak, as they determine the velocity and dimensions of the upward migrating gas phase. For a given groundwater velocity, a more dilute or more slowly migrating gas plume allows for a larger fraction of the total flux to dissolve. Sediment properties that limit upward migration include the vertical permeability, the residual gas saturation, and the capillary pressure-saturation and relative permeability relations. Using the Brooks-Corey model, the latter two are defined by the capillary entry pressure and the pore size distribution index. Indeed, for two simulated scenarios with sediment properties associated with more finely grained sediments (D50 0.16 and 0.09 mm, compared to 0.61 mm in the base case), methane migration resulting from a 0.1 m 3 CH 4 atm .d −1 flow rate stagnated entirely due to the complete dissolution of the migrating gas plume at depth (Fig. 5).
Model outcomes also show that the dimension of the inlet over which the leakage is imposed also has a strong influence. Here, we assumed that gas migrated into the model domain over a 2 × 2 m 2 area, and that the origin of the leak was somewhere below the bottom of the model domain, for example caused by gas circumvention or SCP-induced leakage (Fig. 1). This also explains why retention was greatly enhanced in two simulated case studies where the sandy aquifers where alternated by sediment beds of lower permeability (Fig. 2). Besides their greater thickness and depth, such layers cause significant spreading of the gaseous plume to occur, which in turn enhances dissolutive retention significantly (Figs. 6 and 7). Flow rates of 10 m 3 CH 4 atm .d −1 were either completely dissolved, causing migration to stagnate well below the top of the modelled domain, or took several decades to fully develop (Table S2). Significant lateral spreading of migrating methane was also shown to occur for SCP-induced leakage originating in a low permeable shale below a shallow aquifer (Rice et al., 2018). In such conditions, our results show that more permeable overlying layers can become strong barriers to upward migration due to dissolutive retention. Conversely, the impact of dissolutive retention would be much smaller when methane migration is focused through high permeable pathways such as faults or fractures. As methane solubility is strongly depth dependent (Fig. S3), dissolutive retention is likely to be even more significant when gas migration occurs at greater depths than those considered in this study. However, the effect of increasing solubility due to increasing pressures may be counteracted somewhat by increasing salinities, which decrease methane solubility. Also, groundwater is generally more stagnant at depth than in surficial aquifers. As shown in this work, dissolutive retention is limited in the absence of groundwater flow. The complex interaction of these processes, and their effect on the fate of gas migration occurring at depths greater than those considered in this study (up to 480 mbgl), would be a relevant subject for further research.

Gaseous pooling and permeability clay layers to gas migration
While significant gas phase pooling occurred below the low permeable clay layers (k v up to minimum of 4•10 −15 m 2 , P e up to 7.8 kPa) in the 2nd case study, they were ultimately permeable to the migrating gas (Fig. 7). Invasion of a higher entry pressure layer requires a large enough gas accumulation in the underlying layer such that the combined capillary and gravitational forces exceed the entry pressure barrier. For large enough entry pressure values and horizontal barriers, such conditions would not be obtained as the accumulating gas pool spreads out laterally into a thin, pancake like shape. High resolution permeability measurements have shown that thin, low permeability clay lenses can be present within larger clayey formations, with permeabilities that are orders of magnitude lower than the average permeability of the unit (Rogiers et al., 2014). Given their low permeability and presumably higher entry pressure, these lenses would likely act as impermeable barriers to gas phase flow, essentially turning the entire formation into a cap rock. However, the ubiquity and lateral extent of these lenses is poorly understood. Therefore, further research into how heterogeneity of clay deposits impacts gas migration would be required.

Continuum modelling of gas migration
The low observed saturations call into question whether the injected gas actually forms a continuous phase, and hence whether the widelyused Darcy's law approach is applicable for modelling gas migration through unconsolidated aquifers. Air sparging experiments have shown that, for unconsolidated sediments, the injected gas phase migrates as continuous gaseous channels when grain sizes are ≤1 mm. For larger grain sizes the flow pattern gradually changes to discontinuous bubble flow (Brooks et al., 1999). The grain sizes considered in this study fall in the former category, and hence channel flow may be expected to dominate, although the flow type is also dependent on the flux size. Furthermore, whether such migration takes the form of a dendritic network of small channels or a smaller amount of larger gas channels is poorly understood (Brooks et al., 1999). Continuous dewatered channels through which the gas phase dominantly migrates also occurred in the field experiments of gas migration by Cahill et al. (2018). Observed dissolved methane concentrations below theoretical solubility were attributed to this flow pattern, as the available contact area over which mass transfer between phases can occur is smaller, and the collected groundwater samples aggregate both water that has been in contact with a gas channel and water that has not.
As this process is not captured in our simulations, the resulting amount of dissolution may be overestimated. On the other hand, mechanical dispersion was not included in our simulations. Excluding dispersion actually leads to an underestimation of the dissolutive capacity, as dispersion would spread out both the gas phase and the dissolved methane more rapidly, allowing for more contact with the water and quicker transport of dissolved methane away from the gas phase plume. Ultimately, the suitability of continuum models for modelling gas migration should be further analyzed through reproduction of experimental results. A promising experimental method which could be employed for this purpose was recently published by Van De Ven and Mumford (2018), who injected gaseous CO 2 at the bottom of a 2D flow cell and used pH as a proxy to visually track dissolved CO 2 concentrations. Reproducing such experiments with numerical simulation would be a good model validation strategy.

Implications for methane monitoring
From this study a number of important implications for monitoring of gas migration and leakages in the vicinity of oil and gas wells can be derived. Notably, surficial expressions of gas leakage originating at depth in unconsolidated sedimentary basins may take several years to manifest in settings dominated by relatively coarse grained sands, and may take at least up to several decades when migrating gas also encounters interbedded low permeable layers. Surficial detection of leaking oil and gas wellbores may therefore only become possible long after the onset of leakage, and potentially after wells have already been abandoned. This is particularly relevant as in many oil and gas jurisdictions decommissioned wells are cut and buried below the ground surface (Davies et al., 2014;Schout et al., 2019), in which case direct measurements of well integrity failure at the wellhead are no longer possible. Furthermore, for example in the Netherlands there are currently no regulations in place that require operators to monitor the locations of gas wells for extensive periods post-decommissioning.
Relying on surficial flux measurements or near-surface groundwater wells may lead to the underestimation of the actual occurrence and magnitude of gas leakage, and groundwater contamination at depth may remain unnoticed. Hence, measuring dissolved gas molecular and isotopic compositions in groundwater wells in close proximity to potential leakage points is likely a more reliable method. When surface expressions do occur, the dissolution of methane may still be such that measured fluxes at the surface only represent a fraction of the flow rate at the leakage point. Lastly, our simulations show that regional groundwater flow may cause the center of a migrating gas plume to be transported several tens of meters downstream from the well before reaching the surface. Recent field studies have employed search radii of 15 m (Schout et al., 2019) and 20 m (Forde et al., 2019) from the coordinates of gas wells when carrying out surficial flux measurements. While likely suitable in most cases, widening the investigated area in the direction of groundwater flow may need to be considered depending on the conditions.

Conclusions
Oil and gas well failure leading to gas migration in the subsurface poses both an environmental and safety hazard. Effective mitigation of these hazards and the proper assessment of their impact is reliant on the successful detection and quantification of leaks. In this study, upward methane migration through unconsolidated aquifers was analyzed using two-phase, two-component (H 2 O and CH 4 ) numerical simulations in DuMu x . Results show that the retention of migrating methane due to dissolution into laterally flowing groundwater can become significant at groundwater Darcy velocities as low as 1 m.yr −1 . Retention was shown to be dependent on a complex interaction between the lateral groundwater velocity, the depth of the leak, and the velocity and shape of the upwardly migrating gas plume. The latter is in turn a function of the leakage flux and sediment properties (permeability, capillary pressure-saturation and relative permeability relations, and residual gas saturation). Across a range of conditions representative of unconsolidated aquifers and leak sizes up to 1 m 3 CH 4 atm .d −1 , the time it took before the gas plume propagated through to the top of the modelled domain varied from less than 2 months to 5 years. In some scenarios, the total methane leakage rate was dissolved and transported laterally, causing the gas plume to stagnate at depth. Subsurface methane retention was even more pronounced in additional simulations of migration through stratified sequences based on the hydrogeological conditions at two leaking gas wells in the Netherlands, consisting of a number of alternating sedimentary units ranging in grain size from clays to coarse sands. Under such conditions, depending on the imposed groundwater head gradient, flow rates of up to 10 m 3 CH 4 atm .d −1 were either entirely retained in subsurface or took several decades to fully develop. Not only the presence of fine grained and low permeable layers formed barriers to upward gas migration, but also high permeable sands with large groundwater velocities. These allow for a larger dissolutive capacity, particularly when the migrating gas phase has been spread out laterally by underlying low permeability units. Overall, the results of this study show that for the most commonly observed methane leakage rates (0.1-10 m 3 .d −1 ), unconsolidated aquifer systems with lateral groundwater flow can retain significant amounts of migrating methane due to dissolution. Consequently, resulting atmospheric methane emissions above such leaks may be delayed with decades after the onset of leakage, significantly reduced or prevented entirely. Therefore, groundwater contamination and future explosion hazards may go unnoticed.

Declaration of Competing Interest
None