Large-sample evidence on the impact of unconventional oil and gas development on surface waters

The rise of shale gas and tight oil development has triggered a major public debate about such unconventional development, in which horizontal drilling is combined with hydraulic fracturing (HF). HF is the high-pressure injection of water mixed with chemical additives and propping agents such as sand to create fractures in low-permeability formations, allowing oil or gas to flow. Additives in the HF fluids vary with the geological characteristics of the formation and by operator, but the fluid mix usually contains friction reducers, surfactants, scale inhibitors, biocides, gelling agents, gel breakers, and inorganic acid (1, 2). HF wells produce large amounts of wastewater, which initially consists of flowback of HF fluids but over time increasingly consists of produced water from deep formations. The latter brine is naturally occurring water, into which organic and inorganic constituents from the formation have dissolved, resulting in high salt concentrations (1, 3–5).

Although unconventional oil and gas (O&G) development has been important for energy production (6), we do not fully understand the associated environmental and social risks (7–11). These risks include hydrocarbon emissions, water usage, and pollution, along with potential human and ecological health consequences (7, 10, 12–16). Among these, the impact of unconventional O&G development and HF on water quality remains a key concern (2, 8, 9, 17–20). In the US, one reason for this concern is that HF is exempt from the Underground Injection Control provisions of the Safe Drinking Water Act (8). Around the world, unconventional drilling has either just been introduced, or is being considered, by many countries, with uncertain effects on water quality (16).

The US Environmental Protection Agency (EPA) reviewed and synthetized scientific evidence concerning the impact of HF on US water resources. The final report concluded that HF activities can affect drinking water resources under some circumstances (17), but the report did not identify widespread evidence of contamination. Groundwater studies primarily examine contamination from stray gas or deep formation brines, which could occur because of cementing or casing failures or because of brine migration to shallow aquifers through faults or other preexisting pathways (2, 21, 22). Instances of stray gas contamination have been found in Pennsylvania in connection with shale gas development of the Marcellus Shale (23–29) but not in Arkansas for the Fayetteville Shale (30). Geochemical evidence of gas contamination has also been documented for the Barnett Shale in Texas (24) and the Denver-Julesburg basin in Colorado (31). Studies of brine migration from deep formations, mostly in northeastern Pennsylvania, have provided mixed evidence (23, 32, 33). No evidence of brine contamination of groundwater has been documented for the Fayetteville Shale (30). For Pennsylvania, increases in shale gas–related contaminants have been documented at groundwater intake locations of community water systems that are in close proximity to shale gas wells (18).

For surface water, the evidence is more limited. Instances of contamination have been ascribed primarily to discharges of inadequately treated wastewater, HF fluid leaks, and spills and other mishandling of flowback and produced waters (2, 4, 8, 17, 20, 34, 35). Specifically, increased chloride and bromide concentrations downstream of effluents from wastewater treatment plants have been found in western Pennsylvania up to 2011, when the release of wastewaters from unconventional wells into streams through municipal wastewater treatment plants was not prohibited (2, 4, 36). Further, a high frequency of brine spills in North Dakota has resulted in elevated levels of salts and other contaminants in surface waters up to 4 years after the spills occurred (37). A large-sample statistical examination of the effects of shale gas development activities on surface water in Pennsylvania found higher chloride concentrations downstream of wastewater treatment facilities and an association between gas well density in a watershed and increased total suspended solid (but not chloride) concentrations during the period 2000–2011 (38). The authors of this study suggest that both insufficient wastewater treatment and building infrastructure for unconventional O&G extraction could explain the results. Evidence also exists for barium concentrations in Pennsylvania being higher in areas with unconventional wells than in areas without them, but the authors of this study point out that this evidence cannot be solely ascribed to unconventional wells, as it could also reflect the presence of basin brines or a sulfate decrease in acid rains (39). In sum, prior studies document localized instances of surface water contamination related to unconventional O&G development, with spills and leaks being the most common pathway (20).

We investigated the potential impact of unconventional O&G development on surface water quality using a large-sample statistical approach. We combined a geocoded database of 46,479 HF wells from 24 shales with 60,783 surface water measurements over 11 years (2006–2016) across 408 watersheds (HUC10s; HUC, hydrologic unit code) with HF activity (Fig. 1, fig. S2, and tables S1 to S3). Our analysis focuses on concentrations of bromide (Br⁻), chloride (Cl⁻), barium (Ba), and strontium (Sr) in watersheds exposed to unconventional O&G development. We chose these four ions for the following reasons. First, they are usually found in high concentrations in flowback and produced water from HF wells and hence could indicate surface water impact, if and when it exists (1, 2, 4, 5, 8, 20, 27, 33, 36, 39, 40). Moreover, unlike some organic components of HF fluids, these four ions do not experience biodegradation, and their presence has been measured several years after HF spill events (20, 37). Thus, an analysis of ions is a likely mode of detection (2, 4). Second, many of the water quality concerns associated with HF wastewaters are related to the chemistry of the deep formation brines as well as the salinity in flowback and produced water (8). Third, the four ions are measured in many watersheds with reasonable frequency, which is not the case for other potential signatures of HF wastewater. As high salt concentrations can also occur in surface waters for many natural and anthropogenic reasons, such as brine migration or road deicing (2, 4, 33, 39), sufficient data is needed to estimate reliable, local, and time-varying baselines for the background ion concentrations.

We construct these baselines with regression analysis and then exploit temporal and spatial variation in the spudding of HF wells within and across US watersheds to identify anomalous changes in ion concentrations associated with newly spudded HF wells in the same watersheds. Our regression model includes temperature and precipitation control variables as well as an extensive set of fixed effects that construct local and time-varying baselines for background ion concentrations (41, 42). Specifically, our model allows for arbitrary monthly variation in the average background ion concentrations across subbasins (HUC8) and within a given subbasin over time. This flexible regional baseline controls for subbasin differences in water quality, geochemistry, salinity, climate, water body types, or economic activity and also for over-time changes in subbasin concentrations due to seasons, weather patterns and related road deicing, salinization trends, or economic development. Our model also has a local baseline, using each water monitoring station as its own control, which accounts for arbitrary differences in average local ion concentrations. Our model combines these two baselines and the weather control variables to estimate the association between anomalous concentration changes and new HF wells in the same watersheds (42).

Our model explains >80% (in many cases, >90%) of the variation in ion concentrations across watersheds and through time (table S4), suggesting that the model estimates precise baselines for background ion concentrations. We estimated the association between newly spudded HF wells and ion concentrations at the watershed level using a variable that counts the number of HF wells in a watershed at a given point in time (#wellsHUC10). We estimated our regression model for all US watersheds with HF wells and then separately for Pennsylvania, because Pennsylvania accounts for almost 41% of the sample. We found a robust association between new HF wells in a watershed and elevated ion concentrations in its surface waters (Fig. 2 and table S4). For watersheds in Pennsylvania (PA), the coefficients on #wellsHUC10 are positive for all ions and significant for three of them (table S4, column 2, HUC8 model; Br⁻: 0.00019, P = 0.865; Cl⁻: 0.00071, P = 0.031; Ba: 0.00038, P = 0.086; Sr: 0.00041, P < 0.001). The lack of significance for Br⁻ could reflect measurement issues (42). For watersheds throughout the US (ALL), the coefficients on #wellsHUC10 are generally comparable to those for Pennsylvania in terms of magnitude and significance, except for Ba, which is rarely measured outside of Pennsylvania and for which results are weaker (Fig. 2).

To gauge the magnitude of the estimated effects, we multiply each coefficient by the respective sample mean ion concentration and the average number of wells per watershed to obtain the ion concentration increase in the average HUC10 implied by our estimation (HUC10 impact in Fig. 2). With this approach, and focusing on coefficients from the HUC8 specification, we estimated an average increase of Cl⁻ by 1322.44 μg/liter for PA and 2232.55 μg/liter for ALL; Ba by 1.61 μg/liter for PA; and Sr by 5.19 μg/liter for PA and 8.88 μg/liter for ALL (Fig. 2). These magnitudes are very small but need to be interpreted in the context of our analysis. First, the #wellsHUC10 coefficient is by construction a per-well estimate, but this does not imply that each well is associated with a concentration increase, rather it is an average over all wells. Second, by using measurements from all monitors in a given watershed, the estimated well–ion association reflects the average exposure of monitors in the watershed. However, some monitors could be very far away or upstream from wells, in which case they should not be exposed or affected, thus lowering the average. For these reasons, the impact estimates in Fig. 2 are expected to be small; they should increase when the analysis focuses on the most relevant water measurements, as reported below.

We found robust results in different sensitivity analyses (table S5). (i) Estimating separate #wellsHUC10 coefficients for watersheds in and outside of Pennsylvania shows that the findings for Pennsylvania and the other US states are similar. (ii) Estimating separate effects for different time periods shows similar results over time. The coefficients in later periods tend to be smaller but higher in significance, presumably because the frequency of water measurements increases over time. (iii) Our results are similar when the model is estimated over all watersheds within a subregion (including those without HF wells), albeit in some cases slightly weaker, possibly because this specification uses less relevant baselines for the background ion concentrations. (iv) Adding further controls for snow (to account for related discharge of road salts) or for within-watershed seasonality does not alter the findings, suggesting that the monthly subbasin baselines already control for local weather patterns. (v) Our results are also robust to alternative modeling choices, for example, scaling the cumulative number of wells by watershed size, as in (38); estimating the model with weighted least squares (WLS) to give more weight to observations, for which we have more readings to estimate the monthly subbasin baselines; and using alternative transformations to address skewness in ion concentrations.

We also investigated whether the results could be driven by HF-related patterns in the frequency of water monitoring (e.g., more measurements shortly after spud dates or closer to wells). However, we found no evidence that water measurement is systematically related to new wells in a watershed (table S6). We analyzed three other water quality proxies (dissolved oxygen, phosphorus, and fecal coliforms) that are frequently measured but not as indicative of HF-related impacts. Concentration levels of these proxies could be related to other economic activities with water impacts, such as agriculture (43). We used this analysis to gauge how well our model controls for economic activity and other potential confounds. The estimated #wellsHUC10 coefficients for the three analytes are not different from zero (table S7), which contrasts with the results for the ion concentrations we chose for the main analysis (table S4).

Up to this point, our analysis estimated the long-run association between HF wells and ion concentrations, because we did not restrict the sample and the estimation to a particular period after well spudding. However, concentration increases could be stronger early on and fade over time. Hence, we estimated the association in specific time windows around a new well spud date, allowing us to map out the estimates through time. For this temporal analysis, we modified eq. S1 (see supplementary materials) by replacing #wellsHUC10 with several time-specific well counts, defined for the following time windows measured in days: [−180, −91], [−90, 0], [1, 90], [91, 180], [181, 360], and >360. The coefficients were estimated relative to measurements collected 180 days or more before a new well spudding (table S8). We found increases in ion concentrations 91 to 180 days after the spud date, consistently for all four ions, in Pennsylvania and all US watersheds. In Fig. 3, we plotted the coefficients, estimated over all watersheds, for each window together with the 95% confidence interval. For the [91, 180] window, the well count coefficients are significant for all four ions (Fig. 3 and table S8, panel B; Br⁻: 0.01095, P = 0.036; Cl⁻: 0.00401, P = 0.022; Ba: 0.00347, P = 0.017; Sr: 0.00289, P = 0.015). These ion concentration increases in the [91, 180] window are at least one order of magnitude larger than the long-run estimates in table S4 (shown as a red dot in Fig. 3 for comparison). The coefficients for the [91, 180] window (Fig. 3) correspond to an average (short-run) increase of 178.64 μg/liter for Br, 16,014.30 μg/liter for Cl⁻, 15.46 μg/liter for Ba, and 71.34 μg/liter for Sr per watershed with HF wells. Water measurement is often sparse and, as shown earlier, does not increase around the well spud dates. Thus, most measurements naturally fall into the benchmark period of <−180 days (for which no coefficient is estimated) or in the period beyond 360 days after spudding, which explains why the #wellsHUC10 [>360] coefficients always have the tightest confidence intervals. These coefficients are, as expected, comparable to the long-run estimates from table S4 (red dots in Fig. 3). For all other coefficients, we have far fewer observations, and hence their confidence intervals are wider.

It is useful to interpret the increases in ion concentrations between 91 and 180 days after new well spuddings in the context of the well development and HF process. In our samples, the average time span between the spud date and well completion is 103 days, consistent with (31). Thus, the concentration increases occur after the average HF well is completed and during the early phases of production, when large amounts of flowback and produced water are collected. Although the temporal evidence alone does not identify the exact mechanism for the results, it links elevated concentrations to the unconventional O&G development process.

We conducted two tests that further explore the mechanism by focusing on the role of flowback and produced water at the beginning of production. First, we investigated whether the well–ion associations differ depending on the amount of water a well produces after the spud date. We coded two well count variables at the watershed level: (i) the number of wells with an above-median amount of produced water (#wellsHUC10_High_Prod_Water) and (ii) the number of wells with a below-median amount of produced water (#wellsHUC10_Low_Prod_Water). We defined the medians by subbasin and year to account for regional differences as well as potential changes in HF technology over time. For watersheds where wells generate larger amounts of produced water, we found significant coefficients for Cl⁻ and Sr in all specifications, for Ba in Pennsylvania, and for Br⁻ in one specification (table S9, column 1: Cl⁻: 0.00057, P = 0.010; Ba: 0.00047, P = 0.064; Sr: 0.00041, P = 0.012; column 2: Cl⁻: 0.00071, P = 0.031; Ba: 0.00037, P = 0.084; Sr: 0.00041, P < 0.001; column 3: Br⁻ = 0.00035, P = 0.085; Cl⁻: 0.00057, P = 0.044; Sr: 0.00041, P = 0.002; column 4: Cl⁻: 0.00058, P = 0.080; Sr: 0.00037, P < 0.001). These coefficients are larger in magnitude than those of their respective low-produced-water counterparts. Admittedly, the amount of produced water is likely correlated with other well and watershed characteristics, and hence differences in coefficient magnitudes between the two groups do not reflect produced water levels alone. Even so, this evidence ties the elevated ion concentrations more closely to HF activity.

Next, we performed an analysis that exploits subbasin variation in the regional geochemistry and the salinity of deep formations (1, 5, 40). The idea was to explore whether the estimated associations are stronger in areas where HF wells are expected to generate produced waters with higher salinity. We used data from the US Geological Survey Produced Waters Geochemical Database and identified subbasins for which produced waters in previous drillings exhibited total dissolved solids (TDS) concentrations above (below) the median, indicating higher (lower) salinity of the deep formations. We created two well count variables at the watershed level: (i) the number of wells in subbasins with an above-median TDS level (#wellsHUC10_High_Salinity) and (ii) the number of wells in subbasins with a below-median TDS level (#wellsHUC10_Low_Salinity). The coefficients in the high-salinity group are positive and larger in magnitude than the coefficients in the low-salinity group (table S10; Br⁻: 0.00114, P = 0.189; Cl⁻: 0.00120, P = 0.082; Ba: 0.00064, P = 0.003; Sr: 0.00040, P < 0.001). Thus, the association between new HF wells and ion concentrations is most pronounced in subbasins where deep formations exhibit higher levels of salinity. Importantly, natural brine seepage in high salinity areas cannot explain this result, because the model uses the average concentration levels at the monitoring station as its own control. The result implies a (statistical) link between elevated concentrations and HF wells in areas where produced waters have higher salinity, which is consistent with produced water being part of the mechanism.

As wells can be far from the closest monitor (in our sample, average = 10.3 km and median = 8.2 km), we examined whether the associations between HF wells and ion concentrations are more pronounced when wells and monitors are closer together. However, the sparsity of water measurements makes such distance gradient analyses challenging. We therefore estimate results for different distance bins, [0 to 5 km], [0 to 10 km], and so on up to [0 to 30 km], using either Cl⁻ or pooling observations for all ions in one model. In both cases, we found a negative distance gradient (fig. S5), meaning that the coefficients are largest when we estimated them for wells within 5 km of a monitoring station and then steadily decline when we expanded the widths of the distance bins and included wells that are farther away. Owing to sparse data, our confidence intervals are widest when we include only close-by wells but then become smaller as statistical power increases with the expansion of the distance bins.

To address the data sparsity, and concurrently tighten the analysis, we combined dimensions (such as time, distance, and direction of surface water flows) and then either partitioned or restricted the sample. To partition the sample, we defined watersheds as “high type” when there are relatively (i) more water measurements around new well spuddings, (ii) more monitors in close proximity to new wells, and (iii) more monitors that are likely downstream of new wells (37). We estimated an alternative version of eq. S1, in which we replace #wellsHUC10 with two nonoverlapping well counts, counting the cumulative number of wells spudded in high-type and low-type watersheds, respectively. The associations between HF wells and elevated ion concentrations stem mostly from the high-type watersheds (table S11). This finding is reassuring because it shows that the results come from the types of watersheds that a priori have a higher chance of showing a relation between HF wells and ion concentrations if it exists (e.g., because monitors sit closer to HF wells).

We also conducted a separate analysis restricting the sample to the most relevant water measurements. We paired all wells and monitoring stations within a watershed so that we could make determinations for each well–monitor pair with respect to time, distance, and direction of water flows. We then restricted the sample to water measurements that are (i) either taken before or up to 360 days after the spud date, (ii) from monitors within a 15-km radius of the well, and (iii) from monitors for which the paired well likely sits upstream, and we estimated eq. S2 using WLS regressions (42). The variable of interest, Post Spud, estimates the average change in the respective ion concentration after a new well spudding, using each well–monitor pair as its own control and after controlling for temperature, precipitation, and monthly variation in background ion concentrations at the subbasin level. We found elevated ion concentrations in the year after new well spuddings (Fig. 4). In the ALL sample, we estimated positive and significant coefficients on Post Spud for Cl⁻, Ba, and Sr but not for Br⁻ (table S12, panel C; Br⁻: −0.00577, P = 0.499; Cl⁻: 0.02044, P = 0.045; Ba: 0.01585, P = 0.048; Sr: 0.02014, P = 0.098). To gauge the magnitude of the estimates, we calculated the 360-day impact on the average watershed, multiplying the respective coefficient by the mean ion concentration and the average number of new wells spudded per year in the average watershed. For the ALL sample, the average 360-day concentration increases are 47,338.29 μg/liter for Cl⁻, 49.17 μg/liter for Ba, and 282.67 μg/liter for Sr (Fig. 4). These increases are much larger than the long-run increases (Fig. 2) but still well below the EPA maximum contaminant and health advisory levels noted in the Fig. 2 caption. The estimated coefficients are smaller using OLS regressions (table S12, panel D), but even these are still an order of magnitude larger than the respective coefficients for the overall sample (table S4). As a robustness test, we also estimated regressions using observations from “never-treated” monitors that are far away or upstream from any well in the watershed as a control group (table S12, panel E). The results were similar to those in Fig. 4. Notably, the well–monitor pair analysis (Fig. 4 and table S12) showed that estimated well–ion associations strengthened when we focused on the most relevant water measurements.

Our analysis reveals a robust association between new HF wells and elevated concentrations of Cl⁻, Ba, and Sr in same-watershed surface waters. Our results for Br⁻ are generally weaker and often insignificant. The well–ion association was more pronounced when we estimated it using water measurements from monitors that are close to and likely downstream of wells and from the early phases of production, when wells generate large amounts of flowback and produced water. Although this evidence is based on associations, and as such not causal, the regression analysis controls for a large number of observed and unobserved factors and models background ion concentrations in a flexible and comprehensive manner. Thus, we find it difficult to explain the results by other factors. The confluence of our findings is consistent with unconventional O&G development driving the anomalous surface water ion concentrations identified by our model.

Overall, our evidence complements prior studies on the impacts of unconventional O&G development on water quality [e.g., (2, 8)] and extends studies using a similar large-scale approach for the Marcellus Shale (18, 38). Our statistical approach does not reveal the exact mechanism for the association between new HF wells and concentration increases. There are three potential channels that are particularly relevant in our context (2, 4, 8, 12, 17): (i) on-site accidents; leaks and spills of HF fluids, flowback, and produced water, including those related to pits (overflow, unlined pits, tears in liners); brine trucking; and long-term leaching of soils after spills (4, 12, 20, 35, 37); (ii) direct disposal of untreated wastewater from HF operations (unauthorized or permitted) (4, 8, 12, 35); and (iii) inadequate treatment of wastewater at disposal sites (2, 36, 38). We do not distinguish between these specific channels. However, the consistent increases in ion concentrations shortly after well completion, for wells with large amounts of produced water and for wells located in areas with high-salinity deep formations (Fig. 3 and tables S9 and S10) point to the handling of flowback and produced waters, including leaks and spills, being part of the mechanism for the well–ion association we documented. We explored other explanations and found that our results for HF wells do not reflect the presence of facilities accepting HF wastewater in a watershed (table S13), are not observed for conventional wells (table S14), and are still present when we control for a large number of previously documented spills (table S15). Unknown or undetected spills and leaks could still explain our findings. Independent of the exact mechanism, our results are relevant because we show that the association between unconventional O&G development and elevated ion concentrations extends to many watersheds over several US shales.

We acknowledge that the long-run impact estimates we documented using all watershed monitors are very small in magnitude. However, not all wells are close to surface water and not all monitors are in locations where they could detect an effect. In addition, any impact likely declines over time. Consistent with these arguments, we estimated larger associations when we restricted the analysis to water measurements that are taken within a year of well spudding and from monitoring stations that are closer and likely downstream of wells. But even the 360-day concentration increases implied by these estimates do not exceed EPA limits or health advisory levels for the ions. We also note that water measurements are predominantly from rivers, and hence dilution is another factor to consider when interpreting the magnitude of the associations we estimated.

Our statistical approach was constrained by the ions included in public databases as well as the sparsity of water quality data (2, 4, 20, 27, 41, 44). The former means that we could not examine other analytes in HF fluids or related to unconventional O&G development that are potentially more dangerous than salts. The latter implies that we could not perform more granular analyses that would better identify the mechanism or a causal link. Both limitations of our analysis highlight that investigations of surface water impacts from unconventional O&G development would be greatly facilitated if there were more targeted water measurements of relevant analytes in close proximity to and timed around the development of new HF wells.