CC-BY 4.0
1. Introduction
Canadian remote northern regions lack a reliable energy source. The majority of communities in such arctic to subarctic regions heavily rely on diesel for electricity and space heating [Grasby et al. 2013]. As a matter of fact, in Nunavik (northern Quebec, Canada), 25 million liters of diesel are consumed annually for electricity generation and 28 million for heat production in the 14 Inuit communities of this region [data for 2009; KRG and Makivik Corporation, 2010]. Beyond the environmental issues that such energetic framework entails, it also impacts economic progress. The development of deep geothermal resources in such areas may be a solution for such energetic problems. In fact, a previous regional geothermal potential study undertaken by Grasby et al. [2012] indicated the areas of interest and possible technologies for deep geothermal development (Figure 1). The generated map suggests the Canadian Shield as a possible target for future Enhanced Geothermal Systems (EGS) development.
Distribution of geothermal potential in Canada based on end use and location of the electrical grid and remote communities. Adapted from Arriaga et al. [2017] and Grasby et al. [2012].
Thus, studying the feasibility of engineered geothermal energy systems, or EGS, in remote northern communities settled on the Canadian Shield is of utmost interest. However, prior to design an engineered geothermal energy system, characterizing the thermal, hydraulic and mechanical structure of the target area is essential [e.g., Nakatsuka 1999]. Miranda et al. [2020a] carried out a comprehensive analysis of the thermophysical properties of the outcropping rocks in the study area. Miranda et al. [2021a] evaluated the heat flux in the study area by developing a numerical approach to infer this parameter from shallow boreholes (80 m deep). Their results suggest heat flux in the range 31.8–69.4 mW⋅m−2 depending on paleoclimate and thermophysical properties conditions. Miranda et al. [2020b] used these findings to estimate the subsurface temperature distribution and evaluate the theoretical and technical potential following the Beardsmore Protocol [Beardsmore et al. 2010]. A temperature of 79–88 °C at 4 km depth and a potential heat output with a probability of 98% to fulfil Kuujjuaq’s (Nunavik, Canada) estimated annual average heating demand of 37 GWh was evaluated. Miranda et al. [2018] carried out a first-order analysis of the geometrical and topological (i.e., connectivity) properties of the fractures sampled in Kuujjuaq. This analysis was further improved by Miranda [2021] who provide a comprehensive statistical analysis of the fractures sampled. This research work also presents several stochastic fracture networks built upon that statistical analysis and their influence on the fluid flow and heat transfer. Finally, in Miranda et al. [2021b], the technical feasibility and economic viability of an EGS in Kuujjuaq were assessed and hydraulically stimulated geothermal reservoirs of potential commercial interest were simulated. However, these results are speculative and subject to uncertainty due to the poor knowledge, among others, of the in situ stresses.
Thus, the present manuscript aims at tackling the stress regime in the community of interest. In fact, the state of stress is a fundamental geomechanical issue to design engineered geothermal energy systems [Brown et al. 2012; Evans et al. 1999; Ghassemi 2012; Richards et al. 1994]. The in-situ stress components not only directly control the pressure required for reservoir stimulation but also influence the lateral extent and orientation of the stimulated rock volume [Evans et al. 1999; Ghassemi 2012]. This man-made geothermal reservoir tends to be an elliptical shape, with the major axis parallel to the maximum principal stress and normal to the least principal stress, hence reflecting the tendency of the reservoir to grow in the direction of the least energy configuration [e.g., Brown et al. 2012; Evans et al. 1999; Hubbert and Willis 1957]. Beyond the orientation of the principal stresses, their absolute and relative magnitudes are also key components of a comprehensive geomechanical model, together with the knowledge of the in-situ fluid pressure [e.g., Jaeger et al. 2007; Zoback 2007]. The in-situ stress regime largely determines the in-situ fracture shear strength and, thus, the necessary additional fluid pressure required to activate shear slip in the optimally oriented fractures (i.e., the reduction of effective normal stress along fractures that support shear stress). Moreover, if a fracture is critically stressed (i.e., the ratio between the applied shear and the effective normal effective stress approaches the friction angle—Amonton’s law) and, thus, is on the verge of failing, then less fluid pressure is required to reactivate the structure [e.g., Evans et al. 1999; Jaeger et al. 2007; Zoback 2007]. A linear Coulomb friction law with a friction coefficient ranging between 0.6 and 1.0 is commonly used to provide a first-order approximation for the upper limit of the shear strength necessary to reactivate structures and initiate shearing [e.g., Evans et al. 1999; Jaeger et al. 2007; Zoback 2007]. However, the friction coefficient depends on rock type and fracture toughness.
A number of techniques have been developed over the years to determine the in-situ stress. These include small volume hydraulic fracture tests, analysis of breakouts and drilling-induced fractures and hydraulic testing of pre-existing fractures [e.g., Amadei and Stephansson 1997]. However, when none of these techniques can be applied due to difficulties such as the lack of deep boreholes (as the case of the study area) or due to the high costs of in situ stress measuring tests, a different approach for an a priori estimation of the stress regime in a specific area is needed. According to Amadei and Stephansson [1997], the in-situ stress field can be estimated from stress-depth relationships, observations of past stress measurements made in the region of interest or by extrapolation from regions with similar geologic and tectonic settings. Topography, geology, rock fabric, rock loading history, first motion analysis of earthquakes, occurrence of stress release phenomena, breakouts in boreholes, tunnels and shafts, rock bursts, and the presence of stratification, heterogeneities or geological structures are also useful to derive stress information. Another key source for a first-order assessment of stress regime is the World Stress Map [Heidbach et al. 2010; Zoback et al. 1989]. In this work, the stress regime is evaluated based on a theoretical approach that uses the World Stress Map and other literature data and stress-depth relationships providing a first-order approximation. The orientation of the stress components was inferred based on the World Stress Map and stress compilations and their magnitude was estimated based on empirical correlations.
1.1. Objectives
The objectives of this work are twofold. First, it aims at defining an a priori estimate of the stress regime in Kuujjuaq. The World Stress Map and other literature data are used to estimate the orientation of the stress components and empirical correlations are applied to estimate the magnitude of the principal stress components, in-situ fluid pressure, theoretical boundaries for the principal stresses, and effective stimulation pressure. The estimated magnitude was posteriorly compared with available stress data for the Canadian Shield.
Then, a first-order estimation of the fracture sets optimally oriented to slip is made through Mohr–Coulomb friction and slip tendency analyses. This was done taking into account the a priori stress model previously defined.
The findings of this work are helpful to discuss the best EGS design for Kuujjuaq, keeping in mind the lessons learned from historic sites where hydraulic stimulation has been applied. This work aims at highlighting how the poor knowledge of the stress regime and associated uncertainty has an important impact in designing and developing EGS.
The community of Kuujjuaq was used in this work as a case study and an example for the remaining settlements located in a similar geological context, facing critical energy issues and having geothermal data gaps. Although the potential error associated with the use of literature-based and empirical data is significant, a first-order characterization of the mechanical conditions of the subsurface is crucial to identify the questions that will have to be resolved by future work if possibilities suggested here are sufficient to justify further geothermal exploration in such remote northern regions.
1.2. Geographical, geological and structural setting
The study area is the community of Kuujjuaq, the regional capital of Nunavik (Figure 2a). This Inuit community is home to almost 3000 inhabitants, the majority Indigenous. Kuujjuaq is located above the 55th parallel and, over the course of the year, it experiences temperatures varying from −28 °C to 18 °C. The annual average surface air temperature is −6 °C and the annual average ground surface temperature is −1 °C [ECCC, 2019].
(a) Geographical location of Kuujjuaq and remaining remote communities and of the deep boreholes in Nunavik; (b) geological setting of Kuujjuaq and surrounding area. LP—Lac Pingiajjulik fault, LG—Lac Gabriel fault.
Kuujjuaq is located in the periphery of northern Labrador Trough in southeastern Churchill Province (SECP) of the Canadian Shield. This geological province is an orogenic belt, oriented NNW–SSE with both flanks (Torngat Orogen and Labrador Trough) recording transpressional development associated with the oblique collisions that led to its assemblage [e.g., Wardle et al. 2002]. The collision between the Core Zone (central part of the SECP) and the Superior Province occurred about 1.82–1.77 Ga before present. The Core Zone is a cratonic fragment mostly consisting of Archean tonalitic to granitic gneiss and granitoid rocks and is characterized by a pervasive E-dipping fabric related to westerly thrusting [e.g., Wardle et al. 2002]. A regional structural grain NW–SE to N–S is observed and believed to be associated with the Core Zone-Superior collision [Simard et al. 2013].
The crustal thickness in northern Quebec was estimated by Vervaet and Darbyshire [2016] to vary within 33 to 49 km with Moho at a depth of about 37 km beneath the Kuujjuaq terrane [Hall et al. 2002; Hammer et al. 2010; Telmat et al. 1999]. Seismic profiles acquired in the east coast of Ungava Bay in the scope of the Lithoprobe project suggest P-wave velocities ranging between 5.9 and 6.2 km⋅s−1 for the upper crustal levels (from 0 to 10 km depth), between 6.2 and 6.5 km⋅s−1 for the mid-crust (from 10 to about 17 km depth) and between 6.6 and 7.0 km⋅s−1 for the lower crust [from around 17 to 37 km depth; Hall et al. 2002; Hammer et al. 2010].
The effective elastic thickness of the Canadian Shield estimated by the maximum entropy method suggests a value ranging between 70 and 85 km within Kuujjuaq area [Audet and Mareschal 2004]. Taking into consideration the approximated effective elastic thickness and crustal thickness, the ratio between these parameters is 1.9–2.3, which suggests a strong “dried jelly sandwich” rheology where the lower crust is strong and both crust and mantle are mechanically coupled [Burov 2011].
The study area is characterized by outcrops of the False Suite [migmatized paragneiss and migmatized garnet paragneiss; SIGÉOM, 2019] and Kaslac Complex [amphibolite diorite and quartz diorite and gabbro, gabbronorite and clinopyroxene; SIGÉOM, 2019]. Smaller outcrops belonging to the Ralleau Suite [amphibolized gabbro and diorite; SIGÉOM, 2019], Aveneau Suite [white tonalite and granite; SIGÉOM, 2019], Dancelou Suite [massive pink granite and massive pegmatitic granite; SIGÉOM, 2019] and Falcoz Swarm [subophilitic gabbro; SIGÉOM, 2019] are also present (Figure 2b). A detailed description of these lithologies can be found in Miranda et al. [2020a] and SIGÉOM [2019]. The main structure crossing the community, the Lac Pingiajjulik fault (Figure 2b), was not observed in the study area, but literature suggests a thrust fault striking NW–SE with a strike-slip component indicating transpressional regime [Simard et al. 2013].
2. Stress regime in the Canadian Shield
The World Stress Map [Heidbach et al. 2010, 2019] is a global database of contemporary tectonic stress of the Earth’s crust and a useful tool to infer the orientation of the maximum horizontal stress. However, the World Stress Map does not have information on the study area. The lack of deep exploratory boreholes associated with the absence of earthquake data with magnitudes greater than 3 contributes for this important gap (Figure 3).
(a) Stress map [Heidbach et al. 2019] and (b) seismic hazard map [NRCan 2018] of Quebec province. The arrows in (a) indicate the regional trend of the contemporary stress field and the orientation of the maximum horizontal compression based on Adams [1989].
No deep borehole is available nearby Kuujjuaq to carry out stress measurements, being the nearest borehole placed at a distance of about 400 km (Figure 2a). Moreover, the Coulon borehole was drilled in the Superior Province and the Raglan and Asbestos Hill, although drilled in the Churchill Province, experienced a different geodynamic assembly history compared to the part of the Churchill Province where Kuujjuaq is located [e.g., Whitmeyer and Karlstrom 2007]. Two other deep boreholes, the Voisey Bay’s (Newfoundland and Labrador) and Nielsen Island’s (Nunavut), lie at a distance of more than 400 km and were drilled in the Nain and Superior provinces, respectively. Nevertheless, a compilation of available data was undertaken to overcome this data gap, assuming that Canadian Shield intraplate region has relatively consistent maximum horizontal stress orientations. Several authors carried out stress measurements in mining sites on the Canadian Shield and compiled the available data [Adams 1989; Andrieux et al. 2009; Arjang 1991, 1998; Herget 1982, 1987, 1993; Young and Maloney 2015]. It is important to highlight that these mining sites do not have the same geologic neither tectonic setting as the study area. Their observations suggest a regional compression trend NE–SW for the maximum horizontal stress in the Canadian Shield [Adams 1989; Arjang 1991, 1998; Herget 1982, 1987, 1993; Young and Maloney 2015] except in Raglan mining site where the maximum horizontal stress is oriented N–S [Andrieux et al. 2009]. The available stress data suggests reverse/strike-slip faulting regime [Adams 1989; Andrieux et al. 2009; Arjang 1998; Herget 1993; Young and Maloney 2015] and some measurements revealed a vertical stress greater than the minimum principal stress component but smaller than the intermediate principal stress component [Herget 1993; Young and Maloney 2015]. The gradient of the vertical stress is found to be about 26.0 MPa⋅km−1 [Andrieux et al. 2009; Arjang 1991, 1998; Herget 1982, 1987, 1993; Kaiser et al. 2016; Young and Maloney 2015]. The gradient of the minimum horizontal stress ranges between 23.0 and 37.8 MPa⋅km−1 while the gradient for the maximum horizontal stress between 34.4 and 51.3 MPa⋅km−1 [Andrieux et al. 2009; Arjang 1991, 1998; Herget 1993; Kaiser et al. 2016; Young and Maloney 2015]. However, it is important to highlight that these gradients vary with depth, making it greatly uncertain when used for comparison purposes.
3. Methodology
3.1. Orientation of the principal stresses
The orientation of the principal stresses was estimated based on the available literature data for the Canadian Shield (Table 1). The following hypotheses were used in the analyses described below to assess the impact of the orientation of the principal stress components on the EGS development and performance. The maximum principal stress was assumed to have a trend and plunge of N248°/10° [Herget 1993], of N227°/02° [Young and Maloney 2015] and N–S/horizontal [Andrieux et al. 2009]. The intermediate principal stress was assumed to be N300–340°/0° [Herget 1993], N310°/08° [Young and Maloney 2015] and E–W/horizontal [Andrieux et al. 2009]. The minimum principal stress was assumed vertical [Andrieux et al. 2009; Herget 1993] and N270°/88° [Young and Maloney 2015].
Principal stresses in the Canadian Shield
| Principal stress | Orientation | Magnitude | Observations | Reference | |
|---|---|---|---|---|---|
| 𝜎V | — | (0.0260–0.0324)z | 0 < z < 2200 m | 𝜎V <𝜎 H,average | Herget [1982, 1987] |
| 9.86 + 0.0371z | 0 < z < 900 m | ||||
| 𝜎H,average | — | 33.41 + 0.0111z | 900 < z < 2200 m | ||
| 𝜎V | — | (0.0266 ± 0.008)z | 60 < z < 1890 m | 𝜎V <𝜎 h <𝜎 H | Arjang [1991] |
| 𝜎H,average | — | 5.91 + 0.0349z | |||
| 𝜎H | — | 8.18 + 0.0422z | |||
| 𝜎h | — | 3.64 + 0.0276z | |||
| 𝜎V | — | 0.0285z | 0 < z < 2200 m | 𝜎V <𝜎 h <𝜎 H | Herget [1993] |
| 𝜎1 | N248°/10° | 12.1 + (0.0403 ± 0.0020)z | |||
| 𝜎2 | N300–340°/0° | 6.4 + (0.0293 ± 0.0019)z | |||
| 𝜎3 | Vertical | 1.4 + (0.0225 ± 0.0015)z | |||
| 𝜎V | — | 0.0260z | 0 < z < 6000 m | 𝜎V <𝜎 h <𝜎 H | Arjang [1998] |
| 𝜎1 | NE/horizontal | 13.50 + 0.0344z | |||
| 𝜎2 | NW/sub-horizontal | 8.03 + 0.0233z | |||
| 𝜎3 | Vertical | 3.01 + 0.0180z | |||
| 𝜎1 | N–S/horizontal | 0.0513z | — | 𝜎V <𝜎 h <𝜎 H | Andrieux et al. [2009] |
| 𝜎2 | E–W/horizontal | 0.0378z | |||
| 𝜎3 | Vertical | 0.0270z | |||
| 𝜎V | — | (0.0258–0.0263)z | 12 < z < 2552 m | 𝜎V <𝜎 h <𝜎 H | Young and Maloney [2015] |
| 𝜎1 | N227°/02° | (0.040 ± 0.001)z − (9.185 ± 1.5) | |||
| 𝜎2 | N310°/08° | (0.029 ± 0.001)z + (4.617 ± 1.159) | |||
| 𝜎3 | N270°/88° | (0.021 ± 0.001)z − (0.777 ± 0.872) | |||
| 𝜎V | — | 0.021z | 0 < z < 1300 m | 𝜎V <𝜎 h <𝜎 H | Kaiser et al. [2016] |
| 𝜎1 | — | 0.012z + 42.4 | 660 < z < 1300 m | ||
| 𝜎2 | — | 0.013z + 24.1 | |||
| 𝜎3 | — | 0.007z + 9.7 | |||
𝜎V—vertical stress, 𝜎H—maximum horizontal stress, 𝜎h—minimum horizontal stress, 𝜎H,average—average horizontal stress (𝜎H +𝜎 h)∕2, 𝜎1—maximum principal stress, 𝜎2—intermediate principal stress, 𝜎3—minimum principal stress.
3.2. Magnitude of the principal stresses
The three principal stresses can be assumed as acting vertically and horizontally as a first approximation [Amadei and Stephansson 1997]. Bearing in mind this assumption, then the following empirical correlations can be used as a first-order approximation for the magnitude of the principal stresses and in-situ fluid pressure.
The vertical stress component can be simply estimated based on the weight of the overlying rock at depth as [e.g., Hoek and Brown 1980]:
| (1) |
The horizontal stress components were estimated based on the horizontal to vertical stress ratio as [e.g., Zhang 2017]:
| (2) |
Horizontal to vertical stress ratios inferred for the Canadian Shield
| Stress ratio coefficient | Expression | Reference |
|---|---|---|
| kmax | (357∕z) + 1.46 | Herget [1987] |
| kmin | (167∕z) + 1.10 | Herget [1987] |
| kmax | ((272 ± 8)∕z) + 1.72 | Herget [1993] |
| kmin | ((30 ± 4)∕z) + 0.86 | Herget [1993] |
| kmax | 7.44z−0.198 | Arjang [1998] |
| kmin | 2.81z−0.120 | Arjang [1998] |
kmax—maximum horizontal to vertical stress ratio, kmin—minimum horizontal to vertical stress ratio, z—depth in meters.
A transpression deformation style is favorably developed in tectonic regimes near the transition between compression (minimum principal stress is vertical and the R ratio is small) and wrench (intermediate principal stress is vertical and the R is small). The R ratio is given by the following expression assuming compression tectonic regime:
| (3) |
The in situ fluid pressure can be estimated based on the pore-fluid factor [e.g., Sibson 2004; Zoback and Townend 2001]:
| (4) |
Earth’s crust contains faults, fractures and planar discontinuities at many different scales and orientations and, thus, the stress magnitudes are limited by the frictional strength of these discontinuities [Zoback 2007]. Theoretical boundaries for the maximum and minimum principal stresses were estimated based on the frictional equilibrium theory [e.g., Zoback 2007]:
| (5) |
The maximum and minimum effective principal stresses, assuming reverse faulting regime, are calculated as:
| (6) |
An additional evaluation of the admissible stress state in Kuujjuaq was carried out through stress polygons [Zoback 2007]. Stress polygons are a simple way to illustrate the range of allowable value for horizontal principal stresses in the Earth’s crust for normal, strike-slip and reverse faulting environments [Zoback 2007]. The horizontal stresses can be anywhere in the triangles, indicating a specific stress regime (i.e., normal, reverse and strike-slip). The following formulations are needed to build the stress polygon [Taghipour et al. 2019]:
| (7) |
The stress polygons can help to illustrate the range of possible magnitudes of the minimum horizontal and maximum horizontal stresses at a particular depth for a given in situ fluid pressure and assumed coefficient of friction. The vertical line in the lower left of the polygon indicates the lowest value of the minimum horizontal stress possible in a normal faulting environment. Similarly, the horizontal line defining the top of the polygon corresponds to the value of the maximum horizontal stress at which reverse faulting would occur. The diagonal line bounding the polygon on the upper left corresponds to the value of the maximum horizontal stress at which strike-slip faulting would occur for a given minimum horizontal stress.
Monte Carlo-based sensitivity analysis was undertaken to estimate the magnitude of the principal stress components (Equations (1)–(7)). This approach allows to take into account the uncertainty of each input parameter (Table 3). Density of the geological materials was defined based on bulk density measurements from laboratory tests on outcrop samples [Miranda et al. 2020a]. Density is observed to range between the minimum value of 2465 kg⋅m−3 and the maximum value of 3108 kg⋅m−3, with a median value of 2677 kg⋅m−3. A triangular distribution was used for this parameter based on the assumption that the minimum and maximum density values obtained are least probable than the median value calculated from the dataset analyzed. A single value was assumed for the gravitational acceleration and depth. However, it is convenient to highlight that the calculations were done for depths up to 10 km, but in this work, only the results up to 5 km are shown.
Monte Carlo method input parameters and their uncertainty
| Parameter code | Parameter description | Unit | Variable type | Distribution |
|---|---|---|---|---|
| 𝜌 | Density | kg⋅m−3 | Continuous | Triangular (min, med, max) |
| g | Gravitational acceleration | m⋅s−2 | — | Single value |
| z | Depth | m | — | Single value |
| k | Stress ratio coefficient | — | Discrete continuous | Uniform (min, max) |
| 𝜎V | Vertical stress | Pa | Continuous | Triangular (min, mean, max) |
| Uniform (min, max) | ||||
| Fp−f | Pore to fluid factor | — | — | Single value |
| 𝜎H | Maximum horizontal stress | Pa | Continuous | Triangular (min, mean, max) |
| Uniform (min, max) | ||||
| 𝜇 | Friction coefficient | — | Continuous | Uniform (min, max) |
| Pp | In situ fluid pressure | Pa | Continuous | Triangular (min, mean, max) |
| Uniform (min, max) |
The stress ratio coefficients were assumed to follow a discrete uniform distribution based on the assumption that the probability of each stress ratio coefficient is a finite value and that they have an equal probability of occurrence. However, a continuous uniform distribution was additionally tested to evaluate the impact on the results.
The vertical stress in Equation (2) was assumed to follow a triangular distribution since that is the distribution shape coming from Equation (1). Nevertheless, a continuous uniform distribution was additionally assumed for this parameter to assess the impact of the probability distribution on the outcome.
The pore to fluid factor was assumed as a single value while the vertical stress was assumed to follow continuous triangular and uniform distributions (Table 3). Triangular and uniform probability distributions were used to describe the maximum horizontal and vertical principal stresses and in situ fluid pressure (Table 3). Finally, the friction coefficient was assumed to range between 0.6 and 1.0 [e.g., Zoback 2007] and described by a uniform distribution (Table 3).
The software @RISK [Palisade 2019] was used to carry out the Monte Carlo simulations. Latin Hypercube sampling and the pseudorandom number generator Marsenne Twister were used in the simulations. A total of 10,000 iterations were run per simulation to assure output stability and the initial random number seed was fixed to 1 in all the simulations carried out.
The estimated values for the principal stress components were plotted to find a trend line that best describes the estimated stress with depth. The trend lines were then compared with the literature data to assess the uncertainty of estimating stress magnitude based on this theoretical approach.
3.3. Effective stimulation pressure
Maximum pressures required for reservoir dilation can be estimated following the Hubbert–Willis classical expression for impermeable rocks Hubbert and Willis [1957]:
| (8) |
| (9) |
The critical fluid pressure was also estimated based on Monte Carlo simulations and assuming a reverse faulting regime. Both effective maximum and minimum principal stresses were assumed to follow triangular and uniform distributions (Table 3).
3.4. Mohr–Coulomb friction and slip tendency
Hypotheses tested in the Mohr–Coulomb friction and slip tendency analyses
| Principal stresses | Tested hypotheses | |||
|---|---|---|---|---|
| H1 | H2 | H3 | Effective magnitude (mean value; MPa) | |
| 𝜎1 | N248°/10° | N227°/02° | N–S | 94 |
| 𝜎2 | N340°/0° | N310°/08 | E–W | 66 |
| 𝜎3 | Vertical | N270°/88° | Vertical | 64 |
𝜎1—maximum principal stress, 𝜎2—intermediate principal stress, 𝜎3—minimum principal stress.
Fracture reactivation analysis is commonly done based on the Mohr–Coulomb criterion [e.g., Moeck et al. 2009; Peacock et al. 2021]. However, care must be taken when using Mohr-circle analyses due to a potential underestimation of the risk of fault reactiviation [van den Bogert and van Eijs 2020]. The MohrPlotter software version 3.0.0 [Allmendinger 2020] was used in this study to help determine which fracture sets are optimally oriented to slip taking into account the possibilities for the orientation and magnitude of the principal stresses (Table 4). The trend and plunge of the principal stresses were entered in the software in the geographic coordinate system format assuming a North-East-Down coordinate system for the stress tensor. A linear Coulomb friction law was assumed in this study for the failure criterion and a constant friction coefficient of 0.6 was assumed to simplify the analyzes.
The slip tendency of each fracture plane was evaluated as the ratio of shear to effective normal stress acting on the fracture surface [e.g., Morris et al. 1996; Rutqvist et al. 2007]:
| (10) |
4. Results
4.1. Magnitude of the principal stresses
Monte Carlo simulations applied to Equations (1)–(3) to solve the vertical and horizontal principal stresses and in situ fluid pressure suggest a range of possible values for these parameters at depths up to 5 km (Table 4). The Monte Carlo results also highlight that despite the different dispersion of the results (bimodal versus unimodal), the choice of the probability distribution function as a small impact on the statistics (Table 5; Figure 4). The minimum, mean and maximum values of each stress component for each distribution hypothesis are within the same range of values. For this reason, and to simplify the analysis, the hypotheses were averaged (Table 5; Figure 4).
Monte Carlo simulation results illustrating the distribution of the principal stress components for a depth of 4 km as an example. This depth was chosen since previous research suggested this to be the ideal depth to harness the geothermal resources for direct use of heat.
The sensitivity analysis undertaken also suggests that the stress ratio coefficient is the most influential parameter when estimating the maximum and minimum horizontal stresses. The Spearman correlation is about 80% for the stress ratio coefficient and 50% for the vertical stress. This suggests that increasing both the stress ratio coefficient and vertical stress leads to a greater magnitude of both maximum and minimum principal stresses (Figure 4).
Plotting the values of each principal stress and in situ fluid pressure as a function of depth results in the following relationships:
| (11) |
Comparing the estimated gradients for each of the principal stresses with values from literature (Table 1) reveals a good agreement of the data, with the stress gradient within the literature range. The stress ratio coefficient estimated in this work is also within the literature range (Figure 5). The R ratio (Equation (4); Table 4) is found to be less than 1, varying between 0.2 at 1 km depth and 0.03 at 5 km, suggesting a tectonic regime near the transition between compression and wrench, favorable for transpression deformation style.
Estimated stress components
| Depth (km) | Statistical parameters | Principal stresses | |||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 𝜎V (MPa) | 𝜎h (MPa) | 𝜎H (MPa) | Pp (MPa) | ||||||||||||||
| A | B | C | D | Average | A | B | C | D | Average | A | B | Average | |||||
| 1 | Mean | 27.0 | 30.4 | 31.0 | 29.1 | 29.7 | 30.0 | 43.6 | 44.4 | 39.8 | 40.5 | 42.1 | 10.9 | 11.0 | 10.9 | ||
| St dev | 1.3 | 4.8 | 5.2 | 3.3 | 3.7 | 3.0 | 9.5 | 10.0 | 5.7 | 6.2 | 5.7 | 0.6 | 0.8 | 0.5 | |||
| Median | 26.8 | 32.3 | 32.1 | 29.0 | 29.5 | 30.3 | 48.1 | 47.7 | 39.7 | 40.3 | 42.9 | 10.9 | 11.0 | 10.9 | |||
| [Min–Max] | [24.2–30.4] | [21.6–38.6] | [21.4–39.3] | [22.0–38.3] | [21.4–39.2] | [21.9–38.1] | [27.4–57.6] | [27.2–58.7] | [27.8–54.5] | [27.4–56.2] | [28.1–56.2] | [9.6–12.4] | [9.6–12.4] | [9.7–12.3] | |||
| 2 | Mean | 54.0 | 57.3 | 57.9 | 55.5 | 56.1 | 56.7 | 77.1 | 77.9 | 71.1 | 71.8 | 74.5 | 21.7 | 21.8 | 21.7 | ||
| St dev | 2.6 | 7.8 | 8.4 | 5.6 | 6.2 | 5.0 | 16.9 | 17.6 | 10.6 | 11.3 | 10.2 | 1.1 | 1.5 | 0.9 | |||
| Median | 53.7 | 59.9 | 59.2 | 55.4 | 55.8 | 57.0 | 85.5 | 84.4 | 71.1 | 71.3 | 75.7 | 21.7 | 21.8 | 21.7 | |||
| [Min–Max] | [48.4–60.9] | [42.4–72.1] | [42.0–72.2] | [42.8–71.5] | [42.1–72.0] | [43.5–70.7] | [48.4–100.6] | [47.8–100.8] | [49.2–98.5] | [48.0–99.6] | [49.0–97.9] | [19.2–24.4] | [19.2–24.4] | [19.4–24.3] | |||
| 3 | Mean | 81.0 | 83.7 | 84.8 | 82.0 | 83.1 | 83.4 | 109.4 | 110.8 | 102.4 | 103.7 | 106.6 | 32.7 | 32.8 | 32.6 | ||
| St dev | 3.9 | 10.6 | 11.3 | 7.8 | 8.6 | 6.8 | 23.6 | 24.4 | 15.5 | 16.4 | 14.4 | 1.5 | 2.1 | 1.3 | |||
| Median | 80.5 | 86.3 | 86.4 | 81.7 | 82.8 | 83.6 | 120.6 | 120.3 | 102.1 | 103.2 | 108.1 | 32.6 | 32.8 | 32.6 | |||
| [Min–Max] | [72.6–91.3] | [63.3–105.5] | [63.5–105.2] | [64.0–104.5] | [63.9–104.9] | [64.5–102.5] | [69.3–144.1] | [69.4–143.7] | [70.4–143.2] | [69.5–142.8] | [71.3–139.0] | [29.2–36.3] | [29.2–36.4] | [29.2–36.2] | |||
| 4 | Mean | 107.9 | 109.6 | 111.2 | 108.4 | 110.0 | 109.8 | 140.9 | 143.0 | 133.7 | 135.6 | 138.3 | 43.6 | 43.8 | 43.5 | ||
| St dev | 5.2 | 13.4 | 14.4 | 10.0 | 11.4 | 8.7 | 30.1 | 31.1 | 20.4 | 21.7 | 18.7 | 2.0 | 2.9 | 1.8 | |||
| Median | 107.4 | 111.5 | 112.2 | 108.2 | 109.5 | 109.9 | 154.0 | 154.1 | 133.3 | 135.3 | 140.0 | 43.5 | 43.8 | 43.5 | |||
| [Min–Max] | [96.8–121.8] | [84.0–139.0] | [84.2–139.3] | [84.4–137.9] | [84.4–139.0] | [85.1–136.3] | [89.9–188.6] | [90.0–189.0] | [91.2–186.3] | [90.5–188.1] | [93.1–185.3] | [38.9–48.8] | [38.8–48.8] | [39.0–48.4] | |||
| 5 | Mean | 134.9 | 135.4 | 137.0 | 134.9 | 136.5 | 135.9 | 171.9 | 174.0 | 165.0 | 166.9 | 169.5 | 54.4 | 54.6 | 54.3 | ||
| St dev | 6.6 | 16.2 | 17.4 | 12.3 | 13.8 | 10.7 | 36.2 | 37.6 | 25.4 | 26.7 | 22.7 | 2.5 | 3.6 | 2.2 | |||
| Median | 134.2 | 135.6 | 137.5 | 134.5 | 135.9 | 137.1 | 184.5 | 186.9 | 164.5 | 166.5 | 171.2 | 54.3 | 54.6 | 54.3 | |||
| [Min–Max] | [121.0–152.3] | [104.8–172.4] | [104.8–172.3] | [105.9–170.0] | [104.9–172.1] | [106.7–168.3] | [110.8–233.1] | [110.6–232.8] | [112.5–231.2] | [111.1–232.1] | [112.8–225.8] | [48.5–60.7] | [48.4–60.8] | [48.6–60.5] | |||
𝜎V—vertical stress, 𝜎H—maximum horizontal stress, 𝜎h—minimum horizontal stress, Pp—in situ fluid pressure. A—𝜎V triangular distribution and k discrete distribution, B—𝜎V uniform distribution and k discrete distribution, C—𝜎V triangular distribution and k continuous distribution, D—𝜎V uniform distribution and k continuous distribution.
Literature-based versus estimated (a) minimum and (b) maximum stress ratio coefficients as a function of depth. In (c) and (d) are the distribution of the stress ratio coefficients for a depth of 4 km as an example.
The frictional equilibrium theory criterion (Equation (5)) suggests that the estimated ratio between effective maximum principal stress and effective minimum principal stress is below 4.3, the theoretical frictional limit for a critically oriented discontinuity assuming a friction coefficient varying between 0.6 and 1.0 (Table 6; Figure 6). This suggests that the subsurface beneath Kuujjuaq may not be critically stressed, that the ratio between principal stresses is low for an effective hydraulic stimulation treatment at low fluid overpressure and that high injection pressures may be needed to reactivate the well-oriented pre-existing discontinuities.
Effective stress ratio for a depth of 4 km as an example.
Estimated effective stress ratio coefficient
| Depth (km) | Statistical parameters | |||
|---|---|---|---|---|
| A | B | Average | ||
| 1 | Mean | 1.9 | 1.9 | 1.9 |
| St dev | 0.4 | 0.6 | 0.3 | |
| Median | 1.9 | 1.9 | 1.9 | |
| [Min–Max] | [0.9–3.3] | [0.9–3.6] | [1.0–3.1] | |
| 2 | Mean | 1.6 | 1.6 | 1.6 |
| St dev | 0.3 | 0.5 | 0.3 | |
| Median | 1.6 | 1.6 | 1.6 | |
| [Min–Max] | [0.7–2.7] | [0.7–3.0] | [0.8–2.6] | |
| 3 | Mean | 1.5 | 1.5 | 1.5 |
| St dev | 0.3 | 0.4 | 0.3 | |
| Median | 1.5 | 1.5 | 1.5 | |
| [Min–Max] | [0.7–2.6] | [0.7–2.7] | [0.8–2.5] | |
| 4 | Mean | 1.5 | 1.5 | 1.5 |
| St dev | 0.3 | 0.4 | 0.3 | |
| Median | 1.5 | 1.4 | 1.5 | |
| [Min–Max] | [0.7–2.6] | [0.6–2.7] | [0.7–2.4] | |
| 5 | Mean | 1.4 | 1.4 | 1.4 |
| St dev | 0.3 | 0.4 | 0.3 | |
| Median | 1.4 | 1.4 | 1.4 | |
| [Min–Max] | [0.7–2.5] | [0.6–2.6] | [0.7–2.4] | |
—effective maximum principal stress, —effective minimum principal stress.
Additionally, the sensitivity analysis suggests that the maximum principal stress (in this case the maximum horizontal stress) is the most influential parameter with a Spearman correlation of 90%. Furthermore, the effective minimum principal stress (here assumed the vertical stress) has a negative correlation of 40%. This indicates that decreasing the vertical stress and increasing the maximum horizontal stress leads to an increase in the theoretical frictional limit (Figure 6) and ultimately places the fractures in Kuujjuaq near a critical state of stress. This may lead to more effective hydraulic stimulation treatments. The in situ fluid pressure is observed to have a negligible influence on the results with a Spearman correlation of less than 1%. The variability imposed by the probability distribution functions is once again small (Table 6; Figure 6).
The maximum horizontal stress was recalculated taking into account the stress polygon formulation for a reverse faulting regime (Equation (7)). For the crust to be at critical state of stress and considering the range of friction coefficients assumed and the estimated vertical stress and in situ fluid pressure, the maximum horizontal stress should be about 50% greater than estimated based on the stress ratio coefficient (Table 6; Figure 7). The sensitivity analysis suggests that the most influential parameter is the friction coefficient with a Spearman correlation of 90%. The vertical stress has a Spearman correlation of 40% and the in situ fluid pressure of − 10%. These results suggest that the greater the friction coefficient and vertical stress and the smaller the in situ fluid pressure is, the greater the maximum horizontal stress will be (Figure 7). The impact of the probability distribution function chosen to describe each variable has little impact on the results (Table 7; Figure 7).
(a) Stress polygon, (b) maximum horizontal stress distribution and (c) spider plot highlighting the influence of each parameter on the maximum horizontal stress. In (b) red—maximum horizontal stress estimated based on the stress-depth relationships, blue—theoretical maximum horizontal stress.
Possibilities for the maximum principal stress
| Depth (km) | Statistical parameters | 𝜎H (MPa) | Relative difference (%) | |
|---|---|---|---|---|
| Estimated based on stress ratio coefficient (average) | Estimated based on stress polygon (average) | |||
| 1 | Mean | 30.0 | 82.9 | 49 |
| St dev | 3.0 | 14.0 | ||
| Median | 30.3 | 81.8 | ||
| [Min–Max] | [21.9–38.1] | [54.3–126.0] | ||
| 2 | Mean | 56.7 | 164.5 | 54 |
| St dev | 5.0 | 27.6 | ||
| Median | 57.0 | 162.5 | ||
| [Min–Max] | [43.5–70.7] | [107.5–244.7] | ||
| 3 | Mean | 83.4 | 247.9 | 57 |
| St dev | 6.8 | 41.4 | ||
| Median | 83.6 | 244.7 | ||
| [Min–Max] | [64.5–102.5] | [160.1–367.9] | ||
| 4 | Mean | 109.8 | 330.7 | 58 |
| St dev | 8.7 | 55.2 | ||
| Median | 109.9 | 326.5 | ||
| [Min–Max] | [85.1–136.3] | [209.0–492.7] | ||
| 5 | Mean | 135.9 | 412.3 | 59 |
| St dev | 10.7 | 68.7 | ||
| Median | 137.1 | 406.7 | ||
| [Min–Max] | [106.7–168.3] | [270.1–614.5] | ||
𝜎H—maximum horizontal stress.
4.2. Effective stimulation pressure
In situ fluid pressure should be increased by about 50% to induce slip on an optimally oriented fracture (Table 8; Figure 8). The Monte Carlo-based sensitivity analysis suggests that the maximum principal stress (here the maximum horizontal stress) is the most influential parameter and has a negative impact on the results (Figure 8). The Spearman correlation is − 80%. The vertical stress, on the other hand, has a positive correlation with a Spearman coefficient of 60%. Finally, the in situ fluid pressure has a correlation of − 13%. These results suggest that an effective hydraulic stimulation treatment at lower overpressure may be achieved if the maximum principal stress is near its maximum value and if the vertical stress is near its minimum value (Figure 8). The results also suggest the minimal impact of the probability distribution functions (Table 8; Figure 8).
(a) Effective stimulation pressure to induce slip on optimally oriented fractures at a depth of 4 km as an example. In (b) comparison between hydrostatic and overpressure conditions, and in (c) spider plot illustrating the impact of the input parameters on the effective stimulation pressure. In (b) orange—estimated in situ fluid pressure at hydrostatic conditions, brown—estimated in situ fluid pressure at overpressure conditions.
Stimulation fluid pressure
| Depth (km) | Statistical parameters | Pp (MPa) | Relative difference (%) | ||
|---|---|---|---|---|---|
| Hydrostatic conditions (average) | Stimulation fluid pressure (average) | Overpressure conditions | |||
| 1 | Mean | 10.9 | 9.1 | 20.0 | 45 |
| St dev | 0.5 | 3.1 | 3.3 | ||
| Median | 10.9 | 9.1 | 20.0 | ||
| [Min–Max] | [9.7–12.3] | [0–19.2] | [9.5–30.4] | ||
| 2 | Mean | 21.7 | 23.0 | 44.7 | 51 |
| St dev | 0.9 | 5.6 | 6.0 | ||
| Median | 21.7 | 23.0 | 44.7 | ||
| [Min–Max] | [19.4–24.3] | [4.3–42.4] | [24.4–65.7] | ||
| 3 | Mean | 32.6 | 37.4 | 70.0 | 53 |
| St dev | 1.3 | 7.7 | 8.3 | ||
| Median | 32.6 | 37.4 | 70.0 | ||
| [Min–Max] | [29.2–36.2] | [12.0–62.0] | [43.0–96.8] | ||
| 4 | Mean | 43.5 | 50.7 | 94.2 | 54 |
| St dev | 1.8 | 10.6 | 11.5 | ||
| Median | 43.5 | 50.5 | 94.0 | ||
| [Min–Max] | [39.0–48.4] | [17.4–86.3] | [57.3–133.9] | ||
| 5 | Mean | 54.3 | 65.2 | 119.5 | 55 |
| St dev | 2.2 | 13.1 | 14.2 | ||
| Median | 54.3 | 65.3 | 119.6 | ||
| [Min–Max] | [48.6–60.5] | [25.2–107.0] | [75.5–163.9] | ||
Pp—in situ fluid pressure.
4.3. Mohr–Coulomb friction and slip tendency
Mohr–Coulomb friction and slip tendency suggest that if hypothesis H1 (Table 4) prevails in Kuujjuaq, then the favorable sets to slip at the estimated critical fluid pressure (Equation (9)) are the E–W and NW–SE (Figure 9a). However, if in Kuujjuaq the principal stresses are oriented according to hypothesis H2 (Table 4), then the optimally oriented sets are E–W and N–S (Figure 9b). Finally, if the maximum principal stress in Kuujjuaq is oriented N–S (hypothesis H3; Table 4), then the fracture sets NW–SE and NE–SW have the greatest tendency to slip (Figure 9c).
Optimally oriented fractures to slip at the effective stimulation pressure estimated for a depth of 4 km and for the different hypotheses assumed for the orientation of the principal stresses. Warm colors—high tendency to slip, cold colors—low tendency to slip, dots—poles to planes of each fracture.
5. Discussion
Canada’s northern and remote communities have a critical dependence on energy services for their safety, sustainability, and economic growth [Gilmour et al. 2018]. Unconventional geothermal development, such as EGS, has potential to provide sustainable energy by exploiting deep energy sources of petrothermal systems within the Canada Shield [e.g., Grasby et al. 2012]. Given the opportunity that off-grid renewable solutions may provide to Canadian northern communities, this work carried out in the community of Kuujjuaq aims at providing a first-order approximation to the subsurface state of stress, which is a key geomechanical parameter for the successful development of hydraulically stimulated geothermal systems [e.g., Brown et al. 2012; Evans et al. 1999; Richards et al. 1994]. The implications for EGS design and development and the strengths and weaknesses of this study, together with future envisioned work, are discussed in this section.
5.1. Implications for EGS design and development
A hydraulically stimulated geothermal reservoir tends to grow in the direction of the maximum principal stress [e.g., Brown et al. 2012; Evans et al. 1999; Hubbert and Willis 1957]. Thus, based on the hypotheses assumed for the orientation of the principal stresses, a reservoir developed in Kuujjuaq could have its major axis parallel to NE–SW or N–S. Moreover, based on the field tests carried out in, for instance at Rosemanowes [e.g., Pine and Batchelor 1984; Richards et al. 1994], the performance of the system may be improved if the open hole section of the wells is drilled with its azimuth parallel to the minimum horizontal principal stress. Considering the hypothses studied, it could be NW–SE or E–W. If the azimuths of the open hole section of the wells are drilled at right angles to the minimum horizontal principal stress, this may lead to the interception of few conductive fracture sets. In Kuujjuaq, based on the field data gathered by Miranda [2021] and the stress field estimated in this work, the fracture sets E–W and NE–SW have the most optimal orientation to slip if the maximum horizontal stress is around N248°/10°. However, if the maximum horizontal stress has an orientation N227°/02°, the optimally oriented sets are E–W and N–S. Finally, if the maximum horizontal stress is oriented N–S, then the optimally oriented sets are NW–SE and NE–SW.
However, despite the potential advantages of drilling deviated wells to tap into the elongated stimulated volume and well-conductive fractures, this may bring wellbore stability problems. The collapse of the borehole is something that should be avoided since it will have important financial impacts on the project. Thus, deviated wells have pros and cons and more research is needed to understand the risk of borehole failure if a deviated well is drilled in Kuujjuaq to intercept the conductive structures.
The in-situ stress conditions are also found to influence the development of stimulated rock regions or migration of induced seismicity [e.g., Xie et al. 2015]. Results from sites where hydraulic stimulation has been applied suggest that the seismic cloud appears vertical to sub-vertical in strike-slip faulting regime and horizontal to sub-horizontal in reverse faulting stress conditions [Xie et al. 2015]. Thus, based on the results of this first-order analysis, the seismic cloud could appear horizontal to sub-horizontal. A multiple-well system may be a good configuration for the reverse faulting stress regime [Xie et al. 2015].
Not only does the stress regime influences the stimulated region but also the stress difference impacts the stimulation performance [e.g., Xie et al. 2015]. Cases where large differential stress exists tend to require less additional fluid pressure to activate shear slip of natural fractures [Xie et al. 2015]. For example, the differential stress at Hijiori is smaller compared to Fenton Hill, Soultz or Rosemanowes which led to greater effort to trigger shear slip [Xie et al. 2015]. The analysis carried out in this work suggest that the maximum principal stress should be 50% greater than estimated based on stress-depth relationships for the subsurface to be at critical state of stress. Nevertheless, the stimulation fluid pressure estimated in this work, and for the same depth interval, is comparable with the test sites reviewed by Xie et al. [2015]. However, the results obtained suggest that developing EGS in Kuujjuaq may be easier if the hydraulically stimulated reservoir is developed at shallower depths where the optimally oriented fracture sets may be near critical conditions, considering the distribution of values given by the Monte Carlo simulations. Previous studies carried out by Miranda et al. [2020b] to evaluate the subsurface temperature distribution suggested that at 4 km depth the geothermal energy source in place has a great probability (98%) of meeting Kuujjuaq’s annual average space heating demand (37 GWh). However, the study of Miranda et al. [2020b] highlighted that a large uncertainty exists about the estimated temperature. In fact, space heating may be produced at shallower depths (of 2–3 km) if the subsurface temperature is near the maximum estimated value. This would not only be advantageous from the technical point of view (as explained above) but also from an economic perspective.
Ultimately, the necessary flow rate for the hydraulic stimulation treatment can be approximated using the maximum injection rate versus the associated wellhead pressure plot of Xie et al. [2015] and taking into account the common assumed and typically obtained hydraulic impedance. For a hydraulically stimulated reservoir located at 4 km depth, a fluid stimulation pressure of about 51 MPa may be needed which suggests a maximum flow rate ranging between 45 to more than 200 L⋅s−1 to meet the hydraulic impedance of 1 and 10 L⋅s−1⋅MPa−1. At 2 km depth, the required stimulation pressure diminishes and, consequently, the maximum flow rate.
All these observations, however, are valid if the mean value given by the Monte Carlo simulations prevails. As the simulations revealed, there is a large possible distribution for the values and more effort is needed to reduce the uncertainty highlight throughout this work.
5.2. Strengths and weaknesses of the theoretical approach
Estimating in situ stresses from theoretical and empirical relationships has been a challenging task for many researchers [Taherynia et al. 2016]. Although this approach has certain advantages due to the difficulty and costs of carrying out in situ stress measuring tests, a correct calibration of the estimates with previous stress measurements done in the region is necessary [e.g., Amadei and Stephansson 1997]. Unfortunately, no stress measurements have been made to date in Kuujjuaq which is a challenge for the application of the theoretical approach. In fact, the stress regime estimated in this study is the target of high uncertainty since the orientation of the principal stresses was taken from literature and their magnitude was calculated using stress-depth relationships that are based on observations made in different geologic and tectonic contexts. Thus, planning the design and development of EGS without actual stress measurements is highly speculative and may lead to severe consequences for the geothermal project. The research work presented in this study should thus be seen as a first-order approximation to highlight the need of further geothermal research. Nevertheless, this study suggests the existence of optimally oriented fracture sets that may be at critical state of stress if the magnitude of the maximum principal stress is about 50% greater than estimated based on the stress-depth relationships, or if the in-situ fluid pressure is greater than hydrostatic regime. Although the latter may be unlikely since most in situ fluid pressure assessments in several deep boreholes indicate hydrostatic regime up to 12 km depth [e.g., Townend and Zoback 2000; Zoback and Townend 2001]. Moreover, estimates of optimally oriented sets were undertaken, and these are strongly dependent on the orientation of the maximum horizontal stress. The optimally oriented sets could be E–W and NE–SW, but also E–W and N–S and NW–SE and NE–SW as well. This consequently has a strong influence on the azimuth of the wells. If the optimally oriented sets are E–W and NE–SW or E–W and N–S, then the wells should be drilled oriented NW–SE to intercept the conductive sets. If, on the contrary, the optimally oriented sets are NW–SE and NE–SW, then this orientation for the wells would be parallel to one of the main conductive sets. This highlights how important the knowledge of the stress field is to properly design an EGS. Similarly, the required fluid overpressure to initiate shear slip is also highly dependent on the differential stress. The values obtained in this study may be reduced if the differential stress is greater than here estimated.
In summary, although highly speculative, the theoretical approach put forward in this work is a useful tool for an a priori estimation of the stress regime but cannot replace actual stress measurements that will have to be carried out for the next exploration steps. Nevertheless, the findings of this research work seem to justify the need for further geothermal exploration in order to resolve the uncertainty found and accurately design and develop EGS in remote northern regions.
In situ stress measurements are the next step to calibrate the results obtained in this work. In fact, the findings of this research work should be integrated with stress measurements done in geothermal exploratory boreholes [e.g., Kruszewski et al. 2021]. This will help to decrease the uncertainty here found and design more accurately, not only the reservoir dilatation strategy to employ, but also the drilling of injectors and producers.
An approach strategy for rock stress estimation has been developed by the International Society of Rock Mechanics (ISRM) and it should be followed in future in situ stress measurements [Hudson et al. 2003]. For example, inverse focal plane solutions or seismic shear wave anisotropy can be used to estimate the principal stress directions and the ratio of stress differences. Hydraulic or drilling induced fractures and borehole breakout orientations can help to establish the minimum principal stress orientation. Indirect methods on borehole core, such as the Kaiser effect and differential strain analysis, can be used to find the components of the stress tensor. Reservoir dilation tests in boreholes can establish the minimum and maximum principal stresses. Overcoring tests and hydraulic testing of pre-existing fractures can be used to establish the complete stress state. Finally, numerical analyses can help to establish the variation of the stress state across the site due to different geological strata and fractures.
6. Conclusions
Remote, northern and rural communities in Canada and worldwide can highly benefit from off-grid renewable energy solutions to improve their sustainability and living standards and boost their economic growth. Geothermal energy sources harvested via engineered geothermal energy systems can be one of such viable options if current exploration challenges are overcome. A key parameter for the correct design of such geothermal systems is the stress regime. Thus, in this work, a priori estimation of the stress regime was undertaken in an Arctic off-grid community (Kuujjuaq, Canada) to provide guidelines and an example for the remaining settlements located in a similar geological context.
The stress regime was estimated following a theoretical approach that uses the World Stress Map and other literature data to evaluate the orientation of the principal stresses and empirical correlations to calculate the magnitude of the principal stresses, in-situ fluid pressure and limits for the principal stresses. Monte Carlo simulations were carried out to deal with the uncertainty of each input parameter and the results obtained were compared with literature data for the Canadian Shield. Additionally, Mohr–Coulomb friction and slip tendency analyses were carried out to estimate which fracture sets previously sampled are optimally oriented to slip considering the uncertainty associated with the orientation of the maximum horizontal stress. This is also helpful to determine the wells drilling orientation and stimulation pressure to initiate shear slip.
A transpression deformation style may be possible based on the findings of this work. Considering the orientation of the principal stresses used in this work, the fracture sets optimally oriented to slip can be E–W and NE–SW, or E–W and N–S, or NW–SE and NE–SW. However, these are not at critical state of stress. Given such observation, the optimal azimuth for the wells could be NW–SE or E–W in order to intercept the conductive sets. However, further efforts are needed to understand potential borehole failure. Moreover, developing an EGS in Kuujjuaq may be more effective if if the reservoir is located at depths of 2–3 km, where the optimal oriented fractures may be closer to a critical state of stress.
In conclusion, the results obtained in this work provide a first-order approximation to the state of stress of the subsurface beneath Kuujjuaq, that although highly speculative due to the lack of stress measurements and calibration provide an a priori planning for the design and development of EGS in this remote northern community.
Notation
| Symbol | Definition | Unit |
|---|---|---|
| F | Factor | — |
| FBP | Formation breakdown pressure | Pa |
| g | Gravitational acceleration | m⋅s−2 |
| k | Stress ratio coefficient | — |
| P | Pressure | Pa |
| S | Slip tendency | — |
| T | Tensile strength | Pa |
| z | Depth | m |
| Greek letters | ||
| 𝜇 | Friction coefficient | — |
| 𝜌 | Density | kg⋅m−3 |
| 𝜎 | Stress | Pa |
| 𝜎′ | Effective stress | Pa |
| 𝜏 | Shear stress | Pa |
| Subscript | ||
| 1 | Maximum principal | |
| 3 | Minimum principal | |
| H | Maximum horizontal | |
| h | Minimum horizontal | |
| max | Maximum | |
| min | Minimum | |
| n | Normal | |
| p | In-situ fluid | |
| p − f | Pore to fluid | |
| V | Vertical | |
Conflicts of interest
Authors have no conflict of interest to declare.
Acknowledgments
This study was funded by the Institut nordique du Québec (INQ) through the Chaire de recherche sur le potentiel géothermique du Nord awarded to Jasmin Raymond. The Centre d’études nordiques (CEN), supported by the Fonds de recherche du Québec—nature et technologies (FRQNT), and the Observatoire Homme Milieu Nunavik (OHMI) are further acknowledged for helping with field campaigns cost and logistics. The authors are also grateful to the two anonymous reviewers whose comments helped to greatly improved the manuscript.