# Estimating theoretical stress regime for engineered geothermal energy systems in an arctic community
(Kuujjuaq, Canada)

Mafalda M Miranda, Jasmin Raymond, Chrystel Dezayes

► To cite this version:

Mafalda M Miranda, Jasmin Raymond, Chrystel Dezayes. Estimating theoretical stress regime for engineered geothermal energy systems in an arctic community (Kuujjuaq, Canada). Comptes Rendus. Géoscience, 2023, 355, pp.85 - 108. 10.5802/crgeos.193. hal-04001170

HAL Id: hal-04001170

https://brgm.hal.science/hal-04001170

Submitted on 22 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Géoscience
### Sciences de la Planète

Mafalda M. Miranda, Jasmin Raymond and Chrystel Dezayes

Estimating theoretical stress regime for engineered geothermal energy
systems in an arctic community (Kuujjuaq, Canada)

Volume 355 (2023), p. 85-108

Published online: 1 February 2023

https://doi.org/10.5802/crgeos.193

<figure><img src="image_2.png" /></figure>

©

<figure><img src="image_3.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_4.png" /></figure>

Les Comptes Rendus. Géoscience — Sciences de la Planète sont membres du
Centre Mersenne pour l'édition scientifique ouverte

www.centre-mersenne.org

e-ISSN : 1778-7025
Original Article — Tectonics, tectonophysics, geodynamics

# Estimating theoretical stress regime for engineered geothermal energy systems in an arctic community (Kuujjuaq, Canada)

Mafalda M. Miranda<sup>●, *, a</sup>, Jasmin Raymond<sup>● a</sup> and Chrystel Dezayes<sup>● b</sup>

<sup>a</sup> INRS—Institute national de la recherche scientifique, 490, rue de la Couronne
Québec (QC) G1K 9A9, Canada

<sup>b</sup> BRGM, 3 avenue Claude Guillemin, BP 36009, 45060, Orléans Cedex 2, France

*E-mails*: mafalda_alexandra.miranda@ete.inrs.ca, mfldmiranda@gmail.com
(M. M. Miranda), jasmin.raymond@inrs.ca (J. Raymond), c.dezayes@brgm.fr
(C. Dezayes)

**Abstract.** In remote northern regions, lacking deep geothermal exploratory boreholes, a theoretical approach to provide a first-order estimate of the stress regime seems a useful tool. Literature data was used in this context to evaluate the orientation of the stress components and empirical relationships were applied to calculate their magnitude in a community of Nunavik, northern Quebec. A Monte Carlo-based sensitivity analysis was carried out due to the uncertainty of the input parameters. Mohr-Coulomb friction and slip tendency analyses were additionally undertaken to assess the stress state and potential reactivation of existing fractures. The results highlight how the poor knowledge of the stress field has an important impact on the design and development of engineered geothermal energy systems in the Canadian off-grid community of Kuujjuaq.

**Keywords.** Stress regime, Mohr-Coulomb criterion, Slip tendency, Fluid overpressure, Geothermal energy, Nunavik.

Manuscript received 8 December 2021, revised 25 November 2022, accepted 28 November 2022.

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

<footnote>* Corresponding author.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** 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].</figcaption></figure>

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. [2021b] 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⁻² 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. [2021a], 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].

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⁻¹ for the upper crustal levels (from 0 to 10 km depth), between 6.2 and 6.5 km·s⁻¹ for the mid-crust (from 10 to about 17 km depth) and between 6.6 and 7.0 km·s⁻¹ 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 [subphilititic 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).

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.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** (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.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**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].</figcaption></figure>
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⁻¹ [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⁻¹ while the gradient for the maximum horizontal stress between 34.4 and 51.3 MPa·km⁻¹ [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].

## 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]:

$$ \sigma_V = \rho g z \quad (1) $$

where $\sigma_V$ (Pa) is the vertical stress, $\rho$ (kg·m⁻³) is the density of the geological materials, $g$ (m·s⁻²) is the gravitational acceleration and $z$ (m) is depth.

The horizontal stress components were estimated based on the horizontal to vertical stress ratio as [e.g., Zhang, 2017]:

$$ \begin{cases} \sigma_H = k_{\max} \sigma_V \\ \sigma_h = k_{\min} \sigma_V \end{cases} \quad (2) $$

where $k$ (-) is the stress ratio coefficient and $\sigma_H$ and $\sigma_h$ (Pa) are the maximum and minimum horizontal stresses, respectively. Stress ratio expressions obtained for the Canadian Shield (Table 2) were used in this work.

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:

$$ R = \frac{\sigma_2 - \sigma_3}{\sigma_1 - \sigma_3} \Leftrightarrow R = \frac{\sigma_h - \sigma_V}{\sigma_H - \sigma_V} \quad (3) $$

where $\sigma_1$ (Pa) is the maximum principal stress, $\sigma_2$ (Pa) is the intermediate principal stress and $\sigma_3$ (Pa) is the minimum principal stress.
<table><caption>**Table 1.** Principal stresses in the Canadian Shield</caption><thead><tr><th>Principal stress</th><th>Orientation</th><th>Magnitude</th><th>Observations</th><th>Reference</th></tr></thead><tbody><tr><td>σ<sub>V</sub></td><td>—</td><td>(0.0260–0.0324)z<br>9.86 + 0.0371z</td><td>0 &lt; z &lt; 2200 m<br>0 &lt; z &lt; 900 m<br>σ<sub>V</sub> &lt; σ<sub>H,average</sub></td><td>Herget [1982, 1987]</td></tr><tr><td>σ<sub>H,average</sub></td><td>—</td><td>33.41 + 0.0111z</td><td>900 &lt; z &lt; 2200 m</td><td></td></tr><tr><td>σ<sub>V</sub></td><td>—</td><td>(0.0266 ± 0.008)z</td><td></td><td></td></tr><tr><td>σ<sub>H,average</sub></td><td>—</td><td>5.91 + 0.0349z</td><td>60 &lt; z &lt; 1890 m</td><td>Arjang [1991]</td></tr><tr><td>σ<sub>H</sub></td><td>—</td><td>8.18 + 0.0422z</td><td></td><td></td></tr><tr><td>σ<sub>h</sub></td><td>—</td><td>3.64 + 0.0276z</td><td></td><td></td></tr><tr><td>σ<sub>V</sub></td><td>—</td><td>0.0285z</td><td></td><td></td></tr><tr><td>σ<sub>1</sub></td><td>N248°/10°</td><td>12.1 + (0.0403 ± 0.0020)z</td><td>0 &lt; z &lt; 2200 m</td><td>Herget [1993]</td></tr><tr><td>σ<sub>2</sub></td><td>N300–340°/0°</td><td>6.4 + (0.0293 ± 0.0019)z</td><td></td><td></td></tr><tr><td>σ<sub>3</sub></td><td>Vertical</td><td>1.4 + (0.0225 ± 0.0015)z</td><td></td><td></td></tr><tr><td>σ<sub>V</sub></td><td>—</td><td>0.0260z</td><td></td><td></td></tr><tr><td>σ<sub>1</sub></td><td>NE/horizontal</td><td>13.50 + 0.0344z</td><td>0 &lt; z &lt; 6000 m</td><td>Arjang [1998]</td></tr><tr><td>σ<sub>2</sub></td><td>NW/sub-horizontal</td><td>8.03 + 0.0233z</td><td></td><td></td></tr><tr><td>σ<sub>3</sub></td><td>Vertical</td><td>3.01 + 0.0180z</td><td></td><td></td></tr><tr><td>σ<sub>1</sub></td><td>N-S/horizontal</td><td>0.0513z</td><td></td><td></td></tr><tr><td>σ<sub>2</sub></td><td>E-W/horizontal</td><td>0.0378z</td><td>—</td><td>Andrieux et al. [2009]</td></tr><tr><td>σ<sub>3</sub></td><td>Vertical</td><td>0.0270z</td><td></td><td></td></tr><tr><td>σ<sub>V</sub></td><td>—</td><td>(0.0258–0.0263)z</td><td></td><td></td></tr><tr><td>σ<sub>1</sub></td><td>N227°/02°</td><td>(0.040 ± 0.001)z – (9.185 ± 1.5)</td><td>12 &lt; z &lt; 2552 m</td><td>Young and Maloney [2015]</td></tr><tr><td>σ<sub>2</sub></td><td>N310°/08°</td><td>(0.029 ± 0.001)z + (4.617 ± 1.159)</td><td></td><td></td></tr><tr><td>σ<sub>3</sub></td><td>N270°/88°</td><td>(0.021 ± 0.001)z – (0.777 ± 0.872)</td><td></td><td></td></tr><tr><td>σ<sub>V</sub></td><td>—</td><td>0.021z</td><td>0 &lt; z &lt; 1300 m</td><td></td></tr><tr><td>σ<sub>1</sub></td><td>—</td><td>0.012z + 42.4</td><td></td><td></td></tr><tr><td>σ<sub>2</sub></td><td>—</td><td>0.013z + 24.1</td><td>660 &lt; z &lt; 1300 m</td><td>Kaiser et al. [2016]</td></tr><tr><td>σ<sub>3</sub></td><td>—</td><td>0.007z + 9.7</td><td></td><td></td></tr></tbody></table>

σ<sub>V</sub>—vertical stress, σ<sub>H</sub>—maximum horizontal stress, σ<sub>h</sub>—minimum horizontal stress, σ<sub>H,average</sub>—average horizontal stress (σ<sub>H</sub> + σ<sub>h</sub>)/2, σ<sub>1</sub>—maximum principal stress, σ<sub>2</sub>—intermediate principal stress, σ<sub>3</sub>—minimum principal stress.

The in situ fluid pressure can be estimated based on the pore-fluid factor [e.g., Sibson, 2004, Zoback and Townend, 2001]:

$$P_p = F_{p-f} \sigma_v \quad (4)$$

where P<sub>p</sub> (Pa) is the in situ fluid pressure and F<sub>p-f</sub> is the pore to fluid factor. The latter is approximately 0.4 assuming a hydrostatic regime.

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]:

$$\frac{\sigma'_{1}}{\sigma'_{3}} \le \left( \sqrt{\mu^2 + 1} + \mu \right)^2 \quad (5)$$

where σ'<sub>1</sub> and σ'<sub>3</sub> (Pa) are the maximum and minimum effective principal stresses, respectively, and μ is the friction coefficient. This theoretical criterion was used to assess if the estimated stresses do not exceed the theoretical frictional strength of pre-existing
<table><caption>**Table 2.** Horizontal to vertical stress ratios inferred for the Canadian Shield</caption><thead><tr><th>Stress ratio coefficient</th><th>Expression</th><th>Reference</th></tr></thead><tbody><tr><td>k<sub>max</sub></td><td>(357/z) + 1.46</td><td>Herget [1987]</td></tr><tr><td>k<sub>min</sub></td><td>(167/z) + 1.10</td><td>Herget [1987]</td></tr><tr><td>k<sub>max</sub></td><td>((272 ± 8)/z) + 1.72</td><td>Herget [1993]</td></tr><tr><td>k<sub>min</sub></td><td>((30 ± 4)/z) + 0.86</td><td>Herget [1993]</td></tr><tr><td>k<sub>max</sub></td><td>7.44z<sup>-0.198</sup></td><td>Arjang [1998]</td></tr><tr><td>k<sub>min</sub></td><td>2.81z<sup>-0.120</sup></td><td>Arjang [1998]</td></tr></tbody></table>

k<sub>max</sub>—maximum horizontal to vertical stress ratio, k<sub>min</sub>—minimum horizontal to vertical stress ratio, z—depth in meters.

discontinuities assuming a reverse faulting regime typical for the Canadian Shield (Table 1).

The maximum and minimum effective principal stresses, assuming reverse faulting regime, are calculated as:

$$ \begin{cases} \sigma'_{1} = \sigma_{H} - P_{p} \\ \sigma'_{3} = \sigma_{H} - P_{p} \end{cases} \qquad (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]:

$$ \left\{ \begin{aligned} f(\mu) &= \left(\sqrt{\mu^2 + 1} + \mu\right)^2 \\ \sigma_H &= \sigma_h \\ \frac{(\sigma_V - P_p)}{f(\mu)} + P_p &= \sigma_h - \text{normal faulting regime} \\ \frac{(\sigma_H - P_p)}{(\sigma_h - P_p)} &= f(\mu) - \text{strike-slip faulting regime} \\ f(\mu)(\sigma_V - P_p) + P_p &= \sigma_H - \text{reverse faulting regime}. \end{aligned} \right. \qquad (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⁻³ and the maximum value of 3108 kg·m⁻³, with a median value of 2677 kg·m⁻³. 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.

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
<table><caption>**Table 3.** Monte Carlo method input parameters and their uncertainty</caption><thead><tr><th>Parameter code</th><th>Parameter description</th><th>Unit</th><th>Variable type</th><th>Distribution</th></tr></thead><tbody><tr><td>ρ</td><td>Density</td><td>kg·m<sup>-3</sup></td><td>Continuous</td><td>Triangular (min, med, max)</td></tr><tr><td>g</td><td>Gravitational acceleration</td><td>m·s<sup>-2</sup></td><td>—</td><td>Single value</td></tr><tr><td>z</td><td>Depth</td><td>m</td><td>—</td><td>Single value</td></tr><tr><td>k</td><td>Stress ratio coefficient</td><td>—</td><td>Discrete continuous</td><td>Uniform (min, max)</td></tr><tr><td>σ<sub>V</sub></td><td>Vertical stress</td><td>Pa</td><td>Continuous</td><td>Triangular (min, mean, max)<br>Uniform (min, max)</td></tr><tr><td>F<sub>p-f</sub></td><td>Pore to fluid factor</td><td>—</td><td>—</td><td>Single value</td></tr><tr><td>σ<sub>H</sub></td><td>Maximum horizontal stress</td><td>Pa</td><td>Continuous</td><td>Triangular (min, mean, max)<br>Uniform (min, max)</td></tr><tr><td>μ</td><td>Friction coefficient</td><td>—</td><td>Continuous</td><td>Uniform (min, max)</td></tr><tr><td>P<sub>p</sub></td><td>In situ fluid pressure</td><td>Pa</td><td>Continuous</td><td>Triangular (min, mean, max)<br>Uniform (min, max)</td></tr></tbody></table>

[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]:

$$FBP = 3\sigma_h - \sigma_H - P_p + T \quad (8)$$

where FBP (Pa) is the formation breakdown pressure and T (Pa) is the tensile strength. The latter is assumed zero in this work following the values described by Andrieux et al. [2009]. However, this criterion gives an upper bound on the fluid overpressure

required to hydraulically fracture rock formation. The actual stimulation pressure may be less than the value estimated based on the aforementioned equation due to, among other features, the presence of pre-existing natural fractures. Thus, the critical fluid pressure required to induce slip on an optimally oriented fracture considering a friction coefficient of 0.6 was estimated as [e.g., Rutqvist et al., 2007, Taghipour et al., 2019]:

$$P_{p,critical} = \frac{3\sigma'_{3} - \sigma'_{1}}{2}. \quad (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

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 reactivation [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
<table><caption>**Table 4.** Hypotheses tested in the Mohr–Coulomb friction and slip tendency analyses</caption><thead><tr><th rowspan="2">Principal stresses</th><th colspan="4">Tested hypotheses</th></tr><tr><th>H1</th><th>H2</th><th>H3</th><th>Effective magnitude<br>(mean value; MPa)</th></tr></thead><tbody><tr><td>σ<sub>1</sub></td><td>N248°/10°</td><td>N227°/02°</td><td>N-S</td><td>94</td></tr><tr><td>σ<sub>2</sub></td><td>N340°/0°</td><td>N310°/08</td><td>E-W</td><td>66</td></tr><tr><td>σ<sub>3</sub></td><td>Vertical</td><td>N270°/88°</td><td>Vertical</td><td>64</td></tr></tbody></table>

σ₁—maximum principal stress, σ₂—intermediate principal stress, σ₃—minimum principal stress.

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

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]:

$$S = \frac{\tau}{\sigma'_{n}} \qquad (10)$$

where S (-) is the tendency to slip, τ (Pa) is the shear stress and σ'_{n} (Pa) is the effective normal stress.

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

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:

$$ \left\{ \begin{aligned} \sigma_V (\text{MPa}) &= (0.027 \pm 0.03)z \\ \sigma_H (\text{MPa}) &= (9.5 \pm 3.5) + (0.033 \pm 0.012)z \\ \sigma_h (\text{MPa}) &= (2.5 \pm 2.5) + (0.027 \pm 0.06)z \\ P_p (\text{MPa}) &= (0.011 \pm 0.01)z. \end{aligned} \right. \qquad (11) $$

The estimated vertical stress suggests a gradient of 27 ± 3 MPa·km⁻¹ and the gradient for the maximum and minimum horizontal stress was evaluated as 33 ± 12 MPa·km⁻¹ and 27 ± 6 MPa·km⁻¹, respectively. The gradient for the in situ fluid pressure was estimated as 11 ± 0.1 MPa·km⁻¹. It is important to highlight that these gradients were obtained when plotting the stresses to a depth of 10 km. The gradients change depending on the depth range used. The value on the left side of the ± sign corresponds to the mean while the value on the right side of the ± sign corresponds to the range between minimum and maximum.

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.
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

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.

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).
<table><caption>Table 5. Estimated stress components</caption><thead><tr><th rowspan="3">Depth (km)</th><th rowspan="3">Statistical parameters</th><th colspan="10">Principal stresses</th></tr><tr><th colspan="3">σ<sub>V</sub> (MPa)</th><th colspan="3">σ<sub>H</sub> (MPa)</th><th colspan="3">σ<sub>H</sub> (MPa)</th><th colspan="3">σ<sub>H</sub> (MPa)</th></tr><tr><th>A</th><th>B</th><th>C</th><th>D</th><th>Average</th><th>A</th><th>B</th><th>C</th><th>D</th><th>Average</th><th>A</th><th>B</th><th>Average</th></tr></thead><tbody><tr><td rowspan="4">1</td><td>Mean</td><td>27.0</td><td>30.4</td><td>31.0</td><td>29.1</td><td>29.7</td><td>30.0</td><td>43.6</td><td>44.4</td><td>39.8</td><td>40.5</td><td>42.1</td><td>10.9</td><td>11.0</td><td>10.9</td></tr><tr><td>St dev</td><td>1.3</td><td>4.8</td><td>5.2</td><td>3.3</td><td>3.7</td><td>3.0</td><td>9.5</td><td>10.0</td><td>5.7</td><td>6.2</td><td>5.7</td><td>0.6</td><td>0.8</td><td>0.5</td></tr><tr><td>Median</td><td>26.8</td><td>32.3</td><td>32.1</td><td>29.0</td><td>29.5</td><td>30.3</td><td>48.1</td><td>47.7</td><td>39.7</td><td>40.3</td><td>42.9</td><td>10.9</td><td>11.0</td><td>10.9</td></tr><tr><td>(Min-Max)</td><td>[24.2-30.4]</td><td>[21.6-38.6]</td><td>[21.4-39.3]</td><td>[22.0-38.3]</td><td>[21.4-39.2]</td><td>[21.9-38.1]</td><td>[27.4-57.6]</td><td>[27.2-58.7]</td><td>[27.8-54.5]</td><td>[27.4-56.2]</td><td>[28.1-56.2]</td><td>[9.6-12.4]</td><td>[9.6-12.4]</td><td>[9.7-12.3]</td></tr><tr><td rowspan="4">2</td><td>Mean</td><td>54.0</td><td>57.3</td><td>57.9</td><td>55.5</td><td>56.1</td><td>56.7</td><td>77.1</td><td>77.9</td><td>71.1</td><td>71.8</td><td>74.5</td><td>21.7</td><td>21.8</td><td>21.7</td></tr><tr><td>St dev</td><td>2.6</td><td>7.8</td><td>8.4</td><td>5.6</td><td>6.2</td><td>5.0</td><td>16.9</td><td>17.6</td><td>10.6</td><td>11.3</td><td>10.2</td><td>1.1</td><td>1.5</td><td>0.9</td></tr><tr><td>Median</td><td>53.7</td><td>59.9</td><td>59.2</td><td>55.4</td><td>55.8</td><td>57.0</td><td>85.5</td><td>84.4</td><td>71.1</td><td>71.3</td><td>75.7</td><td>21.7</td><td>21.8</td><td>21.7</td></tr><tr><td>(Min-Max)</td><td>[48.4-60.9]</td><td>[42.4-72.1]</td><td>[42.0-72.2]</td><td>[42.8-71.5]</td><td>[42.1-72.0]</td><td>[43.5-70.7]</td><td>[48.4-100.6]</td><td>[47.8-100.8]</td><td>[49.2-98.5]</td><td>[48.0-99.6]</td><td>[49.0-97.9]</td><td>[19.2-24.4]</td><td>[19.2-24.4]</td><td>[19.4-24.3]</td></tr><tr><td rowspan="4">3</td><td>Mean</td><td>81.0</td><td>83.7</td><td>84.8</td><td>82.0</td><td>83.1</td><td>83.4</td><td>109.4</td><td>110.8</td><td>102.4</td><td>103.7</td><td>106.6</td><td>32.7</td><td>32.8</td><td>32.6</td></tr><tr><td>St dev</td><td>3.9</td><td>10.6</td><td>11.3</td><td>7.8</td><td>8.6</td><td>6.8</td><td>23.6</td><td>24.4</td><td>15.5</td><td>16.4</td><td>14.4</td><td>1.5</td><td>2.1</td><td>1.3</td></tr><tr><td>Median</td><td>80.5</td><td>86.3</td><td>86.4</td><td>81.7</td><td>82.8</td><td>83.6</td><td>120.6</td><td>120.3</td><td>102.1</td><td>103.2</td><td>108.1</td><td>32.6</td><td>32.8</td><td>32.6</td></tr><tr><td>(Min-Max)</td><td>[72.6-91.3]</td><td>[63.3-105.5]</td><td>[63.5-105.2]</td><td>[64.0-104.5]</td><td>[63.9-104.9]</td><td>[64.5-102.5]</td><td>[69.3-144.1]</td><td>[69.4-143.7]</td><td>[70.4-143.2]</td><td>[69.5-142.8]</td><td>[71.3-139.0]</td><td>[29.2-36.3]</td><td>[29.2-36.4]</td><td>[29.2-36.2]</td></tr><tr><td rowspan="4">4</td><td>Mean</td><td>107.9</td><td>109.6</td><td>111.2</td><td>108.4</td><td>110.0</td><td>109.8</td><td>140.9</td><td>143.0</td><td>133.7</td><td>135.6</td><td>138.3</td><td>43.6</td><td>43.8</td><td>43.5</td></tr><tr><td>St dev</td><td>5.2</td><td>13.4</td><td>14.4</td><td>10.0</td><td>11.4</td><td>8.7</td><td>30.1</td><td>31.1</td><td>20.4</td><td>21.7</td><td>18.7</td><td>2.0</td><td>2.9</td><td>1.8</td></tr><tr><td>Median</td><td>107.4</td><td>111.5</td><td>112.2</td><td>108.2</td><td>109.5</td><td>109.9</td><td>154.0</td><td>154.1</td><td>133.3</td><td>135.3</td><td>140.0</td><td>43.5</td><td>43.8</td><td>43.5</td></tr><tr><td>(Min-Max)</td><td>[96.8-121.8]</td><td>[84.0-139.0]</td><td>[84.2-139.3]</td><td>[84.4-137.9]</td><td>[84.4-139.0]</td><td>[85.1-136.3]</td><td>[89.9-188.6]</td><td>[90.0-189.0]</td><td>[91.2-186.3]</td><td>[90.5-188.1]</td><td>[93.1-185.3]</td><td>[38.9-48.8]</td><td>[38.8-48.8]</td><td>[39.0-48.4]</td></tr><tr><td rowspan="4">5</td><td>Mean</td><td>134.9</td><td>135.4</td><td>137.0</td><td>134.9</td><td>136.5</td><td>135.9</td><td>171.9</td><td>174.0</td><td>165.0</td><td>166.9</td><td>169.5</td><td>54.4</td><td>54.6</td><td>54.3</td></tr><tr><td>St dev</td><td>6.6</td><td>16.2</td><td>17.4</td><td>12.3</td><td>13.8</td><td>10.7</td><td>36.2</td><td>37.6</td><td>25.4</td><td>26.7</td><td>22.7</td><td>2.5</td><td>3.6</td><td>2.2</td></tr><tr><td>Median</td><td>134.2</td><td>135.6</td><td>137.5</td><td>134.5</td><td>135.9</td><td>137.1</td><td>184.5</td><td>186.9</td><td>164.5</td><td>166.5</td><td>171.2</td><td>54.3</td><td>54.6</td><td>54.3</td></tr><tr><td>(Min-Max)</td><td>[121.0-152.3]</td><td>[104.8-172.4]</td><td>[104.8-172.3]</td><td>[105.9-170.0]</td><td>[104.9-172.1]</td><td>[106.7-168.3]</td><td>[110.6-233.1]</td><td>[110.6-232.8]</td><td>[112.5-231.2]</td><td>[111.1-232.1]</td><td>[112.8-225.8]</td><td>[48.5-60.7]</td><td>[48.4-60.8]</td><td>[48.6-60.5]</td></tr></tbody></table>

σ<sub>V</sub>—vertical stress, σ<sub>H</sub>—maximum horizontal stress, σ<sub>H</sub>—minimum horizontal stress, p<sub>p</sub>—in situ fluid pressure, A—σ<sub>V</sub> triangular distribution and k discrete distribution, C—σ<sub>V</sub> triangular distribution and k continuous distribution, D—σ<sub>V</sub> uniform distribution and k discrete distribution, C—σ<sub>V</sub> uniform distribution and k discrete distribution.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** 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.</figcaption></figure>

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

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

## 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
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Effective stress ratio for a depth of 4 km as an example.</figcaption></figure>

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

**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
<table><caption>**Table 6.** Estimated effective stress ratio coefficient</caption><thead><tr><th rowspan="2">Depth (km)</th><th rowspan="2">Statistical parameters</th><th colspan="3">σ'<sub>1</sub>/σ'<sub>3</sub></th></tr><tr><th>A</th><th>B</th><th>Average</th></tr></thead><tbody><tr><td rowspan="4">1</td><td>Mean</td><td>1.9</td><td>1.9</td><td>1.9</td></tr><tr><td>St dev</td><td>0.4</td><td>0.6</td><td>0.3</td></tr><tr><td>Median</td><td>1.9</td><td>1.9</td><td>1.9</td></tr><tr><td>[Min-Max]</td><td>[0.9-3.3]</td><td>[0.9-3.6]</td><td>[1.0-3.1]</td></tr><tr><td rowspan="4">2</td><td>Mean</td><td>1.6</td><td>1.6</td><td>1.6</td></tr><tr><td>St dev</td><td>0.3</td><td>0.5</td><td>0.3</td></tr><tr><td>Median</td><td>1.6</td><td>1.6</td><td>1.6</td></tr><tr><td>[Min-Max]</td><td>[0.7-2.7]</td><td>[0.7-3.0]</td><td>[0.8-2.6]</td></tr><tr><td rowspan="4">3</td><td>Mean</td><td>1.5</td><td>1.5</td><td>1.5</td></tr><tr><td>St dev</td><td>0.3</td><td>0.4</td><td>0.3</td></tr><tr><td>Median</td><td>1.5</td><td>1.5</td><td>1.5</td></tr><tr><td>[Min-Max]</td><td>[0.7-2.6]</td><td>[0.7-2.7]</td><td>[0.8-2.5]</td></tr><tr><td rowspan="4">4</td><td>Mean</td><td>1.5</td><td>1.5</td><td>1.5</td></tr><tr><td>St dev</td><td>0.3</td><td>0.4</td><td>0.3</td></tr><tr><td>Median</td><td>1.5</td><td>1.4</td><td>1.5</td></tr><tr><td>[Min-Max]</td><td>[0.7-2.6]</td><td>[0.6-2.7]</td><td>[0.7-2.4]</td></tr><tr><td rowspan="4">5</td><td>Mean</td><td>1.4</td><td>1.4</td><td>1.4</td></tr><tr><td>St dev</td><td>0.3</td><td>0.4</td><td>0.3</td></tr><tr><td>Median</td><td>1.4</td><td>1.4</td><td>1.4</td></tr><tr><td>[Min-Max]</td><td>[0.7-2.5]</td><td>[0.6-2.6]</td><td>[0.7-2.4]</td></tr></tbody></table>

σ'₁—effective maximum principal stress, σ'₃—effective minimum principal stress.

the hypotheses 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
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

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
<table><caption>**Table 7.** Possibilities for the maximum principal stress</caption><thead><tr><th rowspan="2">Depth (km)</th><th rowspan="2">Statistical parameters</th><th colspan="2">σ<sub>H</sub> (MPa)</th><th rowspan="2">Relative difference (%)</th></tr><tr><th>Estimated based on stress ratio coefficient (average)</th><th>Estimated based on stress polygon (average)</th></tr></thead><tbody><tr><td rowspan="4">1</td><td>Mean</td><td>30.0</td><td>82.9</td><td rowspan="4">49</td></tr><tr><td>St dev</td><td>3.0</td><td>14.0</td></tr><tr><td>Median</td><td>30.3</td><td>81.8</td></tr><tr><td>[Min-Max]</td><td>[21.9-38.1]</td><td>[54.3-126.0]</td></tr><tr><td rowspan="4">2</td><td>Mean</td><td>56.7</td><td>164.5</td><td rowspan="4">54</td></tr><tr><td>St dev</td><td>5.0</td><td>27.6</td></tr><tr><td>Median</td><td>57.0</td><td>162.5</td></tr><tr><td>[Min-Max]</td><td>[43.5-70.7]</td><td>[107.5-244.7]</td></tr><tr><td rowspan="4">3</td><td>Mean</td><td>83.4</td><td>247.9</td><td rowspan="4">57</td></tr><tr><td>St dev</td><td>6.8</td><td>41.4</td></tr><tr><td>Median</td><td>83.6</td><td>244.7</td></tr><tr><td>[Min-Max]</td><td>[64.5-102.5]</td><td>[160.1-367.9]</td></tr><tr><td rowspan="4">4</td><td>Mean</td><td>109.8</td><td>330.7</td><td rowspan="4">58</td></tr><tr><td>St dev</td><td>8.7</td><td>55.2</td></tr><tr><td>Median</td><td>109.9</td><td>326.5</td></tr><tr><td>[Min-Max]</td><td>[85.1-136.3]</td><td>[209.0-492.7]</td></tr><tr><td rowspan="4">5</td><td>Mean</td><td>135.9</td><td>412.3</td><td rowspan="4">59</td></tr><tr><td>St dev</td><td>10.7</td><td>68.7</td></tr><tr><td>Median</td><td>137.1</td><td>406.7</td></tr><tr><td>[Min-Max]</td><td>[106.7-168.3]</td><td>[270.1-614.5]</td></tr></tbody></table>

σ<sub>H</sub>—maximum horizontal stress.

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⁻¹ to meet the hydraulic impedance of 1 and 10 L·s⁻¹·MPa⁻¹. 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
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** 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.</figcaption></figure>

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
<table><caption>**Table 8.** Stimulation fluid pressure</caption><thead><tr><th rowspan="2">Depth (km)</th><th rowspan="2">Statistical parameters</th><th colspan="3">P<sub>p</sub> (MPa)</th><th rowspan="2">Relative difference (%)</th></tr><tr><th>Hydrostatic conditions (average)</th><th>Stimulation fluid pressure (average)</th><th>Overpressure conditions</th></tr></thead><tbody><tr><td rowspan="4">1</td><td>Mean</td><td>10.9</td><td>9.1</td><td>20.0</td><td rowspan="4">45</td></tr><tr><td>St dev</td><td>0.5</td><td>3.1</td><td>3.3</td></tr><tr><td>Median</td><td>10.9</td><td>9.1</td><td>20.0</td></tr><tr><td>[Min-Max]</td><td>[9.7-12.3]</td><td>[0-19.2]</td><td>[9.5-30.4]</td></tr><tr><td rowspan="4">2</td><td>Mean</td><td>21.7</td><td>23.0</td><td>44.7</td><td rowspan="4">51</td></tr><tr><td>St dev</td><td>0.9</td><td>5.6</td><td>6.0</td></tr><tr><td>Median</td><td>21.7</td><td>23.0</td><td>44.7</td></tr><tr><td>[Min-Max]</td><td>[19.4-24.3]</td><td>[4.3-42.4]</td><td>[24.4-65.7]</td></tr><tr><td rowspan="4">3</td><td>Mean</td><td>32.6</td><td>37.4</td><td>70.0</td><td rowspan="4">53</td></tr><tr><td>St dev</td><td>1.3</td><td>7.7</td><td>8.3</td></tr><tr><td>Median</td><td>32.6</td><td>37.4</td><td>70.0</td></tr><tr><td>[Min-Max]</td><td>[29.2-36.2]</td><td>[12.0-62.0]</td><td>[43.0-96.8]</td></tr><tr><td rowspan="4">4</td><td>Mean</td><td>43.5</td><td>50.7</td><td>94.2</td><td rowspan="4">54</td></tr><tr><td>St dev</td><td>1.8</td><td>10.6</td><td>11.5</td></tr><tr><td>Median</td><td>43.5</td><td>50.5</td><td>94.0</td></tr><tr><td>[Min-Max]</td><td>[39.0-48.4]</td><td>[17.4-86.3]</td><td>[57.3-133.9]</td></tr><tr><td rowspan="4">5</td><td>Mean</td><td>54.3</td><td>65.2</td><td>119.5</td><td rowspan="4">55</td></tr><tr><td>St dev</td><td>2.2</td><td>13.1</td><td>14.2</td></tr><tr><td>Median</td><td>54.3</td><td>65.3</td><td>119.6</td></tr><tr><td>[Min-Max]</td><td>[48.6-60.5]</td><td>[25.2-107.0]</td><td>[75.5-163.9]</td></tr></tbody></table>

$P_p$—in situ fluid pressure.

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 bore-holes 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 di-latation 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 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

<table><thead><tr><th>Symbol</th><th>Definition</th><th>Unit</th></tr></thead><tbody><tr><td><i>F</i></td><td>Factor</td><td>—</td></tr><tr><td>FBP</td><td>Formation breakdown pressure</td><td>Pa</td></tr><tr><td><i>g</i></td><td>Gravitational acceleration</td><td>m·s<sup>-2</sup></td></tr><tr><td><i>k</i></td><td>Stress ratio coefficient</td><td>—</td></tr><tr><td><i>P</i></td><td>Pressure</td><td>Pa</td></tr><tr><td><i>S</i></td><td>Slip tendency</td><td>—</td></tr></tbody></table>
<table><tr><td><i>T</i></td><td>Tensile strength</td><td>Pa</td></tr><tr><td><i>z</i></td><td>Depth</td><td>m</td></tr><tr><td colspan="3"><i>Greek letters</i></td></tr><tr><td>μ</td><td>Friction coefficient</td><td>—</td></tr><tr><td>ρ</td><td>Density</td><td>kg·m<sup>-3</sup></td></tr><tr><td>σ</td><td>Stress</td><td>Pa</td></tr><tr><td>σ'</td><td>Effective stress</td><td>Pa</td></tr><tr><td>τ</td><td>Shear stress</td><td>Pa</td></tr></table>

Subscript

<table><tr><td>1</td><td>Maximum principal</td></tr><tr><td>3</td><td>Minimum principal</td></tr><tr><td>H</td><td>Maximum horizontal</td></tr><tr><td>h</td><td>Minimum horizontal</td></tr><tr><td>max</td><td>Maximum</td></tr><tr><td>min</td><td>Minimum</td></tr><tr><td>n</td><td>Normal</td></tr><tr><td>p</td><td>In-situ fluid</td></tr><tr><td>p-f</td><td>Pore to fluid</td></tr><tr><td>V</td><td>Vertical</td></tr></table>

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.

References

Adams, J. (1989). Crustal stresses in Eastern Canada. In Gregersen, S. and Basham, P. W., editors, *Earthquakes at North-Atlantic Passive Margins: Neotectonics and Postglacial Rebound*, pages 289–298. Kluwer Academic Publishers, Dordrecht.

Allmendinger, R. W. (2020). MohrPlotter. http://www.geo.cornell.edu/geology/faculty/RWA/programs/mohrplotter.html. (Accessed November 19, 2020).

Amadei, B. and Stephansson, O. (1997). *Rock Stress and its Measurement*. Springer, Dordrecht.

Andrieux, P. P., Li, H., Brummer, R. K., Caumartin, R., and Rispoli, A. (2009). Rock mechanics approach for the recovery of the zone G crown pillar at the Raglan Katinniq Mine. In *Proceedings, 3rd CANUS Rock Mechanics Symposium, Toronto (ON), Canada, May 2009*. https://geogroup.utoronto.ca/wp-content/uploads/RockEng09/PDF/Session3/3965%20PAPER.pdf. (Accessed November 19, 2021).

Arjang, B. (1991). Pre-mining stresses at some hard rock mines in the Canadian Shield. *CIM Bull.*, 84(945), 80–86.

Arjang, B. (1998). Canadian crustal stresses and their application in mine design. In Singhal, R. K., editor, *Mine Planning and Equipment Selection 1998*, pages 269–274. A. A. Balkema, Rotterdam.

Arriaga, M., Brooks, M., and Moore, N. (2017). Energy access—the Canadian context. Waterloo Global Initiative's OpenAccess Energy Blueprint. http://wgsi.org/sites/wgsi-live.pi.local/files/EnergyAccess_Canadian_Context_Infographic_Spread-OpenAccess_EnergyWGSI_2017.pdf. (Accessed November 18, 2019).

Audet, P. and Mareschal, J.-C. (2004). Variations in elastic thickness in the Canadian Shield. *Earth Planet. Sci. Lett.*, 226, 17–31.

Beardsmore, G. R., Rybach, L., Blackwell, D., and Baron, C. (2010). A protocol for estimating and mapping global EGS potential. *GRC Trans.*, 34, 301–312.

Brown, D. W., Duchane, D. V., Heiken, G., and Hriscu, V. T. (2012). *Mining the Earth's Heat: Hot Dry Rock Geothermal Energy*. Springer, Berlin.

Burov, E. B. (2011). Rheology and strength of the lithosphere. *Mar. Petrol. Geol.*, 28, 1402–1443.

ECCC—Environment and Climate Change Canada (2019). Historical Climate Data. Canada: Government of Canada. Updated 2019-09-19. http://climate.weather.gc.ca/. Accessed 20 September 2019.

Evans, K., Cornet, F., Hayashi, K., Hashida, T., Ito, T., Matsuki, K., and Wallroth, T. (1999). Stress and rock
mechanics issues of relevance to HDR/HWR en-
gineered geothermal systems: review of develop-
ments during the past 15 years. *Geothermics*, 28(4-
5), 455-474.

Ghassemi, A. (2012). A review of some rock mechan-
ics issues in geothermal reservoir development.
*Geotech. Geol. Eng.*, 30, 647–664.

Gilmour, B., Oldfield, E., Platis, H., and Wicks, E. (2018). Toward a positive energy future in northern and remote communities. Summary Report. QUEST, Ottawa. 12 pages.

Grasby, S. E., Allen, D. M., Bell, S., Chen, Z., Fer-
guson, G., Jessop, A., Kelman, M., Ko, M., Ma-
jorowicz, J., Moore, M., Raymond, J., and Ther-
rien, R. (2012). Geothermal energy resource po-
tential of Canada. Report No.: Open File 6914.
Geological Survey of Canada, Ottawa, 322 pages.
http://ftp.geogratis.gc.ca/pub/nrcan_rncan/
publications/ess_sst/291/291488/of_6914.pdf
(Accessed July 23, 2019).

Grasby, S. E., Majorowicz, J., and McCune, G. (2013). Geothermal energy potential for northern communities. Report No.: Open File 7350. Geological Survey of Canada, Ottawa, 50 pages. http://ftp.geogratis.gc.ca/pub/nrcan_rncan/publications/ess_sst/292/292840/of_7350.pdf (Accessed July 23, 2019).

Hall, J., Louden, K. E., Funck, T., and Deemer, S. (2002). Geophysical characteristics of the continental crust along the Lithoprobe Eastern Canadian Shield Onshore-Offshore Transect (ECSOOT): a review. *Can. J. Earth Sci.*, 39, 569–587.

Hammer, P. T. C., Clowes, R. M., Cook, F. A., van der Velden, A. J., and Vasudevan, K. (2010). The Lithoprobe trans-continental lithospheric cross sections: imaging the internal structure of the North American continent. *Can. J. Earth Sci.*, 47, 821–857.

Heidbach, O., Rajabi, M., Reiter, K., and Ziegler, M. O. (2019). World stress map. In Sorkhabi, R., editor, *Encyclopedia of Petroleum Geoscience*. Springer, Cham.

Heidbach, O., Tingay, M., Barth, A., Reinecker, J., Kur-
feß, D., and Müller, B. (2010). Global crustal stress
pattern based on the World Stress Map database re-
lease 2008. *Tectonophysics*, 482, 3–15.

Herget, G. (1982). High stress occurrences in the
Canadian Shield. In *23rd US Symposium on*

*Rock Mechanics (USRMS)*, Berkeley, US, August 25,
1982.

Herget, G. (1987). Stress assumptions for under-
ground excavations in the Canadian shield. *Int. J.
Rock Mech. Min. Sci. Geomech. Abs.*, 24(1), 95-97.

Herget, G. (1993). Rock stresses and rock stress mon-
itoring in Canada. In Hudson, J. A., editor, *Rock
Testing and Site Characterization: Principles, Prac-
tice and Projects*, pages 473–496. Elsevier, Berlin.

Hoek, E. and Brown, E. T. (1980). *Underground Ex-
cavation in Rock*. Institute of Mining & Metallurgy,
London.

Hubbert, M. K. and Willis, D. G. (1957). Mechanics of hydraulic fracturing. *Trans. Soc. Petrol. Eng. AIME*, 210, 153–163.

Hudson, J. A., Cornet, F. H., and Christiansson, R. (2003). ISRM suggested methods for rock stress estimation—part 1: strategy for rock stress estimation. *Int. J. Rock Mech. Min. Sci.*, 40, 991–998.

Jaeger, J. C., Cook, N. G. W., and Zimmerman, R. W. (2007). *Fundamentals of Rock Mechanics*. Blackwell Publishing, Oxford.

Kaiser, P. K., Maloney, S. M., and Yong, S. (2016). Role of large scale heterogeneities on in-situ stress and induced stress fields. In *Proceedings, 50th US Rock Mechanics/Geomechanics Symposium, Houston, Texas, USA, 26-29 June 2016*. American Rock Mechanics Association.

KRG—Kativik Regional Government, Makivik Corpo-
ration (2010). Plan Nunavik. Avataq Cultural Insti-
tute, Montreal.

Kruszewski, M., Hofmann, H., Alvarez, F. G., Bianco,
C., Haro, A. J., Garduño, V. H., Liotta, D., Trumpy, E.,
Brogi, A., Wheeler, W., Bastesen, E., Parisio, F., and
Saenger, E. H. (2021). Integrated stress field esti-
mation and implications for enhanced geothermal
system development in Acoculco, Mexico. *Geother-
mics*, 89, article no. 101931.

Miranda, M. M. (2021). *Assessment of the Deep Geothermal Energy Source Potential in Remote Northern Regions: A Study Undertaken in the Sub-arctic Off-grid Community of Kuujjuaq, Nunavik, Canada*. PhD thesis, Institut national de la recherche scientifique. 438 pages. PhD in Earth Sciences.

Miranda, M. M., Dezayes, C., Giordano, N., Kan-
zari, I., Raymond, J., and Carvalho, J. R. B. (2018).
Fracture network characterization as input for
geothermal energy research: preliminary data from Kuujjuaq, northern Québec, Canada. In *Proceedings, 43rd Workshop on Geothermal Reservoir Engineering, Stanford (CA), US, February 12–14, 2018*. https://pangea.stanford.edu/ERE/pdf/IGAstandard/SGW/2018/Miranda.pdf. (Accessed November 19, 2021).

Miranda, M. M., Giordano, N., Raymond, J., Pereira, A. J. S. C., and Dezayes, C. (2020a). Thermophysical properties of surficial rocks: a tool to characterize geothermal resources of remote northern regions. *Geotherm. Energy*, 8, article no. 4.

Miranda, M. M., Raymond, J., and Dezayes, C. (2020b). Uncertainty and risk evaluation of deep geothermal energy source for heat production and electricity generation in remote northern regions. *Energies*, 13(16), article no. 4221.

Miranda, M. M., Raymond, J., Willis-Richards, J., and Dezayes, C. (2021a). Are engineered geothermal energy systems a viable solution for arctic off-grid communities? A techno-economic study. *Water*, 13, article no. 3526.

Miranda, M. M., Velez Márquez, M. I., Raymond, J., and Dezayes, C. (2021b). A numerical approach to infer terrestrial heat flux from shallow temperature profiles in remote northern regions. *Geothermics*, 93, article no. 102064.

Moeck, I., Kwiatek, G., and Zimmermann, G. (2009). Slip tendency analysis, fault reactivation potential and induced seismicity in a deep geothermal reservoir. *J. Struct. Geol.*, 31, 1174–1182.

Morris, A., Ferrill, D. A., and Henderson, D. B. (1996). Slip-tendency analysis and fault reactivation. *Geology*, 24(3), 275–278.

Nakatsuka, K. (1999). Field characterization for HDR/HWR: a review. *Geothermics*, 28(4–5), 519–531.

NRCan (2018). Earthquake hazard. https://earthquakescanada.nrcan.gc.ca/hazard-alea/index-en.php. (Accessed November 23, 2022).

Palisade (2019). @RISK. https://www.palisade.com/risk/default.asp. (Accessed June 18, 2020).

Peacock, D. C. P., Sanderson, D. J., and Leiss, B. (2021). Use of Mohr diagrams to predict fracturing in a potential geothermal reservoir. *Geosciences*, 11, article no. 501.

Pine, R. J. and Batchelor, A. S. (1984). Downward migration of shearing in jointed rock during hydraulic

injections. *Int. J. Rock Mech. Min. Sci. Geomech. Abs.*, 21(5), 249–263.

Richards, H. G., Parker, R. H., Green, A. S. P., Jones, R. H., Nicholls, J. D. M., Nicol, D. A. C., Randall, M. M., Richards, S., Stewart, R. C., and Willis-Richards, J. (1994). The performance and characteristics of the experimental hot dry rock geothermal reservoir at Rosemanowes, Cornwall (1985–1988). *Geothermics*, 23(2), 73–109.

Rutqvist, J., Birkholzer, J., Cappa, F., and Tsang, C.-F. (2007). Estimating maximum sustainable injection pressure during geological sequestration of CO2 using coupled fluid flow and geomechanical fault-slip analysis. *Energy Convers. Manage.*, 48(6), 1798–1807.

Sibson, R. H. (2004). Controls on maximum fluid overpressure defining conditions for mesozonal mineralization. *J. Struct. Geol.*, 26(6–7), 1127–1136.

SIGÉOM—Système d'information géominière du Québec (2019). Carte interactive. http://sigeom.mines.gouv.qc.ca/signet/classes/I1108_afchCarteIntr. (Accessed May 5, 2019).

Simard, M., Lafrance, I., Hammouche, H., and Legouix, C. (2013). Géologie de la Région de Kuujjuaq et de la Baie d'Ungava (SRNC 24J, 24K). Report No.: RG 2013-04. Gouvernement du Québec, Québec. 62 pages. http://gq.mines.gouv.qc.ca/documents/EXAMINE/RG201304/RG201304.pdf (Accessed May 5, 2019).

Taghipour, M., Ghafoori, M., Lashkaripour, G. R., Moghadas, N. H., and Molaghab, A. (2019). Estimation of the current stress field and fault reactivation analysis in the Asmari reservoir, SW Iran. *Petrol. Sci.*, 16, 513–526.

Taherynia, M. H., Aghda, S. M. F., and Fahimifar, A. (2016). In-situ stress state and tectonic regime in different depths of Earth crust. *Geotech. Geol. Eng.*, 34, 679–687.

Telmat, H., Mareschal, J.-C., and Gariépy, C. (1999). The gravity field over the Ungava Bay region from satellite altimetry and new land-based data: implications for the geology of the area. *Can. J. Earth Sci.*, 36, 75–89.

Townend, J. and Zoback, M. D. (2000). How faulting keeps the crust strong. *Geology*, 28(5), 399–402.

van den Bogert, P. A. J. and van Eijs, R. M. H. E. (2020). Why Mohr-circle analyses may underestimate the risk of fault reactivation in depleting
reservoirs. *Int. J. Rock Mech. Min. Sci.*, 136, article no. 104502.

Vervaet, F. and Darbyshire, F. A. (2016). Moho depth and bulk crustal properties in northern Quebec and Labrador. In *American Geophysical Union, Fall Meeting 2016*.

Wardle, R. J., James, D. T., Scott, D. J., and Hall, J. (2002). The southeastern Churchill Province: synthesis of a Paleoproterozoic transpressional orogen. *Can. J. Earth Sci.*, 39, 639–663.

Whitmeyer, S. J. and Karlstrom, K. E. (2007). Tectonic model for the Proterozoic growth of North America. *Geosphere*, 3(4), 220–259.

Xie, L., Min, K.-B., and Song, Y. (2015). Observations of hydraulic stimulations in seven enhanced geothermal system projects. *Renew. Energy*, 79, 56–65.

Young, S. and Maloney, S. (2015). An update to the Canadian Shield stress database. NWMO-TR-2015-18. NWMO, Toronto, 33 pages. https://www.nwmo.ca/~/media/Site/Reports/2016/

08/23/10/52/NWMO_TR_2015_18.ashx?la=en
(Accessed Feb 26, 2021).

Zhang, L. (2017). *Engineering Properties of Rocks*. Elsevier, Oxford. 394 pages.

Zoback, M. D. (2007). *Reservoir Geomechanics*. Cambridge University Press, Cambridge.

Zoback, M. D. and Townend, J. (2001). Implications of hydrostatic pore pressures and high crustal strength for the deformation of intraplate lithosphere. *Tectonophysics*, 336(1–4), 19–30.

Zoback, M. L., Zoback, M. D., Adams, J., Assumpção, M., Bell, S., Bergman, E. A., Blümling, P., Brereton, N. R., Denham, D., Ding, J., Fuchs, K., Gay, N., Gregersen, S., Gupta, H. K., Gvishiani, A., Jacob, K., Klein, R., Knoll, P., Magee, M., Mercier, J. L., Müller, B. C., Paquin, C., Rajendran, K., Stephanson, O., Suarez, G., Suter, M., Udias, A., Xu, Z. H., and Zhizhin, M. (1989). Global patterns of tectonic stress. *Nature*, 341, 291–298.