Alberta No. 1 (ABNo1) is a geothermal project targeting deep carbonate, conglomerates, and sandstone formations in a potential production and injection zone for geothermal energy exploitation within the Municipal District of Greenview south of Grande Prairie, Alberta, Canada. In geothermal systems without a steam fraction (typically systems under 170 °C), rapid widespread pore pressure changes and slow temperature changes have led to increased deviatoric stresses, resulting in induced seismicity. A concern for the ABNo1 Geothermal Project is that anthropogenic seismicity from oil, gas, and well field fluid injection has created felt events in Alberta. Thus, at the beginning of this type of project, it is prudent to review the potential for induced seismicity. In this study, a geomechanical study of the Leduc and Granite Wash Formations, two potential geothermal fluid exploitation zones, has been undertaken based on borehole geophysics and regional injection-induced earthquake data. Determining subsurface properties such as state of stress, pore pressure, and fault properties, however, poses uncertainties in the absence of actual data from the target formations. Geomechanical analysis results (with associated uncertainties) are used to assess the potential for injection-induced earthquakes. A Monte Carlo probability analysis is employed to estimate the likelihood of slippage of the known faults close to the ABNo1 Geothermal Project. A cumulative distribution function of the critical pore pressure on each fault is derived from the local tectonic stress state and Mohr–Coulomb shear parameter analyses. The resultant probabilistic fault stability maps can serve as a baseline for future fluid injection projects in the region including wastewater disposal, hydraulic fracture stimulation, CO2 sequestration, as well as geothermal energy extraction.

Renewable baseload geothermal energy can contribute significantly to meeting energy demand and decarbonizing the world’s energy sources (Frick et al. 2010; Amponsah et al. 2014). In the last few decades, geothermal energy has increasingly gained prominence as a means of reducing greenhouse gas emissions (Frick et al. 2010; Amponsah et al. 2014). In Canada, geothermal energy can help to attain net-zero emission goals of provinces with suitable geothermal assets in the Western Canadian Sedimentary Basin (WCSB) (Downey et al. 2021). In the central regions of the WCSB, Alberta has been identified as a significant source of underground geothermal energy, with an average geothermal gradient of 25–35 °C/km (Grasby et al. 2011; Hofmann et al. 2014; Hickson et al. 2020, 2021; Champollion et al. 2021; Huang et al. 2021a). The WCSB Alberta portion has been subjected to various studies to identify a suitable geothermal project site. It has been determined that the Municipal District of Greenview (MDGV) is likely to have economically viable deep geothermal resources, particularly south of Grande Prairie (Banks 2016; Banks and Harris 2018; Hickson and Colombina 2021; Huang et al. 2021a). The MDGV is bordered by the Wapiti River on the north, just south of the City of Grande Prairie.

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 km2 (Banks 2016). According to Hickson (2022), ABNo1 is expected to produce 80 000 MW·h of electrical power per year, offset 96 000 tCO2/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 ML 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 ML3.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.

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).

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 ML 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 × 106 m3 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 × 104 m3 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.

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.

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 (SHmax) 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 SHmax and relative stresses (Heidbach et al. 2018). The blue arrows in Fig. 1 represent the azimuth of the SHmax 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 SHmax 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 (Shmin). 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 Shmin. 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 Shmin gradient between 17 and 23 MPa/km. With 28 reported DFITs, the Watt Mountain Formation shows an Shmin gradient of 16–25 MPa/km. Weides et al. (2014) use a range of 13.6–19.7 MPa/km for Shmin 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 Shmin 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 (SHmax). However, borehole failure data, together with knowledge of horizontal and vertical stresses, can be utilized to limit the range of SHmax 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 forumlaφ= S2 - S3S1 - S3 (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.

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 SHmax 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.

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 ML 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.

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.

The authors thank geoLOGIC systems ltd. for their contribution of geoSCOUT™ data and software used in this study.

Data regarding formation depths, pore pressure, and minimum principal stress are available from geoLOGIC™ systems. The sources cited in this study provide access to other data sets that support the findings of this study.

A.Y. designed the project, collected data, developed codes, carried out the analysis, and wrote the manuscript with guidance and support from C.H, M.D. and Y.L. C.H, M.D. and Y.L. supervised the project and contributed to the interpretation of the results and provided critical feedback.

The ABNo1 project is partially funded through a grant from Natural Resources Canada’s Emerging Renewable Power Program. The MITACS program ( and Terrapin Geothermics have supported this research through leveraging the research support provided by the ABNo1 project to the first author.

This work is licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

Competing Interests

The authors declare there are no competing interests.