Fate and Transport of Shale-derived, Biogenic Methane

Natural gas extraction from unconventional shale gas reservoirs is the subject of considerable public debate, with a key concern being the impact of leaking fugitive natural gases on shallow potable groundwater resources. Baseline data regarding the distribution, fate, and transport of these gases and their isotopes through natural formations prior to development are lacking. Here, we define the migration and fate of CH4 and δ13C-CH4 from an early-generation bacterial gas play in the Cretaceous of the Williston Basin, Canada to the water table. Our results show the CH4 is generated at depth and diffuses as a conservative species through the overlying shale. We also show that the diffusive fractionation of δ13C-CH4 (following glaciation) can complicate fugitive gas interpretations. The sensitivity of the δ13C-CH4 profile to glacial timing suggests it may be a valuable tracer for characterizing the timing of geologic changes that control transport of CH4 (and other solutes) and distinguishing between CH4 that rapidly migrates upward through a well annulus or other conduit and CH4 that diffuses upwards naturally. Results of this study were used to provide recommendations for designing baseline investigations.

grey, hard, high plasticity silt and clay 14 . The Mannville Gr (~115 m thick; 100.5-145 Ma BP), located at the base of the Cretaceous, is a regionally extensive, sand-dominated aquifer and defines the base of our investigation.
We measured similar CH 4 concentration-depth profiles below the water table at Sites 2, 5, and 6 (Figs 2-4). The upper portion of the CH 4 profiles are concave downwards and increasing with depth from low (typically < 1.0 mg L −1 ) above the till-shale contact at Sites 2 and 6 or slightly below the anoxic-oxic till contact (12 m (Fig. 4). Below the Base Fish Scales Fm, the concentrations decrease through the Viking and Joli Fou Fms to 236 mg L −1 at the top of the Mannville Gr (881 m BTS). Although the datasets are limited in areal extent, the consistency between the CH 4 depth profiles from all sites and the fact we observe elevated CH 4 concentrations in the Second White Specks Fm at a site located 250 km NE of the study area 15 suggest our data represent general CH 4 conditions in the Cretaceous shales across the WB.
Concentrations of CH 4 through the Quaternary and shale deposits are below the solubility limit ( Fig. S1) as estimated from in vitro gas concentrations corrected for in situ temperatures and pressures. From the Second White Specks Fm to the Base Fish Scales Fm, however, CH 4 concentrations are equal to or exceed the estimated solubility limit suggesting the presence of free CH 4 gas in these formations. This determination is supported by the presence of seismic reflectors believed to be associated with free gas in the same formations at other sites in the WB (pers. comm. D Gendzwill; pers. comm. J. Szmigielski). The Second White Specks Fm is a known early-generation biogenic gas reservoir in Cretaceous sediments of Western Canada 8 .
The dissolved chloride (Cl − ) concentration-depth profiles at Sites 2 and 5 (Figs 2 and 3) positively correlate with the CH 4 data; cross-plotted CH 4 vs. Cl − concentrations yield strong linear relationships (y = 0.04x + 3.4, R 2 = 0.90, n = 15 and y = 0.06x + −28.5, R 2 = 0.89, n = 35, respectively; Fig. S2. Given that Cl − is a conservative species, these strong linear relationships suggest, on the whole, the same transport mechanism is operative for both Cl − and CH 4 with no indication of losses of CH 4 from biogeochemical activity (i.e., no production or removal) in the CH 4 concentration-depth profiles in the till and shale. The lack of sorption-desorption reactions for CH 4 is supported by laboratory experiments and consistent with low concentrations of total organic carbon in the Pierre Shale (~0.8% wt) 16 .
The strong linear correlations over the entire data sets, however, do not apply at the upper boundaries for Cl − and CH 4 profiles at all sites (Fig. S2). We compared the upper boundaries for the Cl − and CH 4 profiles at Sites 2 and 5 to identify the controls on CH 4 depletion at these locations (by either a biological or physical mass loss process). The upper boundaries for both the increasing Cl − and CH 4 depth profiles at Site 2 are at the till-shale interface (about 31 m BG). That both the Cl − (240 mg/L) and CH 4 (1.3 mg/L) concentrations were low and constant at and above this depth suggests the dominant control on both solutes at Site 2 is flushing/dilution by lateral water migration, likely through more permeable sediments in the glacial till 16 . In contrast to Site 2, the upper boundaries for Cl − (450 mg/L) and CH 4 (3.4 mg/L) at Site 5 are located at different depths in the till. The upper boundary for the Cl − profile is located at the base of the oxidized, fractured till zone (about 4.5 m BG) whereas the upper boundary for the CH 4 profile is in the underlying anoxic till (about 12 m BG). Unlike at Site 2, the extension of both Cl − and CH 4 profiles well into the till at Site 5 is attributed to the lack of permeable zones in the till above the shale (i.e., no permeable layers were observed in continuous core samples from the till). The upward migration of solutes from the shale into the overlying till at Site 5 is observed in other tills and also attributed to  dots, and mud gas data from Sites 5 and 6 as solid green squares and solid blue triangles, respectively. a lack of permeable zones 17,18 . The presence of the upper boundary for the Cl − profile in the oxidized till at Site 5 is attributed to the presence of fractures (higher K) in this zone that allows solutes to be transported laterally with flowing groundwater 19 . In contrast, the presence of the upper CH 4 boundary in the underlying anoxic till is attributed to CH 4 oxidation associated with the reduction of SO 4 2− generated in the oxidized till by the oxidation of reduced sulphur during the Altithermal Period, a mid-Holocene warming period between 9000 and 5000 years before present [20][21][22] .
Although the shallow data at Site 6 are limited, the trend in CH 4 at the upper boundary is consistent with Site 2. The geologic log for Site 6 also shows the till has a high sand content, indicative of the presence of permeable sand layers. The fact that the upper boundaries for Cl − and CH 4 at Sites 3 and 4 16 , also located in our study area (Fig. 1b), coincide with the presence of a permeable zone at the top of the Pierre Shale Fm suggests that physical control (i.e., lateral flushing) may be more common at the upper boundary for CH 4 than the biogeochemical reactions (i.e., methane oxidation) observed at Site 5.
The overall strong correlation between CH 4 and Cl − profiles at Sites 2 and 5 was used to estimate Cl − concentrations through the Pierre Shale at Site 6. The good visual fits between the estimated Cl − profile and those measured on squeezed samples from Sites 2 and 5 (Figs 2b and 3b) and the estimated Cl − concentrations being consistent with those measured through the lower 275 m of Cretaceous shales at a site ~250 km SE of the study area 23 (Fig. 4b) provide confidence in the calculated Cl − profile for Site 6. As was the case for the CH 4 profiles, the Cl − profiles exhibit well-defined, slight curvilinear (concave down) depth trends to about 200 m BTS after which they increase in a linear manner from 250 m BTS to near the base of the Pierre Shale Fm (600 m BTS) (Fig. 4b).
Based on Cl − concentration profiles with depth through the WB, the source of Cl − in the Cretaceous shales of the WB was determined to be the Prairie Evaporite Fm, located about 1800 m BG in our study area 23 .
The shapes of the CH 4 and Cl − concentration profiles suggest the dominant transport mechanism for these species is diffusion. Diffusive transport in the Pierre Shale Fm and anoxic tills is supported by low measured hydraulic conductivity (K) values of Cretaceous shale (10 −12 -10 −13 m/s) 24 and anoxic tills (10 −10 -10 −11 m/s) 25,26 in the WB. The dominance of diffusive transport in low K argillaceous sediments is consistent with findings from other argillaceous systems in the WB or proximal to it 15,23,27,28 and in other areas of the world [29][30][31][32][33] .
We modeled the CH 4 and Cl − concentration profiles in the Pierre Shale Fm and Quaternary tills at Sites 2 and 5 assuming diffusion to be the dominant solute transport mechanism. Results of modelled scenarios are summarized in Figs 4 and S6 to S11 and the model parameters used are summarized in Table S1. The curvilinear profiles in the upper 200 m of the profiles developed in response to assumed abrupt changes in the upper boundary condition. These abrupt changes are attributed to the effects of glaciation 17,27 . The good fits between the modeling results and measured data for time frames consistent with known ages of shale and glacial deposition confirm that long-term upward diffusive transport over many millions of years can explain the current measured Cl − and CH 4 profiles. Although the times (i.e., onset of Sutherland Gr. glaciation at 1.6 Ma and onset of Saskatoon Gr. glaciation between 0.1 and 0.6 Ma) used to simulate the measured profiles (especially the upper 200 m of Pierre Shale Fm) are consistent with the timing of the glaciations, the results cannot be considered unique because they are based on a limited understanding of the geologic history of the shale and tills. However, as we show below, the estimates of timing can also be constrained by simulating the δ 13 C-CH 4 profiles.
The δ 13 C-CH 4 depth profiles at all sites exhibit a well-defined curvilinear enrichment in 13 C with increasing depth to about 250 m BG (Figs 2-4). Over this depth, the δ 13 C-CH 4 values increase from about −85‰ to about −70‰. Below about 250 m BG depth, the δ 13 C-CH 4 values are constant, ranging from −65.3 to −68.5‰ (depending on sampling method).
The δ 13 C-CH 4 values (−72‰ to −64‰) and lower concentration ratios of CH 4 to the heavier hydrocarbons (C 1 /C 2+ ; >10; Fig. S12) for all samples collected in this study are consistent with early-generation bacterial gas 8 . Late-generation biogenic gases are often thought to be 13 C-depleted relative to these values and have concentration ratios of CH 4 to the heavier hydrocarbons (C 1 /C 2+ ; >1000) [34][35][36] .
Well-defined changes in the δ 13 C-CH 4 values with depth such as those measured across the Pierre Shale are not documented in the literature. A similar isotopic shift to that observed in the upper 250 m of shale was, however, measured on mud gas samples at a site near Brooks, Alberta (600 km NW of our study area) 37 . In that study, the δ 13 C-CH 4 values increased from −80‰ at depths <50 m BG to −62‰ at about 150 m BG and remained constant at this value to 530 m BG. The similarity between the Brooks profile and our profiles, and given the 600-km distance between these sites, suggests the δ 13 C-CH 4 profiles in our study area represent conditions across the WB. Some authors suggest changes in δ 13 C-CH 4 values can be caused by CH 4 transport in the subsurface 16, 38-40 while others discount the effects of δ 13 C-CH 4 migration on δ 13 C values 41 .
The constant δ 13 C-CH 4 values measured below 250 m BG are consistent with near steady-state CH 4 transport through the shale prior to glaciation (as demonstrated in the modelling). An abrupt change in the upper CH 4 boundary condition, such as glaciation, would trigger transient diffusion that would cause diffusive isotopic fractionation. The isotopic shift in the δ 13 C-CH 4 profiles to 250 m BG is consistent with experimental results that show diffusing CH 4 has a significantly 13 C-depleted δ 13 C value than its source 38 . Because the CH 4 and Cl − profiles in our study can be explained via diffusive transport, we assessed whether diffusive transport could produce the measured δ 13 C-CH 4 profiles by conducting diffusive transport modeling of the δ 13 C-CH 4 depth profiles through the shale at Sites 2 and 5. The simulation results using a consistent set of transport conditions to those used in the CH 4 simulations yield reasonable fits to the measured data (Fig. 4) and suggest that isotopic fractionation as a result of differences in the diffusion of 12 C and 13 C in the CH 4 can explain the measured δ 13 C-CH 4 profile. As was the case for the CH 4 and Cl − simulations, the δ 13 C-CH 4 simulations are non-unique.
Based on high-resolution profiling, we gained an understanding of the origin, fate, and transport of CH 4 in the WB. Specifically, the study provided insights into gas migration processes, CH 4 source area delineation and methods to define upper boundary controls. Findings show the source of the CH 4 and δ 13 C-CH 4 in the shallow groundwater is an early-generation bacterial gas play located in the Lower Colorado sand and Second White Specks Fms. This CH 4 migrates upward from its source area through the overlying shale and into the shallow subsurface via non-reactive diffusion over millions of years. Modelling the migration of the δ 13 C-CH 4 profiles can be used to improve the estimates of timing obtained from modelling of the CH 4 and Cl − profiles. Our observations relating Cl − and CH 4 concentrations in the near surface environment can be used to define the controls on CH 4 (i.e., flushing vs. microbiological). The Cl − and CH 4 concentrations data could also provide insight with respect to distinguishing recent CH 4 production from gas migration in near surface aquifers. For example, if the concentrations of both CH 4 and Cl − in aquifer waters are elevated, the CH 4 may be sourced deeper in the basin. However, if Cl − concentrations are low relative to CH 4 concentrations, CH 4 is likely produced in the aquifer. The −15‰ diffusive shift in δ 13 C-CH 4 over the upper 250 m of the WB can generate isotopic values comparable to those commonly associated with more recent methanogenic CH 4 production in shallow potable aquifers (<−80‰) 22,36,42 ; thus complicating our ability to distinguish CH 4 produced in an aquifer from the impacts of CH 4 that fractionated during migration from early-generation bacterial gas plays. This negative shift in carbon isotopes and co-migration of Cl − allows us to distinguish CH 4 contributions from pre-development CH 4 , recent microbial production, and migration from shallow bacterial gas plays ( 13 C-depleted δ 13 C values relative to −65‰) from fugitive gases that rapidly migrate upward through a well annulus or other conduit originating from bacterial gas plays and deeper thermogenic oil-associated production zones in the WB ( 13 C-enriched δ 13 C values relative to −65‰) 11,[42][43][44] . Although the large δ 13 C-CH 4 fractionation resulting from diffusion can complicate identifying bacterial gas sources, our improved understanding of the regional dynamics of gas migration provides a new means to evaluate CH 4 sources in potable groundwater resources in the WB and similar basins in western North America. Data generated in our study can also inform other valuable studies including those applying the model-derived upward flux of CH 4 (6.6 × 10 −2 mg m −1 d −1 ) to estimate long-term gas production from shallow biogenic gas plays.
The challenge to define the environmental impact of anthropogenic CH 4 on shallow aquifers requires differentiation between natural in situ microbial CH 4 from the migration of microbial or thermogenic CH 4 from depth, from thermogenic CH 4 that migrates from an even greater depth due to anthropogenic activities 45 . Although water samples commonly collected from domestic or farm water wells are used to define baseline conditions 46, 47 , Jackson and Heagle 48 state that such receptor monitoring does not yield baseline concentrations of the aquifer or the regional groundwater supply. We suggest that not wells, but high-resolution profiling be used to provide the much-needed information in environmentally sensitive aquifer systems. This profiling only needs to extend deep enough below the aquifer to allow its evolution to and in the aquifer to be characterized. Based on the profile development in the current study, we suggest the profiling terminate 20-50 m below the bottom of the aquifer. We also suggest, based on the profiles generated in the current study, the sample interval could be reduced from 3m in the current study to 5 m. Assuming the base of the aquifer is located about 50 m BG (a common aquifer depth reported) 47 , and the core hole is terminated 50 m below the base of the aquifer, 20 samples would be collected for analysis. This approach can be integrated with existing drilling programs for new assets (production wells, observation wells, groundwater wells, etc.), where the additional geochemical costs are often negligible compared to drilling costs that could be offset by the potential for major reductions in the future costs associated with fugitive gas migration source tracking and potable aquifer impact assessments.

Methods
Drilling, Sampling, and Analyses. Rotary drilling and geologic logging were performed between September 2013 and September 2014 to 150, 200, 1502, and 1432 m below ground surface (BG) at Sites 2, 5, 6, and 7, respectively (Sites 6 and 7 were located ~800 m apart) (Fig. 1). Continuous core sampling was conducted at all four sites and mud gas logging at Sites 6 and 7 only.
Core samples (~100-150 mm long) were collected from Sites 2 and 5 at 3 m intervals and analyzed for: (1) wet and dry densities (ρ m and ρ d ) and gravimetric water contents (w) (used to calculate total porosity, n T ); (2) dissolved anion chemistry (porewater was squeezed from core samples); and (3) CH 4 , ethane (C 2 H 6 ), and propane (C 3 H 8 ) concentrations (ppmV). Gas concentrations were measured in the headspace of Isojars ® flushed with inert N 2 or Ar prior to sealing (~440 cm 3 ). Core samples collected in Isojars ® at Sites 2 and 5 were allowed to equilibrate with the headspace at standard temperature (25 °C) and pressure (1 atm) for 30 and 101 d, respectively, prior to analyses. For analyses, about 10 cm 3 gas from the headspace was collected from each Isojar ® with a 60 cm 3 syringe and injected into the septa port of an Agilent 7890 gas chromatograph (GC) equipped with a flame ionization detector (FID) to measure light hydrocarbons and a thermal conductivity detector to measure O 2 and N 2 . Based on analyses conducted on Scotty TM 17 L calibration gas standards (concentrations of 0.0010, 0.01, 0.1, and 10% CH 4 and 0.0010, 0.01, 0.1, and 1.0% C 2 H 6 and C 3 H 8 ) and on replicate core samples, the accuracy of the analytical method was ± 5%. The limit of detection (LOD) and limit of quantification (LOQ) for CH 4 , C 2 H 6 , and C 3 H 8 were determined to be 1.1 and 11 ppm, respectively. Concentrations measured in the headspace of the Isojars ® were converted to dissolved concentrations (mg L −1 ) using a previously outlined method 33,49 . These calculations used the core-derived n T values.
At Site 2, the mud gas collection cylinder, located in the mud tank, was connected to an Agilent 490 gas chromatograph (GC) and an IsoTube ® gas-sampling manifold. A Los Gatos Research Inc. (LGR) methane carbon isotope analyzer (MCIA)-24-EP CH 4 was connected in line with the outlet of the IsoTube ® gas sampler at Sites 5 and 6. The C 1 -C 5 concentrations were measured on the mud gas using the Agilent GC every 90 s with a repeatability of 0.5% at 1 mol% and an LOD of 1 ppmV. The LGR CH 4 analyzer measured CH 4 concentrations (0-4,000 ppmV) and δ 13 C-CH 4 values over a CH 4 concentration range of 500-10,000 ppmV every 50 s. The analytical error for the δ 13 C-CH 4 values was ± 1.0‰, based on LGR's model specifications; however, field conditions may have raised the analytical error to ± 1.0‰ -2‰. Although the LGR CH 4 analyzer required no internal calibration during sampling, a CH 4 calibration standard (CH 4 = 10,000 ppmV, δ 13 C-CH 4 = −42.3‰) was analyzed every 6 h as a data quality assurance measure to verify instrument accuracy and precision. The δ 13 C-CH 4 values reported from the LGR were confirmed by submitting 10 gas samples covering the range in values (−67 to −80‰) measured in the Isojars ® and analyzed on the LGR at Site 5 to the Isotope Science Lab, University of Calgary. These duplicate results were in good agreement (mean difference = 0.15 ± 0.7‰). δ 13 C analyses were reported in permil (‰) relative to Vienna Pee Dee Belemnite (VPDB), with an accuracy and precision better than ± 0.2‰. The CH 4 concentrations from mud gas profiling at Sites 2 and 5 yielded the same depth trend as the core-jar samples but were lower; this was attributed to dilution by drilling fluid during hole advancement followed by dilution with atmospheric gases during sample collection in the mud tank 32 . The strong linear relationship (y = 13.6x + 65.1; R 2 = 0.83, n = 76) between the dissolved CH 4 concentrations (mg L −1 ) from core samples collected at Sites 2 and 5 vs. CH 4 concentrations measured during mud gas logging (ppmV) was used to convert CH 4 concentrations from mud gas to dissolved gas concentrations at Sites 2 and 5. Because the collection cylinder used to capture mud gas at Site 6 was reconfigured (to minimize atmospheric contamination into the mud gas line) from that used at Sites 2 and 5, calculated dissolved CH 4 concentrations at Site 6 were determined by multiplying the equation to convert CH 4 concentrations from mud gas to dissolved gas concentrations at Sites 2 and 5 (y = 13.6x + 65.1) by 0.25 to account for dilution of the mud gas concentrations during hole advancement. The strong linear relationship (y = 0.05x-17.6; R 2 = 0.88, n = 56) between dissolved CH 4 concentrations (mg L −1 ) vs. Cl − concentrations (mg L −1 ) measured from core samples collected at Sites 2 and 5 was used to convert calibrated dissolved CH 4 concentrations (mg L −1 ) collected during mud gas logging at Site 6 to dissolved Cl − concentrations.
IsoTubes ® collected at Site 2 were shipped to the G.G. Hatch Stable Isotope Laboratory at the University of Ottawa, Canada, where they were stored for 480 d prior to analysis of δ 13 C-CH 4 using a Thermo Finnigan Delta Plus Advantage IRMS coupled with a VarioEL III, with an accuracy and precision better than ± 0.2‰. IsoTube ® samples collected from Site 7 were shipped to Isobrine Solutions Inc., where they were analyzed within 60 days on a Thermo Finnigan MAT252 IRMS, with an accuracy and precision better than ± 0.5‰. Real-time mud gas δ 13 C-CH 4 values were measured using the LGR CH 4 and δ 13 C analyzer at Sites 5 and 6. A slight enrichment in δ 13 C-CH 4 values measured on Isotube ® samples from Site 7 compared to values measured during mud gas isotope logging at Site 6 (+2.9‰ ± 1.6‰; n = 25; Fig. 4) is attributed to experimental error and/or minor secondary microbial reactions (e.g. methane oxidation) during sampling and storage of the Isotube ® and core samples 33 . To test the reliability of the tubes analysis 32 tubes collected at Sites 2, 3, and 4 that were initially run at 480 d after sampling were resubmitted to the University of Ottawa for re-analysis of CH 4 and δ 13 C-CH 4 1160 d after collection. Results showed that over the additional ~2.5 a of storage since the initial analysis (480 d), the concentration of CH 4 decreased by an average of 13.5% (mean = 13.5 ± 7.0 ppm; n = 32) and the δ 13 C-CH 4 shifted by an average of 0.87‰ (mean = 0.87 ± 1.44‰, n = 32). These data show that long-term storage in Isotubes have a minimal impact on the CH 4 concentrations and δ 13 C-CH 4 values.
Additional details on the methods and measured values from Sites 2 and 5 are provided elsewhere 16, 33 . Transport Modeling. Simulation of the diffusive transport of the dissolved CH 4 , Cl − and δ 13 C-CH 4 profiles was based on modelling conducted 16 on the shallow (<165 m BTS) profiles at Site 2. The simulations were performed using a commercially available finite element model (CTRAN/W) 50 . In the current model, we incorporated the deeper profiles through the Pierre Shale Fm obtained from Sites 6 and 7 and the shallow datasets from Sites 2 and 5 (not modeled in ref. 16). The objectives of the modeling were to assess whether the measured profiles can be attributed to diffusion-dominated solute transport and, if so, can the CH 4 , Cl − , and δ 13 C-CH 4 profiles from each site be explained based on our current understanding of the glacial history of the region (Table S1). The simulations of CH 4 , Cl − , 12 C-CH 4 , and 13 C-CH 4 transport were based on transient one-dimensional diffusive transport: where n e is the effective porosity, D e is the effective diffusion coefficient (m 2 a −1 ), C is the mass concentration of the solute (g m −3 ), z is distance (m), and t is time (a). D e in equation (1) was defined according Fick's first law as: d e e where J d is the diffusive mass flux rate (g m −2 s −1 ). The diffusive transport of 12 C-CH 4 and 13 C-CH 4 were determined separately and used to calculate the δ 13 C-CH 4 profile based on Eqn. 3 40 : C CH std 13 4 ( ) 13 4 where − R C CH ( ) 13 4 is the isotope ratio of CH 4 , defined as where R std is the standard isotopic ratio for Pee Dee Belemnite (PDB; 0.011237) and M is the molecular mass. The values of D e for Cl − were assumed to be 2.1 × 10 −10 and 2.3 × 10 −10 m 2 s −1 in the till 28 and Pierre shale 23 , respectively. Values of D e for CH 4 were assumed to be 2.4 × 10 −10 m 2 s −1 for both the till and shale 51 . The n T values were determined on individual core samples. Mean values of n T for the till and shale were 33.2 ± 4.4 and 34.3 ± 2.9, respectively, were used in the modelling. The value of n e for CH 4 was assumed to be equal to the total porosity (n T ) 51 . The n e values for Cl − were determined to be n e = n T for the till and n e = 0.7n T for the Pierre Shale Fm 16 . The value of D e for 12 C-CH 4 was assumed to be the same as that for CH 4 transport in the shale (2.4 × 10 −10 m 2 s −1 ). The difference in molecular mass between 12 CH 4 and 13 CH 4 results in the D e for 12 CH 4 exceeding that of 13 CH 4 by a factor of 1.0195 based on measurements of diffusion in air 52 . This ratio was assumed to apply to our simulations for aqueous diffusion. As such, the value of D e for 13 C-CH 4 was 2.35 × 10 −10 m 2 s −1 .
One-dimensional transport modelling was undertaken in two phases. Phase 1 defined the distributions of CH 4 , Cl − , and δ 13 C-CH 4 across the shale prior to glaciation, and Phase 2 assessed the impact of Sutherland and Saskatoon Gr glacial periods on the final (present-day) profiles for CH 4 , Cl − , and δ 13 C-CH 4 at Sites 2 and 5. Details of the methodologies used in each phase of modelling are presented in Supplemental Information.
The relative root mean square error (RRSME) 53 was used to compare the simulated and measured concentration profiles for each scenario and determine the best-fit model results (data not presented).