If natural gas production in the United States (US) doubles from 2000 to 2050 as expected, shale gas production will account for more than three-quarters of this natural gas production (U.S.
Energy Information Administration, 2018). This development is largely driven by the combined usage of high-volume hydraulic fracturing and horizontal drilling (U.S. Energy Information Administration, 2018). Other countries and regions (e.g., China, Canada, Middle East) are also ramping up their natural gas production from newly discovered technically recoverable shale reserves (U.S. Energy Information Administration, 2017). In the meantime, putative incidents of deterioration of water quality caused by shale gas development have been investigated throughout the US (Brantley et al., 2014;Darrah et al., 2014;Guo et al., 2014;Nicot and Scanlon, 2012;Vidic et al., 2013;Wen et al., 2019bWen et al., , 2019a;;Woda et al., 2018;Yang et al., 2013), where it is not rare to see that hydrocarbon wells show signs of compromised well integrity (Lackey et al., 2021).
Proving that contamination was caused by shale gas development activities is always difficult, and estimating the frequency of problems is even more problematic. However, for the public to make decisions about granting social license for shale gas development requires some information about the frequency of problems. In this paper we explore the use of a machine learning model on large groundwater datasets from one region to see if such data-driven models can be used in other regions to assess incidents of contamination.
We developed the model by application to data from Pennsylvania (PA), the state in the US with the longest history of conventional oil/gas (OG) drilling and coal mining. PA has become a leading area for shale gas production through use of horizontal drilling and high-volume hydraulic fracturing to extract gas and liquids from the Marcellus and Utica shales. The contaminant that has been identified as the most common problem in PA is methane (CH4), the chief component of natural gas (Brantley et al., 2014). Many studies investigating groundwater quality in shale gas basins have focused on contamination by such anomalous CH4 (Hammond, 2016;Hammond et al., 2020;Li and Carlson, 2014;Molofsky et al., 2013;Nicot et al., 2017c;Osborn et al., 2011;Sherwood et al., 2016;Siegel et al., 2015a;Wen et al., 2017Wen et al., , 2016)). Here 'anomalous CH4' is used to denote dissolved or free-phase CH4 found in groundwater that putatively derived from the effects of recent human activities related to shale gas development. In general, naturally occurring CH4 may have long affected an aquifer whereas anomalous CH4 represents a new source at the impacted site.
The identification of new (anomalous) CH4 is difficult at least partly because naturally occurring CH4 is common in many basins, including in PA groundwater (Baldassare et al., 2014;Siegel et al., 2015b;Wen et al., 2019b). In a limited number of areas, the causes of elevated CH4 in groundwater have been attributed to OG production activities (conventional or unconventional) (Darrah et al., 2014;Grieve et al., 2018;Hammond, 2016;Hammond et al., 2020;Jackson et al., 2013;Osborn et al., 2011;U.S. Environmental Protection Agency, 2015;Warner et al., 2012;Wen et al., 2019b;Woda et al., 2018). Identification of anomalous versus natural CH4 in groundwater often relies on measurement of carbon isotopes in the gas (e.g., Baldassare et al., 2014;Jackson et al., 2013), but this is somewhat expensive and seldomly definitive. Researchers have also sometimes used other hydrogeochemical patterns to identify anomalous CH4. For example, the migration of natural CH4 from deeper formations into shallow groundwater in the Appalachian Basin is often accompanied by salt-containing waters, which leads to slightly elevated salt contents in the shallow aquifer and distinctive chemical signatures (Siegel et al., 2015a(Siegel et al., , 2015b;;Warner et al., 2012;Wen et al., 2019b). Therefore, major and trace elements (e.g., calcium, sodium, and bromide) (Brantley et al., 2014;Cantlay et al., 2020b;Grieve et al., 2018;Lu et al., 2015;McMahon et al., 2017b;Schout et al., 2018;Warner et al., 2012;Wen et al., 2019b;Woda et al., 2018) have been proposed to distinguish varying sources of subsurface fluids. In addition, electron acceptors can be reduced (e.g., SO4 can be reduced to sulfide) by bacteria that use CH4 as an electron donor (McMahon et al., 2017b;Schout et al., 2018;Van Stempvoort et al., 2005;Woda et al., 2018). The concentration of sulfate in waters has sometimes been used to detect anomalous CH4 as it is easier to collect and analyze water samples for sulfate than sulfide. Other tracers were also proposed recently to detect anomalous CH4 in groundwater, e.g., strontium isotopes (Warner et al., 2012;Woda et al., 2018), lithium and boron isotopes (Warner et al., 2014), and noble gas isotopes (Darrah et al., 2015(Darrah et al., , 2014;;Wen et al., 2017Wen et al., , 2016;;Woda et al., 2018). These latter isotopic tools require a high level of expert knowledge and advanced (and expensive) analytical facilities to ensure proper data acquisition and interpretation. Such characteristics so far have prevented these tracers from being applied widely for monitoring of baseline water quality in the area of shale gas production. Among these hydrogeochemical tracers, no single parameter is characteristic of a particular source for anomalous CH4 (Cantlay et al., 2020b). In past studies determining the attribution of anomalous CH4, bivariate plots were often used, in which the dependent and independent variables are generally concentrations of hydrogeochemical parameters (e.g., chloride) and/or concentration ratios (e.g., Ca/Na). No single bivariate plot has been suggested that provides definitive answers about whether anomalous CH4 is present. Instead, the source of elevated CH4 is more likely to be identified through a combination of parallel bivariate plots (Cantlay et al., 2020b;Wen et al., 2019b;Woda et al., 2018).
We wanted to test machine learning approaches to identifying anomalous CH4 in multiple shale gas plays. Logistic Regression (LR), introduced in the late 1980s (McCullagh and Nelder, 1989), is a type of machine learning algorithm used to make dichotomous prediction (e.g., whether the dissolved CH4 concentration is higher or lower than a threshold value). LR has been used in environmental science to predict the presence of certain redox-sensitive contaminants in groundwater, e.g., arsenic (Ayotte et al., 2016), nitrate (Nolan et al., 2002), and the redox state (Tesoriero et al., 2015). LR was also used to predict the probability of CH4 occurrence in aquifers in Alberta (Canada) (Humez et al., 2019). In these studies, a LR model was built based on a few selected categorical or continuous variables. No interactions between variables were explicitly considered. For example, Humez et al. (2019) used sulfate and total dissolved solids as input variables to predict the presence of CH4, but no concentration ratios of these variables were ever considered. Recognizing that the concentrations of chemical constituents in groundwater are prone to differences in the baseline effects of local geology, topography, and hydrology, we suspected that the exclusion of concentration ratio variables in these previously used models largely hampered investigators' ability to generalize the models and findings to other areas. In addition, in the above-mentioned studies, all variables were fed into a single LR model to generate the final dichotomous prediction, and this made it difficult to completely understand how each variable contributed to the prediction.
In this study, we developed an ensemble LR model trained on data (Shale Network, 2015) from the northern Appalachian Basin in PA. Each sub-model focusses on a pair of hydrogeochemical parameters. More complex classifiers like Random Forest or XGBoost (Chen and Guestrin, 2016) were not used because their outputs are harder to explicitly interpret and visualize. We considered predictor features themselves and their interactions. We then explored use of the refined ensemble model on additional held-out hydrogeochemistry data from PA and other data from New York, Texas, and Colorado (Figures 1). To the best of our knowledge, this study is the first to apply an ensemble LR model to predict whether CH4 in groundwater is anomalous or natural. With the improved availability of hydrogeochemistry data, such a prediction model could be potentially beneficial for other areas in the US and even for other countries/regions with existing or planned OG production. Although we do not have definitive evidence for which groundwater samples in each basin are affected by anomalous methane, we validated the model against a dataset of 13 sites with high [CH4] presumed to be caused by gas development activity and 1 site with high [CH4] thought to be unaffected by anthropogenic activities. We also explored the model by testing it against the following hypotheses we derived from the literature: i) the frequency of detection of anomalous methane (# impacted water samples / total water samples tested) in the Appalachian Basin is greater where gas wells target unconventional as compared to conventional reservoirs (Brantley et al., 2014;Ingraffea et al., 2014), and ii) the frequency of detection of anomalous methane in groundwater in the Marcellus/Utica shale play is higher than in western shale plays such as those in Texas and Colorado (e.g., Brantley et al., 2014;Hammond et al., 2020;McMahon et al., 2017a;Sherwood et al., 2016;Woda et al., 2018).
We collated groundwater chemistry analyses from shale plays that are some of the largest in the U.S. and that also are associated with large groundwater datasets: the Appalachian Basin within Pennsylvania (n=23,858;Shale Network, 2015); the Fort Worth, TX-LA-MS Salt, and Western Gulf basins within Texas (n=688; Darvari et al., 2017;Nicot et al., 2017aNicot et al., , 2017bNicot et al., , 2017c)); and the Denver-Julesburg, Raton, San Juan, Paradox, Uinta-Piceance, Greater Green, and North Park basins within Colorado (n=10,258; https://cogcc.state.co.us/). We decided not to include some shale plays such as the Fayetteville shale in Arkansas in our study because 1) that play has shown little evidence for gas migration (McMahon et al., 2017a;Warner et al., 2013); 2) we were unaware of a large dataset of groundwater chemistry for analysis; and 3) this region is not one of the major shale gas production areas. The samples we analyzed were compiled from a variety of data sources (Figures 1; Table 1). A relatively small dataset of groundwater quality in the counties of New York (n=203; Christian et al., 2016) that neighbor Pennsylvania was also tested for comparison because it represents a region of long-standing exploitation of conventional but not unconventional hydrocarbon reservoirs. More detailed descriptions of these datasets and sources are included in the supporting information. Upon publication of this work, all datasets discussed in this study will be released at: https://doi.org/XXXXXXXX. We trained the model to a set of 19,102 groundwater sample analyses that included a small set of groundwater quality data compiled from sites that are presumed to have been impacted by anomalous CH4. We then used an additional subset of analyses from such putatively contaminated samples to validate the model performance. Samples are referred to here as putatively contaminated when they were reported to be contaminated in published journal articles or government reports based on field and laboratory investigations (Llewellyn et al., 2015;U.S. Environmental Protection Agency, 2015;Wen et al., 2019b;Woda et al., 2018). These data were divided into a set for training (n=16) [from Granville Road and Paradise Road (Bradford Co.), and Gregs Run (Lycoming Co.)] and a set for validation (n=13) [from Paradise Road (Bradford Co.) and Sugar Run (Lycoming Co.)]. An additional sample of a well-known brine spring that is naturally CH4-rich from northern PA was also used as part of the validation dataset. Sites were largely determined to be putatively contaminated based on temporal changes in methane concentration, carbon isotope data in collected water samples, and geological investigations.
We strove to develop a machine learning-based ensemble model that could use geochemical analyses to quantify the likelihood that a groundwater sample with dissolved CH4 had been impacted by anomalous rather than natural CH4. A total of 118 geochemical parameters were considered as potential "predictor features" in the model. In other words, these features were tested for their ability to predict if CH4 in a given sample was anomalous or natural. These predictor features divide into three main sub-categories discussed in the next paragraphs. Detailed descriptions of these predictor features can also be found in the supporting information.
Data for many of the concentrations were beneath detection (i.e., censored). For example, many of the CH4 concentrations (31.1%) were censored (i.e., below reporting limits). To maximize the number of measurements for each feature, we chose the fifteen geochemical parameters for which the percent of non-censored measurements was >31.1%. These 15 features are referred to as 'Single Geochemical Features' or 'SGF' in Table 2. These 15 water quality parameters were reported in the training data from PA (Shale Network, 2015) with over 10,000 measurements each (Bicarbonate Alkalinity -Alk, Calcium -Ca, Chloride -Cl, Magnesium -Mg, Sodium -Na, Sulfate or SO4 2--SO4, Total Dissolved Solids -TDS, Barium -Ba, Iron -Fe, Manganese -Mn, pH, Hydrogen Ion -H, Specific Conductance -SC, Temperature -Temp, and Turbidity. Specific conductance of a water solution is a measure of its ability to conduct electricity. Turbidity is a measure of the degree to which the water solution has lost transparency due to the presence of suspended particulates. pH was considered in addition to hydrogen ion concentration because pH, the more commonly reported term, is a calculated parameter [pH = -log (hydrogen ion concentration)] that is not directly comparable to the other concentrations.
The second and third sub-categories of predictor features are the reciprocals of the concentrations of all the SGF except for pH and water temperature (n=13; 'Reciprocal Feature'), and the ratios of any two SGF ('Ratio Features'). Ratio features were calculated for all features except for pH, H + , SC, water temperature, turbidity (n=90). The inclusion of features in the second and third categories was motivated by the observation that these two types of arithmetic combinations of geochemical parameters have been widely used by geoscientists in tracing the source of solutes in groundwater because they can represent stoichiometric coefficients in the governing geochemical reactions (e.g., Bau et al., 2004;Brantley et al., 2014;Cantlay et al., 2020dCantlay et al., , 2020aCantlay et al., , 2020c;;Tisherman and Bain, 2019). Furthermore, we reasoned that the stoichiometries of these geochemical reactions need not be restricted to one hydrocarbon basin alone, and thus might be more predictive across basins. Subtraction and multiplication of SGF were not considered because, unlike ratios, there are no known fundamental reasons why such functions should be predictive.
A machine learning model using only geochemical predictor features allows us to focus on assessing the interplay of geochemical parameters as well as their relative importance in determining whether the elevated [CH4] present in groundwater samples is likely to represent anomalous CH4. Therefore, no non-geochemical features (e.g., land use and bedrock geology) were considered in this study.
Our workflow included sequential phases: ensemble model development (which included training), validation, and application (Figure 2). Taking two parameters from the 118 parameters at a time without repetition yields a total number of 6,903 pair combinations. We started with these 6,903 possible geochemical pairs that were each tested in a sub-model. The process of model development described in the next paragraphs identified the subset of 6,903 predictor pairs that were most successful in determining the likelihood of a groundwater sample being impacted by anomalous CH4. The likelihood was calculated as the ratio of the number of sub-models flagging the CH4 as not being naturally derived (and therefore, putatively, "anomalous") divided by the total number of models considered. The collection of sub-models is termed here the ensemble model.
At the stage of ensemble model development (green in Figure 2), sub-models using the full list of 118 geochemical features paired into 6,903 predictors were trained on 19,102 analyses referred to as training data (Table 1; Figure 2). Each of the 6,903 sub-models (Figure 2) considers only one pair of predictor features. In effect, each sub-model imitates the traditional procedure of using bivariate plots of geochemical variables to investigate the origin of dissolved CH4 in groundwater [see, for example, Figure S10 in Woda et al. (2018)]. The performance of each submodel was evaluated with respect to two tasks. First, it was evaluated for its ability to predict whether the [CH4] is above or below the threshold identified as potentially problematic, 10 mg/L, in the 19,086 water samples (Task 1). Given that these samples were collected by professional consultants working for gas companies before drilling gas wells within a few kilometers of the water wells, the waters could be contaminated or uncontaminated by CH4 from previous oil/gas activity, but we have no outside evidence of contamination. We hypothesized that successful submodels would predict [CH4] < 10 mg/L for samples reported to have [CH4] ≥10 mg/L. We inferred these samples are impacted by anomalous CH4. This approach is implicitly built on the assumption that natural high-[CH4] waters are chemically distinct from anomalous high-[CH4] waters.
However, we acknowledge that our method will not detect anomalous methane in waters with very low values of [CH4]. In Task 2, each sub-model was additionally evaluated to see if it successfully identified the CH4 in 16 of the putatively contaminated training samples as anomalous (Task 2).
The concentration of 10 mg/L was selected because it is a threshold above which immediate action is needed (Eltschlager et al., 2001;Wen et al., 2019b).
During ensemble model development (green in Figure 2), a linear classifier -logistic regression (LR) -is used in the two tasks with 5-fold cross-validation. LR is generally used to solve binary classification problems. With LR, the relationship between predictor features and the prediction outcome can be described by a sigmoidal function as illustrated by Equation (1) in which x represents predictor features and y is the binary model output of 0 ([CH4] <10 mg/L) or 1
In the training, the sub-model allows optimization of w and b to best predict the training data.
Once w and b are optimized, each sub-model can be used in Task 1 to predict which groundwater samples are characterized by [CH4] ≥ 10 mg/L. For this task, Area Under the Curve (AUC) is considered as the performance metric where a higher AUC value indicates higher accuracy of prediction. AUC, in particular, was chosen as the accuracy measure as it holistically evaluates the model performance in predicting both positive (i.e., high methane sample) and negative (i.e., low methane sample) results.
In Task 2 the sub-model was applied to 16 putatively contaminated samples. F1 score was calculated following Equation (2) as the performance metric:
Here recall is defined as the proportion of the 16 samples identified correctly while precision is the proportion of predicted anomalous CH4 samples that are actually contaminated. F1 score was chosen to prioritize maximizing true positives (i.e., known problematic samples correctly identified) as it does not consider true negatives (i.e., known natural samples correctly identified).
Both Tasks 1 and 2 were used to identify the best ensemble model to detect anomalous CH4 in impacted waters. Samples reported to have high [CH4] that were incorrectly classified as low [CH4] samples by sub-models in Task 1 were inferred to have been impacted by anomalous CH4.
We chose sub-models with higher accuracy of prediction in Task 1 to yield a lower number of false positives. To evaluate the performance of each sub-model for both Tasks 1 and 2, we defined a synthesized metric, namely, 0.7 x F1 score + 0.3 x AUC. The selection of coefficients 0.7 and 0.3 is not entirely arbitrary. We assigned the higher weight (0.7) to F1 as that task is more critical.
All sub-models (and the associated pair of features; n=6,903) were ranked by the final synthesized score. We also completed a sensitivity test of the coefficients (Table S1a andS1b) and determined that feature members in the top feature list (e.g., top 20) are mostly unchanged (e.g., SC, Na, Cl/SO4), although the exact rank of features by frequency changes to a limited extent.
In the phases of model validation and application, the prediction results for the best performing sub-models were used to estimate the likelihood of impact of groundwater samples by anomalous CH4 (Figure 2). This likelihood was defined as the percentage of considered sub-models that classified the sample as impacted by anomalous CH4. We calculated four types of likelihood for each water sample using the best n sub-models where we explored n = 100, 200, 500 as well as 1,000. We ultimately chose n = 1,000 to be conservative. The selection of n =1,000 is supported by the observation that (1) an ensemble model with fewer sub-models would be more likely to be biased towards a few geochemical parameters (e.g., the ensemble model for n = 100 was observed to be highly biased towards SC as discussed in the text below) and (2) an ensemble model with a larger number of sub-models (e.g., n = 1,000) was expected to assess water chemistry more holistically, improving model applicability and performance for new datasets with less likelihood of overfitting.
In this study, the ensemble model was trained and validated with hydrogeochemistry data in Pennsylvania (i.e., Marcellus shale) and then was tested in other natural gas production regions in Colorado, Texas, and New York. Differences in bedrock geology, regional hydrogeochemistry, topography, and land use between testing and training/validation datasets might in some cases hamper the applicability of the developed ensemble model. We strove to minimize this issue by considering interactions of ratios of concentrations of single geochemical features as these concentration ratios tend to yield less variability across different regions. As more hydrogeochemical data become available, in particular from presumably contaminated sites, the ensemble model can be further validated and even re-trained with new data.
Two statistical methods -Pearson correlation and Spearman's rank correlation -were adopted to assess the pairwise correlation of the 15 single geochemical features in predicting the target feature, [CH4]. The Pearson correlation is used to evaluate linear relationships between continuous variables that are assumed to obey normal distributions, while the Spearman's rank correlation is used to evaluate monotonic relationships between two continuous variables without the requirement that variables are normally distributed. Therefore, the Spearman's rank correlation identifies more pairs of features that report statistically significant correlations (p < 0.05) compared to the Pearson correlation (Figure 3).
The pairwise correlations led us to note three groups for the fifteen single geochemical features: 1) salinity-related (Alk, Ba, Ca, Cl, Mg, Na, TDS, SC); (2) redox-sensitive (Fe, Mn, and SO4); and (3) other (pH, Temp, Turbidity, and H + ). These groupings are used to categorize the best-performing sub-models in subsequent discussions.
The synthesized scores of all sub-models (n=6,903) are listed and ranked in Table S1a. The best performing sub-model (synthesized score = 0.9579, F1 = 0.9677, AUC = 0.9350) corresponds to the pair, 1/SC and Ca/Ba (Table S1a). From this observation we inferred that salinity is an effective tool to detect anomalous CH4 samples (SC, Ca, and Ba are all salinity-related features).
Previous studies (e.g., Brantley et al., 2014;Cantlay et al., 2020c;Tisherman and Bain, 2019;Wen et al., 2019b;Woda et al., 2020Woda et al., , 2018) ) have also identified other salinity-related parameters (Cl, Ca/Na, Mg/Na, Ca/Mg, and Ba/Cl) that are helpful in detecting anomalous CH4 in groundwater samples or distinguishing chemical signatures of produced waters from shale gas development from other types of contamination. These same parameters, i.e., Cl, Ca/Na, Mg/Na, Ca/Mg, and Ba/Cl, show up in sub-models ranked as high as 8 th , 13 th , 19 th , 82 nd , and 191 st (Table S1a), respectively. Furthermore, the frequency of inclusion of these features in the top 1000 sub-models are 87, 4, 2, 11, and 8, respectively (Table S1b andFigure 4). In addition to the salinity-related parameters mentioned by previous authors, a few new ones were identified as important (i.e., high frequency) in the top 1,000 performance sub-models, e.g., Ca/TDS (49) and Ba/Alk (81) (Table S1b). These salinity-related features can effectively distinguish natural migration of thermogenic CH4 from anomalous CH4 presumably because CH4 often migrates naturally with salt-containing Appalachian Basin brine into shallow groundwater (e.g., Warner et al., 2012;Wen et al., 2019b;Woda et al., 2018). Alkalinity may appear because it increases during sulfate reduction, one of the redox processes (see next paragraph) that sometimes couple to methane oxidation (Woda et al., 2018).
In addition to salinity, redox-sensitive features (e.g., Fe/SO4, SO4/Na) are also frequently found in the top performing sub-models (Table S1, and Figure 4). For example, SO4 and Fe concentrations have previously been used to detect anomalous CH4 in groundwater (Wen et al., 2019b;Woda et al., 2020Woda et al., , 2018) ) and the frequency of SO4-containing and Fe-containing features (e.g., Fe/SO4) appearing in the top 100 sub-models are 17 and 10, respectively. It is worth noting that many SO4-or Fe-containing features also incorporate salinity-related features, e.g., SO4/Na and Fe/Ba. Redox features are expected because the presence of CH4 creates an anoxic environment that in turn promotes reduction of redox-sensitive species such as sulfate (SO4 2-) to sulfide (S 2-) and Fe(III) to Fe(II). The formation of sulfide leads to precipitation of metal sulfide and a decrease in SO4 and Fe concentrations over time (Woda et al., 2020(Woda et al., , 2018)). Low concentrations of all of these redox-sensitive parameters are expected when the source of CH4 is natural and the groundwater has had a relatively long time to reach thermodynamic equilibrium.
In contrast, when anomalous CH4 migrates into shallow, oxygenated groundwater, months can pass before the water reaches chemical equilibrium, allowing the co-existence of transiently high methane, sulfate, and iron in the same water samples that can be detected during the transient (Wen et al., 2019b;Woda et al., 2020Woda et al., , 2018)).
Best performing features are summarized for the top 100, 200, 500, and 1,000 models in Figure 4 and the likelihood of the presence of anomalous CH4 in groundwater is calculated for all of these four top lists in the following sections. It is clear that in the first two lists (i.e., top 100, 200), the frequency of SC-containing features (SC and 1/SC) dominate over the other top features.
So, if the likelihood of groundwater samples being impacted by anomalous CH4 is calculated from an ensemble model considering only the top 100 (or 200) features, the calculated likelihood will be largely biased towards SC. Two more features (Na and Cl) are reported with high frequency in the top 500 list, but the calculated likelihood of groundwater samples being impacted by anomalous CH4 still depends strongly on the salinity features SC, Na, and Cl. To holistically assess groundwater chemistry and to be more conservative, we focus the most extensive discussion mainly about the results based on the list for the top 1,000 (the "top 1,000 model", Figure 4).
After training with the use of putatively contaminated sample data (n = 16), the top 1,000 submodels were validated against another groundwater quality dataset of 14 high-[CH4] groundwater samples from northeastern or central PA. All were reported with [CH4] ≥10 mg/L. Of these, 13 were reported for 4 sites putatively impacted by anomalous CH4 (Llewellyn et al., 2015;Wen et al., 2019b;Woda et al., 2018). The likelihood of contamination was >90% for each of these 13 validation samples (Tables 3 andS2) by the ensemble model using the top 1,000 sub-models. This suggests that the ensemble model has high sensitivity (i.e., high recall rate), i.e., the model is very effective in detecting anomalous CH4 in groundwater samples (true positives). The 14 th sample that we used for validation was groundwater from Salt Spring State Park in northern PA. That spring sample is known to contain naturally occurring CH4 with [CH4] > 10 mg/L ('Salt Spring State Park' in Table S2). The predicted likelihood that this groundwater sample is impacted by anomalous CH4 is 14%, which means that the ensemble model predicts that this sample is likely not contaminated by anomalous CH4 (i.e., a true negative). Thus, the proposed ensemble model is very effective in distinguishing whether high [CH4] in groundwater is naturally occurring (e.g., <20% likelihood) or anomalous (e.g., >90% likelihood).
We now apply this ensemble model to the hold-out Pennsylvania data (data not incorporated in the training dataset) as well as groundwater data from other states (Figure 5) to seek to detect anomalous CH4 in other groundwater samples as a way to explore our two hypotheses and to compare to previous research (Christian et al., 2016;Darvari et al., 2017;Nicot et al., 2017aNicot et al., , 2017bNicot et al., , 2017c;;Sherwood et al., 2016;Wen et al., 2019bWen et al., , 2018)).
A total of 4,772 groundwater samples from the Shale Network database are in the Pennsylvania hold-out dataset. Among these, only 64 show both [CH4] ≥ 10 mg/L as well as measurements for all the single predictor features. These 64 groundwater samples are all located in northeastern PA (Figures 5 and6; Table S3). The top 1,000 model predicts that 29 of these 64 groundwater samples show a >50% likelihood and, of these, 13 show ≥80% likelihood of being impacted by anomalous CH4 (Figures 5 and6; Table S3).
To assess the efficacy of the ensemble model for Pennsylvania, we compare the machine learning results to the results based on a more simplified and streamlined approach by Wen et al. (2019b). Based on six selected geochemical parameters (CH4, Cl, Ca, Na, Fe, and sulfate concentration), Wen et al. (2019b) categorized groundwater samples into five geochemical types.
Their type 4 and 5 waters were considered to the most likely to be impacted by anomalous CH4: their characteristics include Ca/Na mass ratio ≥0.52, Cl ≤30 mg/L, [CH4] ≥10 mg/L, and disequilibrated Fe concentrations (≥0.3 mg/L) and/or dis-equilibrated sulfate concentration (≥ 6 mg/L). Using this simplified workflow from Wen et al. (2019b), we identified a total of 17 samples of types 4 or 5 in the hold-out PA dataset. Ten of these 17 samples were used to train the ensemble model and five were not considered by the ensemble model because they lacked values for all of the 15 geochemical features. The remaining two analyses (sample IDs are "PADEP_Predrill_CHK_827" and "PADEP_Predrill_Bradford-Burlington_T_005-well"; Table S3) were detected accurately (true positives) by the top 1,000 model as showing a high likelihood of impact by anomalous CH4.
In addition to the two samples detected by both the ensemble model in this study as well as by Wen et al. (2019b), the ensemble model also identified 27 other samples with a predicted likelihood of >50% (12 samples ≥80%) of anomalous CH4. The workflow of Wen et al. (2019b) did not detect these 27 samples as impacted by anomalous methane because Cl ≥30 mg/L or sulfate ≤6 mg/L. If the new ensemble model approach is correct, these samples are impacted by anomalous methane and were false negatives in the previous study. Lending credence to the conclusion that these sites were impacted by anomalous methane is the observation that all but 11 of these 27 samples were from four areas that lie near locations with known contamination as reported in the literature (Llewellyn et al., 2015;Wen et al., 2019bWen et al., , 2018) ) (Figures 5 and6; Table S3). The last eleven putatively-false-negatives are not close to any sites previously identified as problematic.
Given that all the other 16 sites show some history of local contamination, we conclude that all these putatively-false-negatives were impacted by recent intrusion of anomalous methane (which might be caused by nearby shale gas drilling) at the time of sampling, and should have received in-depth analysis to determine the source and mechanism of the elevated [CH4].
Although the streamlined workflow from Wen et al. (2019b) is easy to implement and effective in detecting anomalous methane samples, it is limited in the number of geochemical parameters considered and it likely returned at least 27 false negatives for the PA hold-out dataset.
In contrast, our ensemble model evaluates a fuller set of groundwater chemistry data and can identify waters that are contaminated that were missed by the streamlined workflow of Wen et al. (2019b).
Although New York is almost geologically identical to Pennsylvania in the counties considered here, it is not producing shale gas because high-volume hydraulic fracturing with horizontal drilling is banned. Instead, conventional oil and gas drilling has long been important in the state. Brantley et al. (2014) summarized reports from the PA regulator that showed that the reporting rate of problems of anomalous CH4 from conventional oil and gas well development in PA was much smaller than that observed for shale gas development after 2004. We thus hypothesized that we would see a very low rate of anomalous CH4 in the small dataset of groundwater quality available for the region of oil and gas production in New York (Christian et al., 2016). This hypothesis is especially compelling because geologically, northern PA is very similar to NY, but the laws are different, and the use of high-volume high-pressure hydraulic fracturing in NY unconventional reservoirs is not allowed. Therefore, unlike PA samples, all of the New York samples were collected in advance of high-volume hydraulic fracturing but might postdate some of the nearby conventional drilling.
In the New York dataset of 203 groundwater samples, only five are reported with high [CH4] (≥10 mg/L). Four of these five samples are predicted by the top 1,000 ensemble model as low likelihood (i.e., <30%) of containing anomalous methane (Tables 3 andS4; Figures 5 and6).
However, one ("ST35B"), with a [CH4] of 13.8 mg/L, yields a 71% likelihood of impact by anomalous CH4 (Table S4 andFigures 5 and6). Sample ST35B was identified in the ensemble model because it has a very low [Cl] of 4.1 mg/L, but relatively high [sulfate] and [Fe] of 8.99 and 0.42 mg/L, respectively. It is thus also detected as a contaminated sample by the streamlined test of Wen et al. (2019b) and is likely, based on these tests, to contain CH4 that has infiltrated the groundwater recently.
Although this sample appears to have been infiltrated with CH4 recently, other lines of evidence suggest that the migration was not caused by oil or gas development. For example, the distance between site ST35B (sampled in 2013) and the nearest oil or gas well (drilled prior to 1980) is ∼6 km. The large distance and long period of time between gas well drilling and water sample collection render the source of methane in water wells unlikely to have originated from the gas wells. But this water well is drilled into the very gas-rich Canadaway Formation near a large number of faults and lineaments, and may nonetheless have been affected by recent migration of naturally occurring CH4 to the water well. Consistent with this conclusion, Christian et al. (2016) pointed out that values of [CH4] in the NY groundwater data (including, implicitly, ST35B) do not statistically vary with proximity to gas wells or faults. They argued that elevated [CH4] in the NY samples is most likely of natural origin. This conclusion was based on the observation that high [CH4] is often associated with Na-rich waters in valleys in NY that are impacted by natural brines (Christian et al., 2016). However, while ST35B had a high [Na] (30.1 mg/L), it had very low [Cl] (4.1 mg/L). Therefore, the high Na in that sample did not derive from contamination by NaCl-rich brine (and was not flagged as such by the ensemble model). On the other hand, water from ST35B was reported to smell like 'rotten egg' (i.e., H2S) and is locally associated with brownish red stains (pers. comm., Laura Lautz), perhaps pointing to ongoing sulfate reduction coupled with methane oxidation (Wen et al., 2019b;Woda et al., 2020Woda et al., , 2018)). Our prediction that ST35B contains anomalous methane might be a false positive: in particular, it could be water that was recently impacted by natural CH4 rather than CH4 migrated from hydrocarbon development activities. In future predictions with the ensemble model, a filter for proximity to oil/gas wells should be included (e.g., waters sampled within 5 km).
A Fisher exact test was used to test our original hypothesis, i.e., to test if the frequency of detection of CH4-impacted groundwater samples (≥50% likely to be impacted by anomalous methane) differ between PA (i.e., 42/2,983) and NY (i.e., 1/78). Whether we assume that ST35B is an example of anomalous methane caused by oil/gas development or not, the Fisher test results
show no statistically significant difference between PA and NY at a confidence level of 95% (p>0.05). Thus, we must reject the hypothesis in that the ensemble model reveals similar frequencies of putative contamination by anomalous methane in the part of the Marcellus shale play utilizing horizontal drilling + hydraulic fracturing as compared to the part with conventional resource development.
We also compiled groundwater quality data from Texas to test the hypothesis that the frequency of identification of anomalous methane in groundwater wells is higher in PA than in other shale plays/states. A total of 688 groundwater samples are included in the compiled Texas data, which cover three major shale gas plays in Texas: Barnett Shale, Eagle Ford Shale, and Haynesville Shale. Among these samples, [CH4] is larger than 10 mg/L in 34 groundwater samples (Table S5 andFigures 5 and6). Of these 34 samples, none were identified to have a likelihood of anomalous CH4 larger than 50%. This is consistent with the findings from previous investigations of the sources of CH4 in these three shale plays (Darvari et al., 2017;Nicot et al., 2017aNicot et al., , 2017bNicot et al., , 2017c)), i.e., no CH4 in these shallow water wells is associated with recent development of nearby shale gas.
The ensemble model was also applied to a large groundwater quality dataset of 10,258 water samples collected across Colorado mainly within the Denver-Julesburg, Raton, San Juan, and Uinta-Piceance basins (Figure 1). Of these samples, high [CH4] (≥10 mg/L) samples (n=58) with complete measurements for all required single geochemical features exist only in the Denver-Julesburg Basin and Raton Basin (Figure 1). These 58 samples were collected from 29 sites (Table S6), of which only four samples from four sites (i.e., 705739, 755481, 752672, 750143) With these data, we test our second hypothesis, namely that the frequency of detection of anomalous methane in the Marcellus shale play is higher than in Texas + Colorado. Treating the Colorado and Texas measurements as one dataset, the combined frequency of detection of groundwater samples showing >50% likelihood is 4/795, statistically lower than the rate in Pennsylvania (i.e., 42/2,983) at a confidence level of 95% (p=0.04) using the Fisher exact test. This finding supports our second hypothesis that the frequency of detection in the Marcellus shale play is higher than in Texas and Colorado. To explain this result, we point to Hammond et al. (2020) who suggested that most of the anomalous CH4 released from gas wells in the US occurs because primary cementation is not completed along the full lengths of production casings from the target shale to intermediate casings (or to surface casings, if intermediate casings are not used).
The higher frequency in the Marcellus play could therefore derive from differences in casings and cementation or differences in gas contents at intermediate depths for boreholes in that play as compared to other plays.
We presented a machine learning ensemble model that shows that salinity-related and redoxrelated measurements are effective geochemical features that can detect anomalous methane in groundwater. One problem with the application of the model to date is that we cannot exclude the possibility that differences in sample sizes contributed to differences in detection frequencies.
Clearly, larger datasets for all shale gas plays could be used to eliminate this problem. Furthermore, a machine learning model that incorporated additional geochemical measurements such as isotopes would presumably be even more adept at finding evidence of anomalous CH4.
For the regions we studied (PA, TX, CO, NY), the frequency of reported water samples with [CH4] ≥10 mg/L, the level often considered to be dangerous, was 2.1%, 5.8%, 7.4%, and 3.4%. In contrast, the frequency of identification of anomalous methane by the ensemble machine learning model was 1.4%, 0%, 0.9%, 1.3% (with ST35B) and 0% (without ST35B). These values were determined by application of the ensemble model where we detected 42/2,983, 0/338, 4/457, and 1/78 (or 0/77) of the water samples were >50% likely to be impacted by newly migrated natural gas in these states, respectively. One NY sample, flagged as likely to have newly migrated methane, was sampled ~6 km from the nearest oil/gas well, and thus may be contaminated by new methane from a different source. Fisher exact tests show statistically significant differences between the results for PA versus TX+CO (i.e., with evidence for a higher frequency of migrated methane sites in PA) but show no statistically significant difference within the Marcellus shale play between regions using versus not using high-pressure high-volume hydraulic fracturing. The new machine learning tool appears to be useful in detecting anomalous methane in multiple shale plays. Future work should include additional training and validation of data-driven models as new data from presumably impacted sites in shale plays other than the Marcellus become available.