Changes in Deep Groundwater Flow Patterns Related to Oil and Gas Activities

Large volumes of saline formation water are both produced from and injected into sedimentary basins as a by‐product of oil and gas production. Despite this, the location of production and injection wells has not been studied in detail at the regional scale and the effects on deep groundwater flow patterns (i.e., below the base of groundwater protection) possibly driving fluid flow toward shallow aquifers remain uncertain. Even where injection and production volumes are equal at the basin scale, local changes in hydraulic head can occur due to the distribution of production and injection wells. In the Canadian portion of the Williston Basin, over 4.6 × 109 m3 of water has been co‐produced with 5.4 × 108 m3 of oil, and over 5.5 × 109 m3 of water has been injected into the subsurface for salt water disposal or enhanced oil recovery. Despite approximately equal values of produced and injected fluids at the sedimentary basin scale over the history of development, cumulative fluid deficits and surpluses per unit area in excess of a few 100 mm are present at scales of a few 100 km2. Fluid fluxes associated with oil and gas activities since 1950 likely exceed background groundwater fluxes in these areas. Modeled pressures capable of creating upward hydraulic gradients are predicted for the Midale Member and Mannville Group, two of the strata with the highest amounts of injection in the study area. This could lead to upward leakage of fluids if permeable pathways, such as leaky wells, are present.


Introduction
Large volumes of fluid have been both produced from and injected into the subsurface by the oil and gas industry over the past century. Produced water represents the largest by-product of oil and gas production, and in 2017 the total volume of produced water in the United States exceeded 3.8 billion m 3 (Veil 2020). Volumes of produced water have increased over the past two decades (Lutz et al. 2013;Tiedeman et al. 2016;Scanlon et al. 2017), in part due to the rise in unconventional oil and gas production (Horner et al. 2016;Kondash et al. 2018). It is estimated that 91.5% of this produced water is managed by subsurface injection (Veil 2020) via salt water disposal (SWD) or for enhanced oil recovery (EOR), which we define broadly to include waterflooding. A portion of this increase has been due to injection of surface water or shallow groundwater into the subsurface for EOR, leading to a surplus of water in the deep groundwater systems of several sedimentary basins (Murray 2013;Scanlon et al. 2017;McIntosh and Ferguson 2019;Veil 2020). Note that here we use the base of groundwater protection, 305 m below-ground surface (Province of Saskatchewan 1966Saskatchewan , 2012, as the dividing point between deep and shallow groundwater systems although a range of other definitions for deep groundwater exist (Alley et al. 2014;Tsang and Niemi 2014;Condon et al. 2020). This depth was chosen to protect shallower groundwater supplies with total dissolved solids less than 4000 mg/L from deeper saline fluids that commonly contain hydrocarbons and other contaminants. Similar regulations exist in other jurisdictions but are sometimes based on water chemistry rather than depth (DiGiulio et al. 2018). This surplus of water in deeper formations can lead to increased reservoir pressures driving solute transport  and, in some cases, induced seismicity (Keranen and Weingarten 2018).
Previous studies have focused on regional water balances as proxies for changes in fluid pressures (National Research Council 2013). However, local pressure changes will not necessarily be correlated with changes in the fluid budget at the sedimentary basin scale. Imbalances in fluid budgets have been noted for NGWA.org Vol. 60, No. 1-Groundwater-January-February 2022 (pages 47-63) individual hydrostratigraphic units in both the Western Canada Sedimentary, including the Canadian portion of the Williston Basin (Ferguson 2015) and the Permian Basin (Scanlon et al. 2019), where fluid budgets are approximately balanced at the basin scale. Changes in pressure are possible even where produced and injected fluid volumes are approximately equal at the sedimentary basin scale because of the spatial distribution of wells. The additional fluid fluxes may exceed background regional (i.e., at scales of 10-100 s of km) groundwater flow patterns in deep aquifers. Numerous studies of deep regional groundwater flow have removed pressure measurements affected by fluid production and injection prior to the construction of potentiometric surfaces (Bair et al. 1985;Toth and Corbet 1986;Barson 1993;LeFever 1998;Anfort et al. 2001;Palombi 2008). Most of these studies have used the cumulative interference index developed by Toth and Corbet (1986) to account for the radial proximity of a drillstem test to a production or injection well and remove tests affected by production or injection wells. Other approaches to removing pressure measurements affected by production and injection, including visual inspection (Bair et al. 1985;LeFever 1998), have been used. While such approaches may be useful in understanding background conditions, they may provide little insight into present flow patterns in deeper strata in sedimentary basins with extensive oil and gas development and their potential to create vertical hydraulic gradients necessary for deeper groundwater interact with overlying fresh groundwater supplies.
Here, we examine the spatial distribution of produced and injected fluids in the Canadian portion of the Williston Basin (Figure 1), where fluid budgets at the basin-scale have suggested a net gain in water (Ferguson 2015). This important oil and gas producing region has a relatively long history of fluid production and injection in multiple geologic units with records available back to the 1950s. We show that there are substantial differences between produced and injected volumes at local scales, which may have important implications to groundwater flow and contaminant transport. Induced seismicity has also been related to changes in fluid budgets (Rubinstein and Mahani 2015) and oil and gas development has been linked to seismicity in the Alberta Basin (Atkinson et al. 2016), which is often grouped together with the Canadian portion of the Williston Basin to form the Western Canada Sedimentary Basin (WCSB). The changes in the fluid budget and associated pressure changes related to oil and gas activities may impact the management of produced water and other emerging subsurface uses, such as carbon sequestration, storage of hydrogen and geothermal energy production.

Geology and Hydrogeology of the Williston Basin
The Williston Basin is an intracratonic sedimentary basin made up of near-continuous alternating layers of Late Cambrian to Late Cretaceous age sandstone, carbonate, and shale dominated formations (Carlson and Anderson 1965;Gerhard et al. 1982;Kent and Christopher 1994). The basin is bordered to the west by the Sweetgrass Arch, which separates it from the Alberta Basin, as well as the Black Hills uplift located further south and the Sioux Uplift to the east (Kent and Christopher 1994). The basin has an area of approximately 250,000 km 2 , covering parts of Saskatchewan, Manitoba, Montana, North Dakota, and South Dakota. The Williston Basin has topographic highs in Montana and lows in Manitoba and a maximum stratal thickness of 4900 m (Gerhard et al. 1982).

Hydrogeology
The major regional aquifer systems in the Williston Basin include: the (1) Lower Paleozoic aquifers, (2) Mississippian aquifers, and (3) Mesozoic aquifers (Palombi 2008) (Figure 2). These aquifers also act as oil and gas reservoirs in some regions of the study area (Kent and Christopher 1994). Important regional aquitards include: an Ordovician shale unit that separates the basal Cambro-Ordovician aquifer from the overlying carbonates; the Prairie Evaporite within the lower Paleozoic carbonates; the Bakken Formation, which separates the lower Paleozoic and Mississippian aquifers; and the Cretaceous shales of the Colorado Group and Bearpaw Formation.
Regional flow in all deep aquifer systems down to the Precambrian basement is predominantly southwest to northeast, with recharge occurring in the Black Hills and other areas of high elevation along the southwest and western edges of the Williston Basin (Bachu and Hitchon 1996;Grasby and Betcher 2002). An exception to this is a zone of high density brines in Paleozoic strata near the center of the basin that have stagnated due to negative buoyancy (Palombi 2008;Ferguson et al. 2018). Discharge zones for the basin are found along Lake Winnipeg, Lake Winnipegosis and Lake Manitoba (Grasby and Betcher 2002;Ferguson et al. 2007).
Formation waters in the Williston Basin are dominantly Na-Cl type waters and range widely from 2000 to 350,000 mg/L total dissolved solids (TDS) (Grasby et al. 2000;Woroniuk et al. 2019). TDS values progressively increase with depth and jump significantly below the Prairie Evaporite. Past studies of fluid chemistry and isotopes indicated the presence of a paleo evaporated sea water-derived brine in deeper portions of the basin (Spencer 1987;Wittrup and Kyser 1990;Grasby et al. 2000) and the presence of a brine that originated NGWA.org K. Jellicoe et al. Groundwater 60, no. 1: 47-63 from dissolution of halite closer to the basin edges (Grasby and Chen 2005). Potable waters are generally restricted to within a few 100 m of ground surface and are found in Quaternary sands and gravels along with the shallower extents of Cretaceous and Tertiary sandstones in some regions. Brackish groundwaters, which could be of strategic importance for water security through the use of desalination technologies or to grow salt resistant crops, are likely present in many Cretaceous aquifers at depths of a few 100 m but are poorly mapped in Saskatchewan due to a lack of data (Ferris et al. 2017).

Oil and Gas Production in the Williston Basin
The first commercial gas well in the Williston Basin was established in 1913 in Montana, and the first commercial oil well was established in 1951 in North Dakota (Anna 2013). This marked the beginning of consistent oil and gas production in the U.S. portion of the basin with modern methods, with Canada to follow shortly after with several wells drilled in 1953 (IHS Energy 2020). Since then, oil and gas production has increased as new fields have been discovered, and new technologies were developed.
Between the 1950s and 1990s, oil and gas production was completed solely through vertically drilled wells in conventional plays initially using primary recovery followed by EOR in the early 1960s. Much of the development during this time focused on Mississippian strata, including the Midale Member. During this time, oil production underwent several cycles, a peak in the mid-1960s before a dip in the 1980s, then a steady rise in production until the mid-2010s. This increase in production after the 1980s can be attributed to the expansion into unconventional oil reserves that relied on the development and utilization of horizontal drilling. Horizontal drilling provided more contact with the reservoir reducing the number of wells required and immensely increased well productivity. By 2008, oil producers in the Williston Basin were drilling twice as many horizontal wells than vertical wells, with that number reaching 13 times more in 2016. This activity was primarily in the Bakken Formation and coincided with the rise in the use of high volume hydraulic fracturing (HVHF).

Data Collection
A database of 54,643 oil and gas wells and injection wells was created for a study area that encompassed the bulk of oil and gas production in the Canadian portion of the Williston Basin using data available from AccuMap (IHS Energy 2020). This database included well locations, producing zones, well types, well modes, cumulative and monthly production and injection volumes, well operation dates, and fluid chemistry values. To address the use of different names for the same stratigraphic unit within the study area, producing zones were reclassified based on the stratigraphic column for southeastern Saskatchewan (Saskatchewan Ministry of Economy 2014).
To assess the spatial distribution of produced and injected volumes, the study area was broken down into a grid of 5 km × 5 km cells. This grid size was used to allow comparison of different formations chosen based on well spacing in the region. In some areas in the Midale Member where infill drilling has occurred, spacing can be as little as <0.2 km 2 , while unconventional reservoirs typically have well spacing of ∼2 km 2 and spacing can exceed 10 km 2 for disposal wells (IHS Energy 2020). The number of wells present within each cell was counted and produced oil and water and injection volumes for every well in the cell were summed to calculate the net injected volume. To facilitate comparison to background regional groundwater flow, net injected fluid values are presented as a millimeter.
Injection rates were also compared to the background regional groundwater for selected hydrostratigraphic units. Background flow rates occurring in each hydrostratigraphic unit were estimated using Darcy's Law for fluid flow through a porous media: where Q (m 3 /s) is the volumetric flow rate, k is the permeability of the reservoir, A (m 2 ) is the cross-sectional area of the reservoir and grid cell, h (Pa s) is the change in pressure head, and L (m) is the length.

Reservoir Pressure Changes
Analytical models of pressure changes in the 5 km × 5 km cells with the largest cumulative injected volumes were created for selected hydrostratigraphic units. This was done to understand the spatial variability of the pressures of fluid production and injection at a local scale in a manner that cannot be resolved from the gridded water budget or a numerical model based on those fluxes. The models are necessary to estimate pressure changes in the region, which are not available from government regulators in the study area or IHS Accumap, which largely draws its data from government sources. This approach also allows for estimation of pressure changes where wells are not available for monitoring.
Simulations were produced with Aqtesolv (Duffield 2007), which uses the Theis equation: where s(t) (m) is the drawdown over time, Q (m 3 /d) is the well pumping or injection rate, r (m) is the radial distance from the well, S (dimensionless) is the storativity of the reservoir, and T (m 2 /s) is the transmissivity of the reservoir. Aqtesolv allows for variations in pumping over time through the use of superposition of transient responses.
These simulations assume single phase flow, fully penetrating wells, homogenous transmissivity and storativity, and that each aquifer is a uniform thickness. While it is clear that assuming single phase flow may not be strictly appropriate due to the presence of oil and gas, it is assumed that the much smaller volume of oil will not have a significant impact on the overall formation pressures created by the injection of waters. In reality, this will result in underestimation of pressure responses due to the decreases in permeability linked to the relative amounts of oil and water in the pore space (Parker 1989). Horizontal wells have also been treated as producing or injecting from a single point instead of along the entire length of the horizontal perforations.
Transmissivity can be calculated by converting permeability data into a hydraulic conductivity using an Equation (4) developed by Muskat (1937) that relates Darcy's permeability and the weight of a fluid: where k (m 2 ) is the permeability, K (m/s) is the hydraulic conductivity, μ (Pa s) is the dynamic viscosity of the fluid, ρ w (kg/m 3 ) is the density of the fluid, and g (9.81 m/s 2 ) is the gravitational acceleration constant. Once a hydraulic conductivity has been calculated, it is possible to calculate transmissivity using the following equation: where T (m 2 /s) is the transmissivity of the reservoir, K (m/s) is the hydraulic conductivity of the reservoir, and b (m) is the thickness. Storativity values were unavailable for the study area because this parameter is not typically used by the oil and gas industry and cannot be easily determined from single well tests. For a confined aquifer, storativity can be calculated with the following equation: where S (dimensionless) is storativity, S s (m −1 ) is the specific storage, and b (m) is the thickness. The specific storage (S s ) (m −1 ) is the volume of water removed from a unit volume of a confined aquifer per unit drop in hydraulic head. It is related to the compressibilities of the aquifer and the fluid.
where ρ w (kg/m 3 ) is the density of the reservoir fluid, and g (9.81 m/s 2 ) is the gravitational acceleration constant, α is the aquifer compressibility (Pa −1 ), n (dimensionless) is the aquifer porosity, and β is the compressibility of the reservoir water. The aquifer compressibility was based on values presented in related literature on the characteristics of target strata (Beliveau 1989).
Changes in hydraulic head ( h) for individual production and injection wells were calculated separately. Overall changes were calculated by invoking the principle of superposition and adding the individual results together. These hydraulic head changes were converted to reservoir pressures ( P ) (Pa) changes, using the following equation: where g is the gravitational acceleration constant and ρ w (kg/m 3 ) is the density of the reservoir fluid.

Distribution of Wells and Produced and Injected Fluids Well Counts and Fluid Volumes
Within the study area, there are more than 54,623 wells total, including 28,100 active oil and gas wells, 2890 injection and disposal wells, 15,339 abandoned wells and 7104 suspended wells, with the remainder having a range of other classifications. Four hydrostratigraphic units contain the majority of the wells, led by the Midale with 8428, the Bakken with 6718, the Frobisher with 6763, and the Tilston with 2703 wells.
While historically vertical wells were common, the use of HVHF and horizontal wells has become prevalent, beginning in 2005, as the rapid adoption of new technologies made it possible to produce from lower permeability strata ( Figure 3). Within the study area there are 28,651 vertical wells, 22,214 horizontal wells, and 1153 directional/deviated (dir/dev) wells ( Figure 4).
Over the last 60 years within the study area, a total of 540 million m 3 of oil, and 51,000 million m 3 of gas production has been reported. The quantity of produced water is nearly 10 times the amount oil produced, at nearly 4600 million m 3 , and 5500 million m 3 of water was injected into hydrostratigraphic units for water flooding or as SWD. Injected volumes provided for the region do not include those associated with hydraulic fracturing. These were not part of the data available during this study.
The bulk of fluid production and injection has occurred within the Mannville Group, the Madison Group (Poplar, Ratcliffe, Midale, Frobisher, Kisbey, Alida, Tilston and Lodgepole) and the Bakken Formation. Activities in these strata are related to production of oil and associated management of flowback and produced waters. Considerable quantities of water are also being injected into the Interlake, Stonewall, and Deadwood formations to dispose of brines produced at potash mines in the region.  almost no net change of fluid volume within those strata. The Bakken Formation is notable due to its relatively recent production history and its use of HVHF as a primary extraction method.
Between 1952 and mid-2019, the Mannville Group has produced a total of 6.1 million m 3 of oil, 376 million m 3 of water and a negligible volume of gas. Additionally, there have been 726 million m 3 of produced water injected into this hydrostratigraphic unit, creating a surplus of 344 million m 3 of fluid ( Figure 5). Widespread use of the Mannville Group as a disposal reservoir did not begin until the late 1990s, after which it only took a short period of time for injected water volumes to surpass monthly produced water rates. Within 10 years of injecting, the cumulative volume of injected water began to exceed the volume of produced water. The volume of produced water from the Mannville Group has stayed relatively constant over the active lifetime of the formation.
The Midale Member has produced a total of 215 million m 3 of oil, 860 million m 3 of water and 26,000 million m 3 of gas between 1953 and mid-2019. Additionally, there have been 1050 million m 3 of water injected, resulting in a loss of 25 million m 3 of fluid within the Midale Member during the same time period ( Figure 5). Oil production rates in the Midale Member have remained relatively stable since the beginning of production in the 1950s. While there was significantly more water injected than produced between the 1960s and 1980s, after this period, these rates become comparable and follow a similar trend. A flattening of injection and production rates in the Midale Member occurred around 2010, coinciding with the rapid increase in oil production from the Bakken Formation.
The Bakken Formation has produced a total of 50 million m 3 of oil, 110 million m 3 of water and 6400 million m 3 of gas between 1956 and mid-2019, with the bulk of this production occurring since 2008. Additionally, there have been 22 million m 3 of produced water injected into the formation, resulting in a loss of 140 million m 3 of fluid within the formation ( Figure 5). The volumes of produced and injected water within the Bakken Formation are significantly lower than those found in the Mannville Group and the Midale Member, however, produced oil rates are comparable to rates in the Midale Member.

Spatial Variability of Production and Injection
Well densities exceed 100 wells/25 km 2 over much of southeastern Saskatchewan and in a small area of southwestern Manitoba (Figure 3). Well densities exceeding 300 wells per km 2 are found locally within this region.
Wells within the Mannville Group are fairly spread out when compared to well densities in other hydrostratigraphic units (Figure 3). The average well density per cell is only 2.8 wells/25 km 2 , while the maximum density is 147 wells/25 km 2 . The sparseness of wells can be attributed to the fact that 67% of active wells in the Mannville Group are disposal wells, with most of the rest being source water wells. These disposal wells often service many surrounding production wells requiring a smaller number of wells with larger spacing.
Wells in the Midale Member are tightly grouped within the middle of the study area and extend toward the Canada-U.S. border (Figure 3). The Midale Member has the highest average density of wells of all strata examined, at 31.1 wells/25 km 2 , as well as the highest number of wells in one cell at 299 wells. This high density of wells is due to the use of water flooding within this hydrostratigraphic unit, leading to a nearly equal number of injection and production wells in many cells. However, production wells (3558) are more common than injection wells (861) in the study area. EOR wells are not classified separately from other injection wells within AccuMap and some injection wells are likely SWD wells operating in nonproductive or previously productive areas of the reservoir. This is reflected by the mean cumulative fluid production of 187,000 m 3 per well and a mean cumulative injection of 738,000 m 3 per well, which deviates from the nearly equal production and injection rates typical of waterflooding.
The wells in the Bakken Formation primarily cluster into three main groupings, one in the southern portion of the study area, one in the center, and one toward the eastern edge extending into Manitoba (Figure 3). The average well density in the Bakken is 17.3 wells/25 km 2 , with a maximum of 191 wells. These wells are nearly all production wells. Wells are primarily focused in the center of each one of these groupings, with the density of wells decreasing outward.
In the study area, the cumulative (i.e., over the history of oil and gas development in the basin) maximum increase in fluid volume per unit area per 25 km 2 cell was 8995 mm (Figure 6). The maximum cumulative decrease in fluid volume per unit area was 3315 mm, and the average change across all cells was 3.9 mm.

Injection Compared to Regional Flow Rates
To quantify the volume of fluid being injected into reservoirs, the annual net fluid budget increase for individual cells was compared to the estimated natural flow rates for each hydrostratigraphic unit. Hydraulic gradients were estimated from hydraulic head maps published by Palombi (2008), who used an algorithm to omit measurements affected by production and injection. The background hydraulic gradient in the Midale Member was 8 × 10 −5 and 2 × 10 −4 in the Mannville Group. Using the background hydraulic gradients and estimated transmissivities, natural flow rates for each cell in the Midale Member were between 1.3 and 130 m 3 /year, while the Mannville Group was estimated at between 320 and 32,000 m 3 /year. The significantly larger value for the Mannville Group compared to the Midale Member is due to it having a higher permeability and being 10 times thicker. In the Midale Member, cumulative net injection values exceeding an absolute value of 50 mm per unit area are found in some areas reaching over 100 km 2 in size (Figure 6), which is equivalent to more than 18,000 m 3 /year over the period from 1950 to 2019. These fluxes, which arise from a combination of EOR and SWD, greatly exceed background flows and will be the strongest determinant of groundwater flow direction. Cumulative injected volumes per unit area exceeding 100 mm are common in the Mannville Group ( Figure 6). This corresponds to an average annual flow rate of over 36,000 m 3 /year over the period from 1950 to 2019, suggesting that disposal has become the dominant fluid flux over part of the study area. The average cumulative injection for all 25 km 2 in the Mannville is 1430 m 3 /year, placing it within the range of estimated values for background fluxes of groundwater.

Estimation of Reservoir Pressure Changes
Reservoir pressures were simulated for selected regions of the Midale Member and Mannville Group within the larger study area using the Theis equation (Equations 2 and 3). Fluid flow in the Bakken Formation was not simulated because long-term pressure increases will not occur due to a lack of EOR and SWD. These

Figure 6. Cumulative differences in produced and injected volumes. (A) Cumulative production and injected volumes for every well showing the complexity of the system. (B) While wells in the Mannville Group are primarily used for injection, there are still many areas that have produced significantly more water. The cell with the maximum cumulative injected volume (1026 mm) is located near the northern extent of the area where injection in the Mannville Group is common. (C) Produced and injected volumes within the Midale Member are almost identical, however there are still regions where the difference in these volumes is quite large, reflecting the presence of both production and injection wells associated with EOR. The cell with the maximum injected volume (497 mm) is located in the western portion of the developed area of the Midale Member. (D) Since most Bakken Formation wells are hydraulically fractured there is little injected water resulting in every cell having a negative fluid volume change.
transient simulations were centered on areas with the highest cumulative injected volumes and considered the cumulative effects of multiple wells.

Midale Member
A model was created for the region of the Midale Member containing the injection well with the largest injection volume. This well injected a total of 6.6 × 10 6 m 3 of water between 1963 and 2019 at an average rate of 320 m 3 /day (Figure 7), primarily for the purpose of creating a pressure gradient to drive oil toward production wells. An area of 5 km by 5 km was chosen to allow for an adequate number of wells to interact with the primary injection well. Within the modeled area, there was a total of 124 wells, 110 have been primarily for oil production, and 14 are used primarily for injection. Since 1957 a total of 50.2 × 10 6 m 3 of fluid was produced, and 42.4 × 10 6 m 3 of fluid was injected in the modeled area. While the total volume of fluid produced is greater than the volume injected, the daily rates follow similar trends (Figure 7). Changes in hydraulic head were calculated for all wells using Aqtesolv, which uses a deconvolution approach to allow for variable pumping and injection rates. The overall change in hydraulic head from all wells was then determined using the principle of superposition.
Transmissivity and storativity values were estimated from compressibility, permeability and thickness values (Table 1).
Injection wells in the modeled area were arranged in a grid pattern along lines that ran NW-SE and SW-NE ( Figure 8). Spread among these injection wells are production wells that follow no organized pattern. Toward the SW corner of the modeled area is the outer edge of the oil field; due to this, there are no injection wells in this area.
For a transmissivity of 10 −7 m 2 /s, local pressure changes exceeding 100 MPa ( h > 10,000 m for fluid density of 1000 kg/m 3 ) were present around injection wells by 2019, with decreases of similar magnitudes around production wells (Figure 8). These pressures are unrealistically high and would likely lead to hydraulic fracturing. At transmissivity of 10 −6 m 2 /s, pressure increases of ∼2 MPa ( h ∼ 200 m) were estimated over much of the northeastern portion of the area simulated. Assuming near hydrostatic initial conditions, this would result in hydraulic head values well above the ground surface and large upward hydraulic gradients (>0.1) would be present in this region. A similar pattern emerges for a transmissivity of 10 −5 m 2 /s however pressure increases are limited to ∼0.2 MPa ( h ∼ 20 m) over most of the simulated area.

Mannville Group
A model for an area of the Mannville Group centered on the injection well with the largest injection volume. This well had injected a total of 9.4 × 10 6 m 3 of water between 1997 and 2019 at an average rate of 1172 m 3 /day. Within the modeled area, the Mannville Group is almost exclusively used for disposal of produced waters and contains a total of 45 disposal wells. In addition, there are two water production wells for use in EOR and HF operations. Since 1997 a total of 0.3 × 10 6 m 3 of fluid was produced, and 82.3 × 10 6 m 3 of fluid was injected in the modeled area (Figure 9). There was zero produced fluid in the study area until late 2018 (Figure 9). Transmissivity and storativity values were estimated from existing hydraulic conductivity, compressibility, permeability and thickness values (Table 1).
In the model with transmissivity of 1 × 10 −5 m 2 /s pressure increases were largely restricted to the southwestern portion of the simulated area ( Figure 10). By 2019, an area where pressures increases of >5 MPa (or equivalent to h of approximately 510 m for freshwater) were predicted near a cluster of disposal wells. When transmissivity was increased to 1 × 10 −4 m 2 /s, maximum pressure increases greater than 0.4 MPa ( h ∼ 43 m) were present over most of simulated areas and increases greater than 2 ( h ∼ 204 m) MPa were restricted to a few km 2 in the southwest. For the simulation with a transmissivity of 1 × 10 −3 m 2 /s, we find a maximum pressure change of 0.44 MPa for 2019 and increases in pressure > 0.1 MPa ( h ∼ 10.2 m) throughout the modeled area. The resulting upward hydraulic gradients are similar to those estimated for the Midale but would be present over larger areal extents.

Fluid Volume Changes
Examining fluid budgets at the basin scale or even individual hydrostratigraphic units is not sufficient to provide an adequate understanding of the impacts of fluid production and injection on fluid pressures and groundwater flow. Net changes of fluid volumes can be near zero at these scales, masking important local changes in fluid budgets and pressures. When different hydrostratigraphic units are lumped together, regions with substantial increases in fluid volumes can exist directly next to those with large decreases in fluid volumes (Figure 8).
The Mannville Group is used primarily for SWD and thus shows substantial increases in fluid volumes for most of the cells. Areas where the Mannville Group is used as source water for water flooding, can be seen in a sizeable region toward the southern extent of the Mannville Group wells where there are net cumulative water losses ( Figure 8B). This area of decreased fluid coincides with a location of the Midale Member that has increased fluid volumes, suggesting that within this area water from the Mannville Group is used for water flooding in the Midale Member. This provides an example of the movement of water between hydrostratigraphic units that may be common both in the Williston Basin and other similar environments.
Although the difference in Midale Member produced and injected fluids is insignificant for the study area as a whole, regions of increased fluid volumes, as well as decreased fluid volumes exist. This is partly due to water flooding, as some cells see more water injected into them to drive oil into the neighboring cell and such differences may disappear if larger cells were used or were shifted slightly in space. The majority of cells have a minor amount of fluid volume change, suggesting that in most cases the amount of water injected for water flooding is similar to the amount of water produced, effectively canceling out the net effect of either. Areas with larger increases in fluid volumes likely reflect the presence of SWD wells in the Midale Member.
Oil and gas production in the Bakken Formation is primarily conducted using HVHF which uses substantially less water for production compared to the water flooding occurring in the Midale Member. This has created a decrease in fluid volume for every cell within the Bakken Formation. Fluid decreases in the Bakken Formation are lower and more consistent than those found in the Mannville Group and the Midale Member.

The Effect of Production and Injection on Reservoir Pressures
While assessing produced and injected water volumes on a hydrostratigraphic unit basis is important, it is also imperative to understand how the spatial distribution of production and injection wells affects the subsurface flow and formation pressures. Reservoirs where the volumes of produced water equal injected water can still have regions where there are differences in produced and injected volumes spatially. These differences can be caused by areas where only conventional production is occurring, and no water is being injected, areas where there is no oil present and SWD is being utilized, and during EOR where patterns of injection wells are used to push water and oil toward a series of producing wells. These localized differences in fluid budgets can be seen in the Midale Member where there are negligible differences in produced and injected water, yet it has a grid cell with a water surplus per unit area of 416 mm directly next to a grid cell with a deficit of −255 mm (Figure 8).
Areas with extensive water injection will see increases in formation pressures, while areas with more production will see a decrease in pressure. The pressure will fluctuate with the distribution of production and injection wells, sometimes over as little as 100 m. These pressure changes can potentially act as drivers of fluid flow and could lead to contamination of overlying freshwater resources where high permeability pathways, such as leaky wells, are present (McIntosh and Ferguson 2019). These increases in pressure are necessary to drive upward flow in the study area, where underpressures and downward hydraulic gradients are common (Palombi 2008). Modeled pressures in the Midale Member showed an increase of >8 MPa up to 250 m away from the injection well, and 2 MPa increases at up to 1.5 km away (Figure 7) where transmissivities less than 10 −6 m 2 /s were simulated. Larger pressure increases, such as those estimated using a transmissivity of 10 −7 m 2 /s, are unlikely because they would result in hydraulic fracturing. Simulated pressure increases are enough to drive the hydraulic head in the Midale Member far above the ground surface, creating the possibility of near-surface groundwater contamination. These upward gradients will be focused in areas of the Midale Member where injection exceeds production. Large areas of the Midale Member have cumulative losses of fluid (Figure 8), and the associated reductions in fluid pressures may lead to downward hydraulic gradients. Injection rates were greater in the Mannville Group where the combination of a lower well density and higher transmissivity resulted in slightly lower estimated increases in pressure, although these increases propagate over larger areas. The importance of hydrogeologic properties and distribution of wells emphasizes that net injection will not be able to provide reliable estimates of increases in fluid pressures or changes in flow direction.
Any contamination of shallow groundwater systems will require the presence of a high permeability pathway. Over much of the study area, low permeability Cretaceous shales will impede the upward movement of fluids and protect overlying groundwater systems (Shaw and Hendry 1998;Hendry et al. 2013). There is some evidence of vertical faulting through these strata (Smith and Pullen 1967;Gendzwill and Stauffer 2006;Szmigielski and Hendry 2017) but there are no permeability data available for these structures. Leaky well bores may also provide a pathway for contaminants to reach shallow groundwater systems. Wells that leak methane are known to occur in the study area (MacKay et al. 2019) and are common in other areas of western Canada with similar well construction and abandonment regulations (Bexte et al. 2008;Watson and Bachu 2009). However, there are no comprehensive studies from the study areas that have evaluated the potential for transport of water from below the base of groundwater protection to shallow aquifers. There are no documented cases of contamination of this type, which could be attributable to a lack of events or the sparse groundwater observation network of ∼70 wells in Saskatchewan (Saskatchewan Water Security Agency 2019) or a combination of the two. Given the concern about the possible impact of HVHF on potable groundwater resources (Birdsell et al. 2015;Brownlow et al. 2016;DiGiulio et al. 2018;McIntosh et al. 2018), the presence or absence of contamination of potable groundwater supplies from injection wells operating over much longer time periods should provide some insight into the relative risk of HVHF and other subsurface activities involving injection.
While our analysis focuses on the Canadian portion of the Williston Basin, similar changes in the fluid budget and pressures maybe occurring in the Dakota Group, which is equivalent to the Mannville Group, in the United States. 1.23 million m 3 of flowback and produced water was injected into the Dakota Group between 2005 and 2014, resulting in concerns about rising fluid pressures . The Madison Group has produced over 160 million m 3 of oil in the United States (Gaswirth et al. 2010), which is less than the over 400 million m 3 of oil produced from equivalent strata in Canada. However, conventional developments may still have important effects on water budgets and changes in groundwater flow patterns locally on the American side of the Williston Basin.
Injection associated with EOR and SWD are more likely to drive solute transport into overlying strata due to the long time periods involved. Unlike hydraulic fracturing, which only increases formation pressures for a few days, EOR and SWD wells can operate for decades . SWD will be associated with increases in formation pressure that are higher and more spatially extensive than EOR, due to the fluid production associated with EOR that will act to balance out pressure increases at larger scales. Elevated pressures due to SWD can persist for >10 years even after considerable reductions in injection rates (Pollyea et al. 2019). The resulting increased hydraulic gradients will also increase contaminant transport distances where high permeability pathways are present. In areas where natural fractures or faults or other breaches in aquitards are absent, leaking abandoned oil and gas wells may act as conduits for contaminant transport. Contamination from leaking wells is not a new phenomenon (Eger and Vargo 1989;Dusseault et al. 2000), and potential increases in contaminant transport distances associated with prolonged injection will only exacerbate the problem.

The Role of Pressure in Induced Seismicity
Over the last decade, there has been an increase in induced seismicity due to activities associated with oil and gas development Keranen and Weingarten 2018). Induced seismicity has been directly tied to several hydraulic fracturing operations, including instances in the western portion of the WCSB (Atkinson et al. 2016), but no such incidents have been documented in the Canadian portion of the Williston Basin. Many induced seismic events in the United States are due to SWD into deep reservoirs (Rubinstein and Mahani 2015). A substantial change in the net fluid budget is thought to be the largest influence of seismicity (NRC 2013). It is estimated that an increase as little as 0.01 to 0.2 MPa along faults or tectonically stressed features can induce seismicity (Keranen et al. 2014;Hornbach et al. 2015). Our findings emphasize that net changes in fluid budgets cannot be related to fluid pressure changes without considering the distribution of injection and production wells and the hydrogeologic properties of the reservoir.
Induced seismicity as a result of large-scale SWD has been extensively studied in Kansas and Oklahoma. In this region, high volumes are injected into the Arbuckle Group, a thick sedimentary reservoir overlying the crystalline basement. An annual total of 16 million m 3 of wastewater was injected into the Arbuckle Group in 2015, leading to recorded reservoir pressure increases of up to 0.4 MPa by 2016 (Peterie et al. 2018). This pressure led to a record number of Magnitude 3 earthquakes recorded by the US Geological Survey (USGS) in 2014. Similarities can be found in the Mannville Group reservoir, which is primarily used for SWD. In the Mannville Group in 2015, there were more than 43 million m 3 of wastewater injected, compared to 163 million m 3 in the Arbuckle Group over an area of similar size (Scanlon et al. 2019). Our estimated increase in the pressure of up to 0.4 MPa around injection wells in the Mannville Group is also similar to modeled pressure increases in the Arbuckle Group (Keranen et al. 2014). Despite these parallels, there are no recorded events of induced seismicity associated with injection into the Mannville Group. Reservoirs resembling the Mannville Group are unlikely to exhibit induced seismicity due to the following factors: (1) aquifer properties allow for large influxes of wastewater without a significant increase in reservoir pressure; (2) a lack of faults and low-stress levels; and (3) distance from the crystalline basement rocks. The stability of the underlying Precambrian basement and low frequency of seismicity historically may also contribute to a lack of induced seismicity in the region (R. Horner and Hasegawa 1978). We do note that there have been seismic events in the northeastern portion of the study area, with eight events exceeding M3.0 occurring since 2009 (Geological Survey of Canada 2020). The locations of these events coincide with an area where injection into the Interlake Group occurs without any fluid production.

Anthropogenic Evolution of Flow
Most maps of potentiometric surfaces for deep aquifers assume steady-state conditions and are commonly constructed by using pressures measured from drillstem tests over periods of decades and remove data that is considered influenced by injection and production wells (Toth and Corbet 1986;Barson 1993). As the results of this study show, fluxes from production and injection wells may exceed those associated regional groundwater flow. Surplus produced water is being reinjected into some reservoirs at rates significantly higher than natural flow rates, with rates as high as 6000 times natural flow rates present in the Midale Member, and up to 30 times in the Mannville Group. While the analytical models presented here examined the areas with the highest cumulative net injected volumes, areas of elevated pressure are likely common in both the Midale Member and the Mannville Group. The pressure anomaly scales linearly with production and injection rates where material properties are the same, indicating that hydraulic head anomalies on the order of 10 to 100 s of m will be present locally around injection wells used for EOR in the Midale and that anomalies of 1 to 10 m will be widespread in the Mannville Group. Regional hydraulic head patterns in areas of extensive oil and gas development are unlikely to resemble background conditions. Previously mapped underpressures (Palombi 2008) that screened out effects of production and injection could be replaced by conditions where upward hydraulic gradients may allow for migration of saline waters and hydrocarbons to overlying Tertiary and Quaternary strata hosting domestic and agricultural supply wells. Contamination would also require that a sufficiently high permeability pathway would allow for transport during the time period (∼decades) when these pressure anomalies would exist.
The shift in fluid pressures will make it difficult for projects, such as carbon sequestration or geothermal energy production that rely on hydraulic head maps to estimate pressures or groundwater flow velocities. To accurately predict reservoir conditions, modeling will need to be conducted at each project site, or new basinwide hydraulic head maps will need to be created with the effect of oil production included. To create basin-wide models of potentiometric surfaces, it will be necessary to model every oil and gas well in the basin and changes in pressure over time will need to be considered. To improve long-term forecasting of these potentiometric surfaces, well pressures and volumes are required. Currently, reservoir pressure measurements are only reported before production starts creating a sparse dataset for the Canadian portion of the Williston Basin. Improved characterization of long-term changes in fluid pressures and groundwater movement in sedimentary basins are required to facilitate emerging uses of the deep subsurface.

Conclusions
Fluid budgets in the Canadian portion of the Williston Basin have experienced large changes locally despite the presence of similar produced and injected volumes at the regional scale. Use of nonproducing formations for SWD has resulted in increases in the fluid budgets for some hydrostratigraphic units and reductions for others. In cases where EOR is common and fluid budgets are approximately balanced at the regional scale, notable local variations can occur due to the distribution of production and injection wells. These local variations will lead to changes in fluid pressure that will drive solute transport and potentially lead to induced seismicity. We expect to see similar fluid pressure changes in other sedimentary basins containing stacked reservoirs with complex development histories.
To further improve the understanding of the changes in fluid flow rates and directions in deep regional aquifers due to oil and gas production, it is recommended that additional data need to be collected to supplement currently available data. Increasing the required number of fluid pressure measurements for each well and requiring the reporting of the source of injected waters will improve the ability to predict changes in subsurface pressures. Expanding the number of available hydrogeological measurements (permeability, compressibility, porosity, etc.) can increase the accuracy of pressure predictions. Implementing the collection of multiple fluid chemistry measurements over a well's lifespan instead of only at the time of completion could be useful in tracking the movement of injected fluids. Collection of these data will require significant participation from industry and will likely require additional regulations.
The shifts in groundwater flow from production and injection in deep strata may have unknown consequences. Deep groundwater flow systems are generally poorly characterized (Alley et al. 2014;Tsang and Niemi 2014) and anthropogenic effects are not well integrated into current characterization efforts. In many instances, the saline fluids in these deep environments are effectively disconnected from the rest of the hydrologic cycle under background conditions (Palombi 2008;Ferguson et al. 2018;McIntosh and Ferguson 2021). Our results demonstrate that fluid pressures and groundwater flow directions in areas with extensive oil and gas development will experience substantial deviations from background conditions, which have often been screened out in some previous studies of deep aquifers. The resulting vertical hydraulic gradient could allow for connection of deep and shallow groundwater systems where high permeability pathways, such as leaky wells are present. Future use of the deep subsurface will need to consider these changes in pressure, as well as those that may arise from future uses. are grateful to Matthew Lindsay, Christopher Hawkes and Roger Beckie for reviews of an earlier version this manuscript. Feedback from two anonymous reviewers greatly improved the quality of this manuscript.