Numerical modelling of the impacts of water abstraction for hydraulic fracturing on groundwater-surface water interaction: a case study from Northwestern Alberta, Canada

ABSTRACT The temporal dynamics of groundwater–surface water interaction under the impacts of various water abstraction scenarios are presented for hydraulic fracturing in a shale gas and oil play area (23 984.9 km2), Alberta, Canada, using the MIKE-SHE and MIKE-11 models. Water-use data for hydraulic fracturing were obtained for 433 wells drilled in the study area in 2013 and 2014. Modelling results indicate that water abstraction for hydraulic fracturing has very small (<0.35%) negative impacts on mean monthly and annual river and groundwater levels and stream and groundwater flows in the study area, and small (1–4.17%) negative impacts on environmental flows near the water abstraction location during low-flow periods. The impacts on environmental flow depend on the amount of water abstraction and the daily flow over time at a specific river cross-section. The results also indicate a very small (<0.35%) positive impact on mean monthly and annual groundwater contributions to streamflow because of the large study area. The results provide useful information for planning long-term seasonal and annual water abstractions from the river and groundwater for hydraulic fracturing in a large study area.


Snowmelt
The modified degree-day method was used to calculate snowmelt. The governing equation of this method is as follows: where M is the snow melting rate, is the degree-day melting coefficient, P is the precipitation rate, is the current air temperature and 0 is the threshold melting temperature (0ºC was used in this study).

1.2
Overland flow Overland flow was simulated by using finite difference method. This method solves 2-D diffusive wave approximation of the St Venant Equations to estimate the overland flow in the x-and ydirections and the water depths in each grid cell of the model domain. During using the diffusive wave approximations of the St Venant equations, the last three terms of the momentum equations are dropped to reduce the complexity of fully dynamic version of the St Venant equations. The discharge per unit length along the cell boundary in the x-and y-directions are estimated by Eqs (2) and (3), respectively, which are based on conservation of mass and the momentum equations.
where ℎ and ℎ are discharge per unit length along the cell boundary in the x-and y-directions, respectively; h is the flow depth above ground surface, and are Manning's M (inverse of Manning's roughness coefficient n) or Stickler roughness coefficient in the x-and y-directions, respectively; = + ℎ, is the ground surface level, u and v are flow velocities in the x and y directions, respectively. By using Eqs (2) and (3), the surface flow across any boundary between square grid cells is determined by: where Q is the flow across any boundary between square grid cells, K is the appropriate Stickler roughness coefficient, ∆ is the length of a square grid cell, ℎ is the water depth that can freely flow into the next cell, and are the maximum and minimum water levels in the grid cell, respectively.

1.3
Unsaturated flow Unsaturated flow was computed using two-layer water balance method, which is based on a simple mass balance approach. This two-layer method includes interception, surface ponding, water content in the root zone, evapotranspiration (ET), infiltration, and groundwater recharge. The output of this method is the actual ET and the ground water recharge. This method uses meteorological data (i.e., temperature, humidity, radiation, wind speed), vegetation characteristic (i.e., leaf area index and root depth) data, and soil moisture to calculate ET from canopy, ponded water, unsaturated zone (UZ), and saturated zone (SZ). The vegetation properties are described by leaf area index and root depth, and they vary from season to season in a year depending on vegetation types. The soil properties contain a constant infiltration capacity, and the soil moisture contents at the saturation, wilting point, and field capacity. First, this method simulates ET from interception storage in the canopy, and then it simulates ET from the ponded water. This method divides the unsaturated zone (UZ) into (i) a root zone (i.e., ET extinction depth), from which ET can occur and (ii) a zone below the root zone, where ET does not occur. The root zone consists of root depth and ET surface depth, which is the thickness of capillary fringe. If the water table drops below the bottom of the Layer 1, then ET from the saturated zone (groundwater) will stop. This does not mean that ET stops. At that time, ET will continue to remove water from the available water in Layer 1 until soil moisture at the wilting point is reached. The actual ET calculated by this method is the summation of ET from the canopy, ponded water, unsaturated zone and saturated zone: where ET is the actual evapotranspiration; ET can , ET pon , ET uz and ET sz are evapotranspiration from the canopy, ponded water, unsaturated zone, and saturated zone, respectively. If the calculated average water content exceeds the maximum water content in the unsaturated zone, the groundwater recharge starts. The infiltration is calculated by Eq. (6): Infiltration = min {ponded water, (saturated conductivity  duration of a time step), (water content at saturationactual water content)  layer depth} (6)

Saturated flow
Saturated groundwater flow was simulated by using a finite difference representation of the following 3-D saturated groundwater flow equation: where , , are the hydraulic conductivity of soil along the x, y and z axes of the model, which are assumed to be parallel to the principle axes of hydraulic conductivity tensor; h is the hydraulic head, S is the storage term, and W is the flux term for sources and sinks.

2.
Simulation of streamflow in MIKE-11 model Streamflow and stream water level along the channel was computed by using an implicit finite difference scheme to solve dynamic wave version of the St Venant equations (DHI 2017). The governing equations used for this method are as follows: where Q is the streamflow or discharge, A is the flow area, q is lateral inflow, h is the stream water level, C is Chezy resistance coefficient, R is hydraulic radius, is the vertical velocity distribution coefficient. The relationship between the Chezy resistance coefficient and Manning's roughness coefficient n is as follows (DHI 2017

Simulation of exchange of water between streamflow and saturated groundwater flow
The exchange of water between a saturated zone grid cell and the river is calculated by using Darcy's law and the following equation (DHI 2009): where Q is the exchange flow between a saturated zone grid cell and the river, C is conductance, and ∆ℎ is the head difference between saturated zone and river. Conductance is calculated by: where is the leakage coefficient of the bed material, w is the wetted parameter of the crosssection, and dx is the grid size of the model. Figure S1. MIKE-SHE model domain along with MIKE-11 model (stream network) for the study area.

Locations of boundary conditions of coupled MIKE-SHE and MIKE-11 model
Because of unavailability of proper information in the study area for setting appropriate boundary conditions around the perimeter of the study area, a no-flow boundary condition was assumed for the developed MIKE-SHE model (Fig. S2) for model simplicity. In the developed MIKE-11 model, the no-flow boundary condition was considered for all unconnected ends of the river branches (Fig. S3) except, Peace River, Smoky River, and Little Smoky River. Since the upstream parts of these three large rivers (which are situated outside the study area) contribute 88%, 6.1%, and 2.3% of the flow at the outlet of the study area, respectively, based on observed data at the monitoring stations of those area; the boundary conditions at those three upstream parts were considered as inflow boundary. At downstream (i.e., Peace River) of the model (i.e., outlet of the study area), water level boundary was considered as based on the relationship between stream flow vs. water level.