1. Introduction
The Alberta No. 1 (ABNo1) geothermal project, located within MDGV, just south of the City of Grande Prairie, will be the province’s first conventional geothermal electrical energy producing facility (
Hickson et al. 2021;
Huang et al. 2021b). Grande Prairie is the largest city in the region with a population of 55 000 spread over 72.8 km
2 (
Banks 2016). According to
Hickson (2022), ABNo1 is expected to produce 80 000 MW·h of electrical power per year, offset 96 000 tCO
2/year, and produce 985 TJ/year of thermal energy. Nearby facilities include the Grovedale Light Industrial Park and Greenview Industrial Gateway. In its first phase, the project plans to drill three production wells and two injection wells. In ABNo1, the primary target formation is a Devonian carbonate unit consisting of interbedded dolomite and sandstone units (
Hickson et al. 2021;
Huang et al. 2021a;
Huang et al. 2021b).
One of the technical/environmental issues that needs to be addressed before starting a deep geothermal development in Alberta is the potential for induced seismicity due to local fault reactivation as a result of fluid injection (
Atkinson et al. 2016;
Bao and Eaton 2016;
Schultz et al. 2017;
Mahbaz et al. 2021). Upon injection of fluid into the fractured rock mass, pore pressure is increased along pre-existing faults and fracture planes; if the resolved shear stress on the plane exceeds the normal effective stress limit, the fault slips and thereby triggers an earthquake (
Yaghoubi 2019). Anthropogenic seismic events in Alberta arising from oil and gas operations are among the largest Mw (moment magnitude) events reported globally (
Atkinson et al. 2016), including events near Fox Creek at Mw 4.1 on 12 January 2016 (
Schultz et al. 2017), and Musreau Lake at M
L 3.94 (local magnitude) on 25 December 2019 (
Li et al. 2021). In terms of distance, Fox Creek is approximately 200 km southeast of the ABNo1 development site. Musreau Lake is approximately 100 km south (
Fig. 1) of ABNo1. Thus, at the beginning of this type of project in Alberta a precautionary assessment is prudent.
Prior to the development of a deep geothermal project, an injection-induced earthquake assessment should be performed. For example, in Switzerland, The Deep Heat Mining Basel Project implemented in a densely populated suburb of Basel was shut down due to an M
L3.4 injection-induced seismic event that occurred on 8 December 2006. The project was initiated for geothermal reservoir enhancement from the granitic basement (“hot dry rock”). During 6 days of hydraulic fracturing (HF) stimulation, approximately 13 000 induced seismic events were detected at a depth of 5–6 km (
Häring et al. 2008;
Mukuhira et al. 2013;
Kraft and Deichmann 2014). Even though the geologic and operational characteristics of the ABNo1 and Basel projects are different, it is still important to apply the lessons learned from Basel to the assessment of the likelihood and risk of injection-induced seismicity in ABNo1.
Two kinds of parameters impact the rate and magnitude of injection-indued seismicity. First, extrinsic parameters that can be controlled, such as injection pressure, flow rate, viscosity, volume, and type of fluid, and second, subsurface intrinsic parameters that are uncontrollable. These intrinsic parameters may include stress state and pore pressure, size and density of fissures and fractures, fault orientation and frictional strength, steady-state coefficient of friction, rock permeability, and compressibility, as well as other geomechanical parameters such as brittleness and deformation properties (
Yaghoubi et al. 2022). However, substantial levels of inherent uncertainty are associated with the value of each intrinsic parameter. In fluid injection projects (well field brine disposal, for example), understanding the uncertainty associated with any intrinsic parameter can assist the user in making better decisions concerning user-controlled parameters such as injection pressure.
Herein, using borehole geophysics data and earthquake focal mechanism reports from nearby locations, we establish a constraint for the region’s state of stress in two target formations, the Leduc and Granite Wash Formations. We have conducted a geomechanical characterization of both formations and identified a range of uncertainties with respect to stress tensors, pore pressure, and fault properties. We then employ a probabilistic approach to determine the potential fault slip tendency for fluid injection in the formations, incorporating the uncertainty distributions associated with Mohr–Coulomb strength parameters (friction angle, normal and shear stress). A benefit of this approach is the generation of probabilistic fault stability maps that may provide a basis for future fluid injection projects in the region, such as carbon sequestration projects or well field fluid disposal.
2. Geological setting
Because of extensive drilling and public data availability, the stratigraphic columns in northwest Alberta are well established, but data are still sparse for deep formations such as the Granite Wash (
Porter et al. 1982;
Glass 1990;
Dec et al. 1996). From the oldest to youngest in northwest Alberta, the Devonian stratigraphy is made up of the Lower Devonian strata, the Elk Point Group, the Beaverhill Lake Group, and the Woodbend Group (
Porter et al. 1982;
Glass 1990;
Dec et al. 1996).
Figure 2 shows the stratigraphic column near the ABNo1 geothermal project. Based on the geothermal gradient, lithology, and hydrogeologic properties,
Banks (2016) and
Banks and Harris (2018) reported that two formations in Grande Prairie County and the MDGV region could be suitable for geothermal energy extraction: the Leduc and Granite Wash Formations. Formation selection is determined by a variety of factors, including formation temperature gradient and depth, as well as the rock physics characteristics of the formation, such as porosity and permeability (
Banks and Harris 2018). This work was corroborated by the studies for ABNo1 (
Hickson et al. 2020).
The Leduc Formation (Woodbend Group), is approximately 175–300 m thick in the Grande Prairie area, contains dominantly limestone and dolomitized limestone regions (
Banks and Harris 2018). The formation is highly porous because of widespread dolomitization during diagenesis (
Amthor et al. 1994). The Granite Wash Formation is another target geothermal fluids zone and consists of sandstone and conglomeratic sandstone of 100–200 m thickness lying directly on the Precambrian crystalline basement. This formation is also known for its high porosity and permeability (
Banks 2016;
Banks and Harris 2018). The Leduc and Grant Wash Formations lie at depths of approximately 3500 and 3800 m, respectively, in the study area.
Over 2920 wells have been drilled in Grande Prairie County (
Fig. 1) to extract fossil fuel, primarily from the Montney and Duvernay Formations (geoSCOUT™). Geomechanical parameters are available only for the Montney and Duvernay Formations, although drilled wells provide information on depth, thickness, lithology, and rock physics properties of subsurface layers such as the Leduc and Granite Wash Formations. The Duvernay and Montney Formations are “tight” shales that require HF to release hydrocarbons. HF is a process used to develop oil and gas wells using high-pressure injection of water, sand, and chemicals into rock formations.
HF is rarely used in conventional geothermal operations where steam or extremely hot fluids are present, but in the case of “hot dry rock” (Enhanced or Engineered Geothermal Systems, “EGS”), operations this technique is used to create a fractured reservoir from which to extract the heat. Basel, Switzerland, is an example of an HF EGS system. In conventional geothermal operations, injection takes place directly, at lithostatic pressure, and no proppants or other chemicals are added other than scale inhibitors which might be necessary. Nonetheless, in geothermal projects similar to the ABNo1 project, limited HF may be used to increase the flow capacity of individual wells upon analysis of initial pumping tests.
In this region (
Fig. 3), only eight wells have been drilled through to the Granite Wash Formation and an additional 28 have reached the top of the Leduc Formation.
Figure 5 shows the location of wells with the color denoting the top (depth) of the Leduc (
Fig. 3a) and Granite Wash Formations (
Fig. 3b). The depth of both formations increases from northeast to southwest as shown in the figure. For geothermal projects to be commercially viable,
Huang et al. (2021b) proposed that the target formation should be no deeper than 4500 m and its temperature should be no lower than 120 °C. Thus, the southwest area of Grande Prairie County could be a promising site for injection and production wells (
Champollion et al. 2021;
Hickson et al. 2021;
Huang et al. 2021b).
3. Seismicity in the region
A positive attribute of the Grande Prairie area is that despite the large number of HF operations carried out as part of oil and gas operations (mostly in the Duvernay and Montney Formations, stratigraphically above the Leduc and Granite Wash formations), no triggered seismic activity has been recorded. The number of seismic stations in a region will affect the level of seismic activity observed/reported in that region. The triangles in figures show the location of regional seismic stations. One of these stations has been operating since 2009, 30 km west of Grande Prairie and north of the Wapiti River (
Stern et al. 2011). In addition, there are two other seismic stations within 100 km of the area of interest (
Fig. 1). The seismic network in this area is sufficiently dense that it is certain that no significant seismic activity M ≥ 3 has resulted from fluid injection into different geological formations in the Grande Prairie area. However, in the region south and east of Grande Prairie, there are two seismogenic regions: Fox Creek and Musreau Lake (
Fig. 1), which, before the fluid injection associated with oil and gas development projects, were both quiescent.
There has been a noticeable increase in seismicity rates in Fox Creek since December 2013. More than 200 earthquakes of magnitude >2.5 have occurred in this area in association with HF operations in the Duvernay Formation, including Mw 4.1 on 12 January 2016, and Mw 3.9 on 13 June 2016 (
Bao and Eaton 2016;
Schultz et al. 2017). Most earthquakes in the Fox Creek area occurred during HF treatments and were spatially and temporally restricted to the region around the horizontal wells. Seismicity around Musreau Lake resulted from oil field wastewater injection into the Winterburn Formation at a depth of 3–4 km, the largest event being M
L 3.94 on 25 December 2019. Like Fox Creek, seismicity in Musreau Lake was spatially and temporally confined to the injection site activity. A point that should be emphasized is that induced seismicity near Fox Creek is due to HF occurring in parts of the Duvernay Formation that lay near critically stressed faults unknown before the operation began (
Fig. 1). This implies that there may be other, as-yet-unrecognized, critically stressed faults which are also susceptible to HF-induced earthquakes.
From 706 multistage HF wells, 9 × 10
6 m
3 cumulative fluid volume (geoSCOUT™) was injected into both the Duvernay and Montney Formations in the Grande Prairie region (
Fig. 4). That is important to highlight because, in a multistage HF well in Fox Creek, fluid injections of 2.04 × 10
4 m
3 triggered 115 earthquakes with a moment magnitude from 1 to 2.9 (
Yaghoubi et al. 2020). Although both regions have a comparable level of HF activity, there is a question as to why Fox Creek is a seismogenic region whereas Grande Prairie is a quiescent area. This paper addresses this question as one of its objectives.
4. Pre-existing faults in the Grande Prairie region
When evaluating fault slip, it is necessary to know the dip and dip directions of the faults that underlie the region.
Berger et al. (2009) mapped faults around the Peace River Arch, a geological structure most prominent north of Grande Prairie using regional high-resolution aeromagnetic data. Specifically,
Berger et al. (2009) explain that the Grande Prairie tectonic state is controlled by a northwest-southeast down-to-the-basin listric fault, as well as a northeast-southwest basement-involved strike-slip fault (
Berger et al. 2009). The latter may be responsible for the induced seismicity in Musreau Lake arising from wastewater injection. Because strike-slip faults may have a small or negligible vertical movement component, they are difficult to identify on seismic surveys.
In this study, the range of dip angle of each fault is described as a uniform probability distribution with the value of 80° ± 10°. It should also be noted that no laboratory study or in situ experiment has been conducted to estimate the coefficient of friction of regional faults in the study area.
Byerlee (1978) has demonstrated, based on experimental studies, that the coefficient of friction of faults varies between 0.6 and 1.0 for different rock types. In our study, we assumed that the coefficient of friction is in the range 0.7 ± 0.1. A further important point to emphasize is that our study is based on known faults in the region. As mentioned above, injection-induced seismicity around the Fox Creek area occurred on previously unknown critically stressed faults.
5. State of stress around Grande Prairie
An integral part of a region’s stress state is the pore pressure. Formation pore pressure can be directly measured from direct wellbore tests such as Drill Stem Test, Repeat Formation Test, or predicted from borehole geophysics data. For this study, we determine pore pressures within two target formations, the Leduc and Granite Wash Formations, using the data sets provided by geoSCOUT™ data services. Of 856 direct pore pressure measurements in the Leduc Formation in the WCSB, six wells including 30 tests are in the region of study.
Figures 5a and
5b display the Leduc Formation’s pore pressure gradient and its histogram presented from wells drilled in the WSCB and around Grande Prairie, respectively (
Figs. 5c and
5d). Similar to the Leduc Formation,
Fig. 6 provides pore pressure gradient variations from wells that have available direct measurements in the Granite Wash Formation. Both formations, as can be seen from the figures, are almost at hydrostatic pressure. In our study, we use a pore pressure (Pp) gradient of 9 ± 2 MPa/km for the Leduc and Granite Wash Formations. Note that pore pressure and other stress parameters are presented as gradients in this section since the depth of the target formations differs at different locations (
Fig. 3).
Since the late 1970s, extensive studies have been conducted on the principal stress orientations in British Columbia and Alberta (
Dusseault 1977;
Gough and Bell 1981;
Bell and McCallum 1990;
Bell and Bachu 2004;
Bell and Grasby 2012;
Haug and Bell 2016). Various methods have been used to determine the maximum horizontal compressive (S
Hmax) orientations in the region, including borehole failures (borehole breakouts and tensile-induced fractures) and earthquake focal mechanisms. The 2018 edition of the World Stress Map (WSM) databases include compilations of S
Hmax and relative stresses (
Heidbach et al. 2018). The blue arrows in
Fig. 1 represent the azimuth of the S
Hmax within the region and are based on borehole breakouts and tensile-induced fracture data from the WSM. We have assigned a mean of 45° and a standard division of 5° to S
Hmax azimuth values in all areas.
The vertical stress (
Sv) can be calculated by integrating the density logs from the surface to the depth of interest. Most drilled wells in the region have density logs that can be used to compute
Sv. The vertical stress component in stress tensors is less uncertain due to the availability of these density logs. The Western Canada Sedimentary Basin has been the subject of several studies that have investigated vertical stress variations, which indicate vertical stress gradients ranging from 24.6 to 25.5 MPa/km at depths greater than 2000 m (
Bell and Gough 1979;
Bell et al. 1990;
Bell and Grasby 2012;
Haug and Bell 2016;
Shen et al. 2018). In our study, we consider an
Sv range of between 24 and 26 MPa/km.
There are 706 multistage HF wells in the Grande Prairie region that have been subjected to Diagnostic Fracture Injection Testing (DFIT) or mini-frac, which can determine the minimum in situ stress magnitude (S
hmin). Only one DFIT well has been completed in the Granite Wash Formation, indicating a fracture closure pressure of 16.55 MPa/km. The closure pressure is approximately equal to S
hmin. The majority of the HF wells were completed in the Montney and Duvernay Formations in the Grande Prairie region. There are DFIT tests available in the Watt Mountain and Muskeg Formations, both of which are situated between the Leduc and Granite Wash Formations, but available data are laterally remote from the Grande Prairie area (
Fig. 7). In the Muskeg Formation, 31 measurements of closure pressure indicate an S
hmin gradient between 17 and 23 MPa/km. With 28 reported DFITs, the Watt Mountain Formation shows an S
hmin gradient of 16–25 MPa/km.
Weides et al. (2014) use a range of 13.6–19.7 MPa/km for S
hmin as input for the likelihood of fault slip due to fault injection in the Granite Wash Formation around the town of Peace River. In our study, we consider an S
hmin range of between 16 and 24 MPa/km in both the Leduc and Granite Wash Formations.
When determining a strike-slip (or thrust) stress state tensor, the most difficult parameter to determine is the maximum principal stress magnitude (S
Hmax). However, borehole failure data, together with knowledge of horizontal and vertical stresses, can be utilized to limit the range of S
Hmax magnitude (
Yaghoubi et al. 2021). Additionally, earthquake focal mechanisms provide helpful information on the relative stress magnitude as well as the maximum principal stress. To constrain the maximum principal stress magnitudes, we used injection-induced earthquake focal mechanisms recorded around the Grande Prairie region. The data set includes 11 HF-induced earthquakes around Fox Creek and 39 wastewater-induced earthquakes near Musreau Lake, Alberta (
Li et al. 2021).
Using the inversion of the focal mechanism, the Angelier’s shape parameter can be determined by

(
Angelier 1979), in which
S1 represents the maximum principal stress magnitude,
S2 represents the intermediate principal stress magnitude, and
S3 represents the minimum principal stress magnitude.
Simpson (1997) generalized the
φ values to determine the relative stress magnitudes in each stress regime by the equation:
Aφ = (
n + 0.5) + (−1)
n (
φ − 0.5) where
n = 0, 1, 2, for normal, strike-slip and reverse faulting respectively.
Aφ ranges continuously from 0 to 1 for normal faults, 1 to 2 for strike-slip faults, and 2 to 3 for reverse faults (
Hurd and Zoback 2012;
Yaghoubi et al. 2021). Simpson’s approach was applied to the 50 compiled focal mechanisms, showing that a strike-slip fault system dominates the area, with an average Anderson fault parameter of
Aφ ≈ 1.19 around Musreau Lake and
Aφ ≈ 1.5 around Fox Creek. In this study, we consider
Aφ ≈ 1.2–1.5 for the slip tendency of faults located around Grande Prairie. That means the vertical stress is closer to the minimum horizontal stress than to the maximum horizontal stress.
6. Assessment of fault-slip potential
In accordance with Coulomb faulting theory, fault slip depends on the relative stress magnitudes, the angle between the principal stress directions and the fault plane, and the coefficient of friction
μ. The slip tendency in a pre-existing cohesionless fault can be calculated in terms of the effective normal stress across the fault,
σn, to the shear stress along the fault,
τ. The likelihood for fault plane slip increases when the resolved shear stress, equals or approaches the frictional resistance of the fault surface; the fault is then referred to as being “critically stressed”. Deterministic fault-slip tendencies can be defined as the ratio between the effective normal stress and the shear stress on a potential sliding surface (
τ/
σn ≥
μ).
Figure 8 illustrates the lower hemisphere stereonet for a strike-slip fault regime with an average S
Hmax orientation of N45°E and a hydrostatic pore pressure. A pole perpendicular to a fault plane corresponds to each location within the stereonet graph. The red color indicates faults that are prone to slipping or that require less pressure to slip. Thus, nearly vertical faults striking at azimuths of 75° and 15° are critically stressed.
Deterministic analyses consider only one analysis as finite and, therefore, significantly underestimate potential risks (
Fig. 8). Based on the evidence provided in the previous section, each geomechanics parameter is associated with some level of uncertainty. A probabilistic analysis, on the other hand, examines slip tendencies by considering uncertainties inherent to each input variable, such as stress magnitudes and directions, fault dip directions, angles, and friction coefficients (
Jones and Hillis 2003;
Wang et al. 2010;
Walsh and Zoback 2016). Various input variables can be assigned as random samples with specific statistical parameters in a Mohr–Coulomb shear failure assessment. The probabilities of failure
Pf can be described as follows:
Pf =
P (
τ −
μσn ≤ 0).
Probabilistic slip tendency analysis is, therefore, more comprehensive and more appropriate for assessing risk in a variety of scenarios. This study examined the slip tendency of faults in injection formations using a Monte Carlo simulation for each fault segment. We evaluated 5000 scenarios within each fault segment as a function of pore pressure perturbation using the Mohr–Coulomb faulting theory. A sample size of 5000 was determined based on the sensitivity analysis of slip probability (with two-digit precision) versus the number of realizations in one segment fault. Uncertainties associated with intrinsic subsurface parameters, such as the state of stress, pore pressure, fault/fracture orientation, and frictional strength, are involved in the analysis. Considering that both the Leduc and Granite Wash Formations exhibit relatively similar geomechanical behavior, one analysis has been conducted to examine fault slip potential in both formations. We present the statistics that were used as inputs to the Monte Carlo simulation for fault slip tendency in the Grande Prairie region in
Fig. 9.
Some points should be noted regarding
Fig. 9. First, pore pressure and principal stress magnitudes in the target formations are assumed to follow a Gaussian distribution. Second, we assume that both the Leduc and Granite Wash Formations exhibit similar geomechanical properties. Finally, as mentioned above, since the depth of both formations varies at different locations, pore pressure and principal stress magnitudes are expressed as gradients (MPa/km). In the fault tendency analysis, however, a depth of 4000 m was applied as the injection depth.
The result is the cumulative distribution function of the critical pore pressure on each fault derived from the local tectonic stress state and the Mohr–Coulomb shear parameter analyses (
Fig. 10). As shown in the Mohr–Coulomb diagram, the mean of the principal stress magnitudes is accompanied by the associated uncertainty (green error bar). We calculated the shear stress and normal stress acting on each fault plane based on the state of stress and fault dip direction and angle. We then calculated the probability of slip for each fault segment in response to injection pressures up to 100 MPa (
Fig. 10b). The vertical black straight line in the figure indicates an injection pressure of +2 MPa greater than the pore pressure (Δ
P (
Pinj − Pp) = 2 MPa). In
Fig. 11, faults are colored according to their slip probability when injection pressure is greater than mean formation pore pressure by 2 MPa. According to the results, for this case (Δ
P = 2 MPa), there is a low probability of slippage of the faults as a result of fluid injection. The northeast-southwest basement-involving strike-slip fault in the study region has a higher probability of slip. The current result is consistent with a low level of seismicity in the region. As pore pressure increases, the probability of each fault slipping increases. For instance, the red vertical line in
Fig. 10b illustrates an injection pressure of 45 MPa in which the most probable fault has a probability slip of 50%. An increase in injection pressure to 60 MPa is certain to trigger a strike-slip fault in the northeast-southwest basement.
7. Discussion
The results of our study are consistent with the lower level of induced seismic activity observed around Grande Prairie when compared with other areas in Alberta and British Columbia that have been affected by HF (
Visser et al. 2020). This theoretical calculation is corroborated by numerous multistage hydraulic fractures having been conducted in this area that have not triggered seismic events. However, Fox Creek, about 200 km southeast of Grande Prairie, experienced some large HF-induced earthquakes (
Schultz et al. 2017;
Yaghoubi et al. 2020). HF seismicity occurred in Fox Creek as a result of multistage HF operations in the Duvernay Formation (
Schultz et al. 2017). The Duvernay Formation differs significantly from other formations in the region in terms of pore pressure gradients, one of its primary geomechanics properties. Pore pressure measurements indicate a gradient of 18 MPa/km in the Duvernay Formation around Fox Creek, which is twice as great as that in the Leduc and Granite Wash Formations in the ABNo1 study area.
Another question that might be raised is about seismicity in Musreau Lake, located less than 100 km south of Grande Prairie. The M
L 3.9 earthquake occurred in Musreau Lake as a result of
long-term wastewater being disposed into the low-pore pressure Winterburn Group. First, the Musreau Lake area is located closest to a northeast-southwest strike-slip fault, and our results indicate that the fault is more likely to slip than other faults in the region, but the likelihood is as low as 20%. Second, it is suggested by
Yu et al. (2021) that an aseismic loading slip mechanism is a triggering mechanism for the earthquake swarm around Musreau Lake. As proposed by
Li et al. (2021), long-term fluid injection in the Musreau Lake region leads to slip of faults striking in an unfavorable orientation along the northwest-southeast direction. The sequence of events at Musreau Lake may be a sign that stress accumulation has triggered a hidden fault striking in an unfavorable orientation in the region. If this is the case, long-term fluid injection in the Grande Prairie region might also cause slow-slip injection-induced seismicity on northwest-southeast down-to-the-basin listric faults.
As pore pressure increases, each fault is more likely to slip.
Figure 12 illustrates the impact of injection pressure (pore pressure) on fault slip in the study region. As can be seen, the northwest-southeast fault patches require an injection pressure of 80 MPa to slip, which is twice the current pore pressure in the Leduc and Granite Wash Formations. Another possible explanation for Grande Prairie’s quiescence is that pre-existing local faults are not in critically stressed orientations in the Grande Prairie region compared to faults in the Duvernay Formation in the Fox Creek area.
8. Conclusions
We developed and presented a method for calculating the likelihood of fault slip caused by fluid injection in a section of the WCSB. The method can be applied generally, wherever sufficient information is available to evaluate parametric uncertainty, and where faults can be identified from data sources. The major characteristics of the method are:
•
data-driven parametric uncertainty in the Mohr–Coulomb slip criteria is included;
•
data-driven parametric uncertainty in stress state, including pore pressure, is included; and
•
we evaluate various injection scenarios to generally assess the probability of slip.
This approach is applied to known mapped faults in the region that surrounds the ABNo1 geothermal site, a proposed geothermal energy development project south of Grande Prairie. Major analysis characteristics and results are:
It is necessary to have mapped faults to use this approach, and this may be a challenge in strike-slip and reverse fault domains. At Fox Creek and Musreau Lake, more distant from the ABNo1 site, events were triggered on faults that were previously unmapped, demonstrating the importance of fault identification prior to a project.