Plan
Comptes Rendus

Tectonique, tectonophysique, géodynamique
Estimating theoretical stress regime for engineered geothermal energy systems in an arctic community (Kuujjuaq, Canada)
Comptes Rendus. Géoscience, Volume 355 (2023), pp. 85-108.

Résumé

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.

Métadonnées
Reçu le :
Révisé le :
Accepté le :
Publié le :
DOI : 10.5802/crgeos.193
Mots clés : Stress regime, Mohr–Coulomb criterion, Slip tendency, Fluid overpressure, Geothermal energy, Nunavik

Mafalda M. Miranda 1 ; Jasmin Raymond 1 ; Chrystel Dezayes 2

1 INRS—Institute national de la recherche scientifique, 490, rue de la Couronne Québec (QC) G1K 9A9, Canada
2 BRGM, 3 avenue Claude Guillemin, BP 36009, 45060, Orléans Cedex 2, France
Licence : CC-BY 4.0
Droits d'auteur : Les auteurs conservent leurs droits
@article{CRGEOS_2023__355_G1_85_0,
     author = {Mafalda M. Miranda and Jasmin Raymond and Chrystel Dezayes},
     title = {Estimating theoretical stress regime for engineered geothermal energy systems in an arctic community {(Kuujjuaq,} {Canada)}},
     journal = {Comptes Rendus. G\'eoscience},
     pages = {85--108},
     publisher = {Acad\'emie des sciences, Paris},
     volume = {355},
     year = {2023},
     doi = {10.5802/crgeos.193},
     language = {en},
}
TY  - JOUR
AU  - Mafalda M. Miranda
AU  - Jasmin Raymond
AU  - Chrystel Dezayes
TI  - Estimating theoretical stress regime for engineered geothermal energy systems in an arctic community (Kuujjuaq, Canada)
JO  - Comptes Rendus. Géoscience
PY  - 2023
SP  - 85
EP  - 108
VL  - 355
PB  - Académie des sciences, Paris
DO  - 10.5802/crgeos.193
LA  - en
ID  - CRGEOS_2023__355_G1_85_0
ER  - 
%0 Journal Article
%A Mafalda M. Miranda
%A Jasmin Raymond
%A Chrystel Dezayes
%T Estimating theoretical stress regime for engineered geothermal energy systems in an arctic community (Kuujjuaq, Canada)
%J Comptes Rendus. Géoscience
%D 2023
%P 85-108
%V 355
%I Académie des sciences, Paris
%R 10.5802/crgeos.193
%G en
%F CRGEOS_2023__355_G1_85_0
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, Volume 355 (2023), pp. 85-108. doi : 10.5802/crgeos.193. https://comptes-rendus.academie-sciences.fr/geoscience/articles/10.5802/crgeos.193/

Version originale du texte intégral (Proposez une traduction )

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.

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

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

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.

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

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

Table 1.

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

𝜎 V = 𝜌 g z (1)
where 𝜎V (Pa) is the vertical stress, 𝜌 (kg⋅m−3) is the density of the geological materials, g (m⋅s−2) 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]:

𝜎 H = k max 𝜎 V 𝜎 h = k min 𝜎 V (2)
where k (–) is the stress ratio coefficient and 𝜎H and 𝜎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.

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

R = 𝜎 2 𝜎 3 𝜎 1 𝜎 3 R = 𝜎 h 𝜎 V 𝜎 H 𝜎 V (3)
where 𝜎1 (Pa) is the maximum principal stress, 𝜎2 (Pa) is the intermediate principal stress and 𝜎3 (Pa) is the 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 pf 𝜎 V (4)
where Pp (Pa) is the in situ fluid pressure and Fp−f 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]:

𝜎 1 𝜎 3 𝜇 2 + 1 + 𝜇 2 (5)
where 𝜎 1 and 𝜎 3 (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 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:

𝜎 1 = 𝜎 H P p 𝜎 3 = 𝜎 h P p (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]:

f ( 𝜇 ) = 𝜇 2 + 1 + 𝜇 2 𝜎 H = 𝜎 h ( 𝜎 V P p ) f ( 𝜇 ) + P p = 𝜎 h normal faulting regime ( 𝜎 H P p ) ( 𝜎 h P p ) = f ( 𝜇 ) strike-slip faulting regime f ( 𝜇 ) ( 𝜎 V P p ) + P p = 𝜎 H reverse faulting regime.   (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.

Table 3.

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

FBP= 3 𝜎 h 𝜎 H P p + T (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 = 3 𝜎 3 𝜎 1 2 . (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

Table 4.

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

S = 𝜏 𝜎 n (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).

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:

𝜎 V ( MPa ) = ( 0 . 0 2 7 ± 0 . 0 3 ) z 𝜎 H ( MPa ) = ( 9 . 5 ± 3 . 5 ) + ( 0 . 0 3 3 ± 0 . 0 1 2 ) z 𝜎 h ( MPa ) = ( 2 . 5 ± 2 . 5 ) + ( 0 . 0 2 7 ± 0 . 0 6 ) z P p ( MPa ) = ( 0 . 0 1 1 ± 0 . 0 1 ) z . (11)
The estimated vertical stress suggests a gradient of 27 ± 3 MPa⋅km−1 and the gradient for the maximum and minimum horizontal stress was evaluated as 33 ± 12 MPa⋅km−1 and 27 ± 6 MPa⋅km−1, respectively. The gradient for the in situ fluid pressure was estimated as 11 ± 0.1 MPa⋅km−1. 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 correspond 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.

Table 5.

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.

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.

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.

Figure 6.

Effective stress ratio for a depth of 4 km as an example.

Table 6.

Estimated effective stress ratio coefficient

Depth (km) Statistical parameters 𝜎 1 𝜎 3
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]

𝜎 1 —effective maximum principal stress, 𝜎 3 —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).

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.

Table 7.

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

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.

Table 8.

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

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.

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.


Bibliographie

[Adams, 1989] J. Adams Crustal stresses in Eastern Canada, Earthquakes at North-Atlantic Passive Margins: Neotectonics and Postglacial Rebound (S. Gregersen; P. W. Basham, eds.), Kluwer Academic Publishers, Dordrecht, 1989, pp. 289-298

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

[Amadei and Stephansson, 1997] B. Amadei; O. Stephansson Rock Stress and its Measurement, Springer, Dordrecht, 1997 | DOI

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

[Arjang, 1991] B. Arjang Pre-mining stresses at some hard rock mines in the Canadian Shield, CIM Bull., Volume 84 (1991) no. 945, pp. 80-86

[Arjang, 1998] B. Arjang Canadian crustal stresses and their application in mine design, Mine Planning and Equipment Selection 1998 (R. K. Singhal, ed.), A. A. Balkema, Rotterdam, 1998, pp. 269-274

[Arriaga et al., 2017] M. Arriaga; M. Brooks; N. Moore Energy access—the Canadian context, 2017 Waterloo Global Initiative’s OpenAcess Energy Blueprint. http://wgsi.org/sites/wgsi-live.pi.local/files/Energy_Access_Canadian_Context_Infographic_Spread-OpenAccess_Energy_Blueprint_WGSI_2017.pdf (Accessed November 18, 2019)

[Audet and Mareschal, 2004] P. Audet; J.-C. Mareschal Variations in elastic thickness in the Canadian Shield, Earth Planet. Sci. Lett., Volume 226 (2004), pp. 17-31 | DOI

[Beardsmore et al., 2010] G. R. Beardsmore; L. Rybach; D. Blackwell; C. Baron A protocol for estimating and mapping global EGS potential, GRC Trans., Volume 34 (2010), pp. 301-312

[Brown et al., 2012] D. W. Brown; D. V. Duchane; G. Heiken; V. T. Hriscu Mining the Earth’s Heat: Hot Dry Rock Geothermal Energy, Springer, Berlin, 2012

[Burov, 2011] E. B. Burov Rheology and strength of the lithosphere, Mar. Petrol. Geol., Volume 28 (2011), pp. 1402-1443 | DOI

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

[Evans et al., 1999] K. Evans; F. Cornet; K. Hayashi; T. Hashida; T. Ito; K. Matsuki; T. Wallroth Stress and rock mechanics issues of relevance to HDR/HWR engineered geothermal systems: review of developments during the past 15 years, Geothermics, Volume 28 (1999) no. 4–5, pp. 455-474 | DOI

[Ghassemi, 2012] A. Ghassemi A review of some rock mechanics issues in geothermal reservoir development, Geotech. Geol. Eng., Volume 30 (2012), pp. 647-664 | DOI

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

[Grasby et al., 2012] S. E. Grasby; D. M. Allen; S. Bell; Z. Chen; G. Ferguson; A. Jessop; M. Kelman; M. Ko; J. Majorowicz; M. Moore; J. Raymond; R. Therrien Geothermal energy resource potential of Canada, 2012 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 et al., 2013] S. E. Grasby; J. Majorowicz; G. McCune Geothermal energy potential for northern communities, 2013 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 et al., 2002] J. Hall; K. E. Louden; T. Funck; S. Deemer Geophysical characteristics of the continental crust along the Lithoprobe Eastern Canadian Shield Onshore-Offshore Transect (ECSOOT): a review, Can. J. Earth Sci., Volume 39 (2002), pp. 569-587 | DOI

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

[Heidbach et al., 2010] O. Heidbach; M. Tingay; A. Barth; J. Reinecker; D. Kurfeß; B. Müller Global crustal stress pattern based on the World Stress Map database release 2008, Tectonophysics, Volume 482 (2010), pp. 3-15 | DOI

[Heidbach et al., 2019] O. Heidbach; M. Rajabi; K. Reiter; M. O. Ziegler World stress map, Encyclopedia of Petroleum Geoscience (R. Sorkhabi, ed.), Springer, Cham, 2019

[Herget, 1982] G. Herget High stress occurrences in the Canadian Shield, 23rd US Symposium on Rock Mechanics (USRMS), Berkeley, US, August 25, 1982 (1982)

[Herget, 1987] G. Herget Stress assumptions for underground excavations in the Canadian shield, Int. J. Rock Mech. Min. Sci. Geomech. Abs., Volume 24 (1987) no. 1, pp. 95-97 | DOI

[Herget, 1993] G. Herget Rock stresses and rock stress monitoring in Canada, Rock Testing and Site Characterization: Principles, Practice and Projects (J. A. Hudson, ed.), Elsevier, Berlin, 1993, pp. 473-496 | DOI

[Hoek and Brown, 1980] E. Hoek; E. T. Brown Underground Excavation in Rock, Institute of Mining & Metallurgy, London, 1980

[Hubbert and Willis, 1957] M. K. Hubbert; D. G. Willis Mechanics of hydraulic fracturing, Trans. Soc. Petrol. Eng. AIME, Volume 210 (1957), pp. 153-163 | DOI

[Hudson et al., 2003] J. A. Hudson; F. H. Cornet; R. Christiansson ISRM suggested methods for rock stress estimation—part 1: strategy for rock stress estimation, Int. J. Rock Mech. Min. Sci., Volume 40 (2003), pp. 991-998 | DOI

[Jaeger et al., 2007] J. C. Jaeger; N. G. W. Cook; R. W. Zimmerman Fundamentals of Rock Mechanics, Blackwell Publishing, Oxford, 2007

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

[KRG—Kativik Regional Government, Makivik Corporation, 2010] KRG—Kativik Regional Government, Makivik Corporation, 2010 (Plan Nunavik. Avataq Cultural Institute, Montreal)

[Kruszewski et al., 2021] M. Kruszewski; H. Hofmann; F. G. Alvarez; C. Bianco; A. J. Haro; V. H. Garduño; D. Liotta; E. Trumpy; A. Brogi; W. Wheeler; E. Bastesen; F. Parisio; E. H. Saenger Integrated stress field estimation and implications for enhanced geothermal system development in Acoculco, Mexico, Geothermics, Volume 89 (2021), 101931 | DOI

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

[Miranda et al., 2020a] M. M. Miranda; N. Giordano; J. Raymond; A. J. S. C. Pereira; C. Dezayes Thermophysical properties of surficial rocks: a tool to characterize geothermal resources of remote northern regions, Geotherm. Energy, Volume 8 (2020), 4 | DOI

[Miranda et al., 2020b] M. M. Miranda; J. Raymond; C. Dezayes Uncertainty and risk evaluation of deep geothermal energy source for heat production and electricity generation in remote northern regions, Energies, Volume 13 (2020) no. 16, 4221 | DOI

[Miranda et al., 2021a] M. M. Miranda; J. Raymond; J. Willis-Richards; C. Dezayes Are engineered geothermal energy systems a viable solution for arctic off-grid communities? A techno-economic study, Water, Volume 13 (2021), 3526 | DOI

[Miranda et al., 2021b] M. M. Miranda; M. I. Velez Márquez; J. Raymond; C. Dezayes A numerical approach to infer terrestrial heat flux from shallow temperature profiles in remote northern regions, Geothermics, Volume 93 (2021), 102064 | DOI

[Miranda, 2021] M. M. Miranda 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, Ph. D. Thesis, Institut national de la recherche scientifique (2021) (438 pages. PhD in Earth Sciences)

[Moeck et al., 2009] I. Moeck; G. Kwiatek; G. Zimmermann Slip tendency analysis, fault reactivation potential and induced seismicity in a deep geothermal reservoir, J. Struct. Geol., Volume 31 (2009), pp. 1174-1182 | DOI

[Morris et al., 1996] A. Morris; D. A. Ferrill; D. B. Henderson Slip-tendency analysis and fault reactivation, Geology, Volume 24 (1996) no. 3, pp. 275-278 | DOI

[Nakatsuka, 1999] K. Nakatsuka Field characterization for HDR/HWR: a review, Geothermics, Volume 28 (1999) no. 4–5, pp. 519-531 | DOI

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

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

[Peacock et al., 2021] D. C. P. Peacock; D. J. Sanderson; B. Leiss Use of Mohr diagrams to predict fracturing in a potential geothermal reservoir, Geosciences, Volume 11 (2021), 501

[Pine and Batchelor, 1984] R. J. Pine; A. S. Batchelor Downward migration of shearing in jointed rock during hydraulic injections, Int. J. Rock Mech. Min. Sci. Geomech. Abs., Volume 21 (1984) no. 5, pp. 249-263 | DOI

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

[Rutqvist et al., 2007] J. Rutqvist; J. Birkholzer; F. Cappa; C.-F. Tsang Estimating maximum sustainable injection pressure during geological sequestration of CO2 using coupled fluid flow and geomechanical fault-slip analysis, Energy Convers. Manage., Volume 48 (2007) no. 6, pp. 1798-1807 | DOI

[Sibson, 2004] R. H. Sibson Controls on maximum fluid overpressure defining conditions for mesozonal mineralization, J. Struct. Geol., Volume 26 (2004) no. 6–7, pp. 1127-1136 | DOI

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

[Simard et al., 2013] M. Simard; I. Lafrance; H. Hammouche; C. Legouix Géologie de la Région de Kuujjuaq et de la Baie d’Ungava (SRNC 24J, 24K), 2013 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 et al., 2019] M. Taghipour; M. Ghafoori; G. R. Lashkaripour; N. H. Moghadas; A. Molaghab Estimation of the current stress field and fault reactivation analysis in the Asmari reservoir, SW Iran, Petrol. Sci., Volume 16 (2019), pp. 513-526 | DOI

[Taherynia et al., 2016] M. H. Taherynia; S. M. F. Aghda; A. Fahimifar In-situ stress state and tectonic regime in different depths of Earth crust, Geotech. Geol. Eng., Volume 34 (2016), pp. 679-687 | DOI

[Telmat et al., 1999] H. Telmat; J.-C. Mareschal; C. Gariépy 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., Volume 36 (1999), pp. 75-89 | DOI

[Townend and Zoback, 2000] J. Townend; M. D. Zoback How faulting keeps the crust strong, Geology, Volume 28 (2000) no. 5, pp. 399-402 | DOI

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

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

[Wardle et al., 2002] R. J. Wardle; D. T. James; D. J. Scott; J. Hall The southeastern Churchill Province: sysnthesis of a Paleoproterozoic transpressional orogen, Can. J. Earth Sci., Volume 39 (2002), pp. 639-663 | DOI

[Whitmeyer and Karlstrom, 2007] S. J. Whitmeyer; K. E. Karlstrom Tectonic model for the Proterozoic growth of North America, Geosphere, Volume 3 (2007) no. 4, pp. 220-259 | DOI

[Xie et al., 2015] L. Xie; K.-B. Min; Y. Song Observations of hydraulic stimulations in seven enhanced geothermal system projects, Renew. Energy, Volume 79 (2015), pp. 56-65 | DOI

[Young and Maloney, 2015] S. Young; S. Maloney An update to the Canadian Shield stress database, 2015 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, 2017] L. Zhang Engineering Properties of Rocks, Elsevier, Oxford, 2017 (394 pages)

[Zoback and Townend, 2001] M. D. Zoback; J. Townend Implications of hydrostatic pore pressures and high crustal strength for the deformation of intraplate lithosphere, Tectonophysics, Volume 336 (2001) no. 1–4, pp. 19-30 | DOI

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

[Zoback, 2007] M. D. Zoback Reservoir Geomechanics, Cambridge University Press, Cambridge, 2007 | DOI


Commentaires - Politique