Next Article in Journal
Application of Machine Learning Techniques in Rainfall–Runoff Modelling of the Soan River Basin, Pakistan
Previous Article in Journal
Analysis of Long-Term Shoreline Observations in the Vicinity of Coastal Structures: A Case Study of South Bali Beaches
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Are Engineered Geothermal Energy Systems a Viable Solution for Arctic Off-Grid Communities? A Techno-Economic Study

1
INRS—Institut National de la Recherche Scientifique, 490 Rue de la Couronne, Québec, QC G1K 9A9, Canada
2
Fracture and Reliability Research Institute, Tohoku University, 2 Chome-1-1 Katahira, Aoba Ward, Sendai 980-8577, Japan
3
BRGM—Bureau de Recherches Géologiques et Minières, 45060 Orléans, France
*
Author to whom correspondence should be addressed.
Water 2021, 13(24), 3526; https://doi.org/10.3390/w13243526
Submission received: 21 October 2021 / Revised: 2 December 2021 / Accepted: 6 December 2021 / Published: 9 December 2021
(This article belongs to the Special Issue Subsurface Hydrothermal Modeling in the Arctic)

Abstract

:
Deep geothermal energy sources harvested by circulating fluids in engineered geothermal energy systems can be a solution for diesel-based northern Canadian communities. However, poor knowledge of relevant geology and thermo-hydro-mechanical data introduces significant uncertainty in numerical simulations. Here, a first-order assessment was undertaken following a “what-if” approach to help design an engineered geothermal energy system for each of the uncertain scenarios. Each possibility meets the thermal energy needs of the community, keeping the water losses, the reservoir flow impedance and the thermal drawdown within predefined targets. Additionally, the levelized cost of energy was evaluated using the Monte Carlo method to deal with the uncertainty of the inputs and assess their influence on the output response. Hydraulically stimulated geothermal reservoirs of potential commercial interest were simulated in this work. In fact, the probability of providing heating energy at a lower cost than the business-as-usual scenario with oil furnaces ranges between 8 and 92%. Although the results of this work are speculative and subject to uncertainty, geothermal energy seems a potentially viable alternative solution to help in the energy transition of remote northern communities.

1. Introduction

The lack of an environmentally benign energy supply for electricity, space heating, domestic hot water and cooking is still a reality, not only in the off-grid communities of Canada, but worldwide [1]. Fossil fuels have been the main sources of electricity, space heating and cooking fuels in the developed world. However, current climate change concerns, people’s well-being and the environment are changing this predominance. Clean energy supply and sustainability are buzzwords nowadays. The worldwide consumption of renewable energies has increased from 40 EJ in 1990 to 64 EJ in 2017 [1]. Although smaller than the consumption of non-renewables (304 EJ in 2017; [1]), this is an important increase. Off-grid renewable energy solutions, including stand-alone systems and microgrids, are often a viable electrification solution [2].
In Canada, for example, from the approximately 280 off-grid communities, 239 rely exclusively on diesel for electricity, space heating and domestic hot water [3]. Within the diesel-based settlements without road access, the fuel must be imported during summer and stored for year-round use [4]. Such an energy situation entails significant costs, low energy security and a high probability of damaging an already vulnerable ecosystem. Therefore, interest in assessing the potential for renewable sources to feed microgrids in remote communities has increased and several studies have been conducted (e.g., [5,6,7,8,9,10,11,12,13,14,15,16,17]). Among these options, deep geothermal energy sources can play a key role to provide baseload power and/or heat to the off-grid settlements (e.g., [18,19,20,21,22,23,24,25,26,27]). In fact, a first-order community-scale geothermal assessment undertaken in Kuujjuaq (Nunavik, Canada) suggested that the deep geothermal energy source can fulfil the community’s annual average heating demand of 37 GWh [28]. This community is settled on the Canadian Shield, a physiographic region that has been considered a target for geothermal exploration through engineered/enhanced geothermal systems [25,29]. The feasibility of such systems has been studied in the Western Canadian Sedimentary Basin and Arctic Lands [18,19,25,30,31,32], but few studies have been conducted in the Canadian Shield [24,25], where there are hundreds of off-grid communities relying heavily on fossil fuels [4,21]. Therefore, there is a need for a comprehensive estimation of the possible performance of deep engineered geothermal energy systems across the Canadian Shield. In this context, work has been carried out in the community of Kuujjuaq (Nunavik, Canada) using the limited local surficial and regional data available to provide first-order answers to the following key questions:
  • Will the hydraulic stimulation technique applied in crystalline basement rocks elsewhere develop a well-connected flowing system in Kuujjuaq? How can this be done? What further local geological and thermo-hydro-mechanical data is required for more accurate predictions?
  • Are the deep geothermal energy sources harvested by engineered geothermal energy systems in Kuujjuaq likely to be cost-competitive compared to fossil fuels?
The answers to these questions are, however, subject to high uncertainties due to the current poor knowledge of both geology and thermo-hydro-mechanical data. Unfortunately, no hydraulic stimulation field experiments have been carried out to date in Kuujjuaq, nor were they within the scope of the present study. Thus, no history matching is available to calibrate the numerical simulations carried out in this work. Nevertheless, a “what-if” approach was used to provide a range of possibilities and help to design an engineered geothermal energy system for each of the uncertain scenarios. Each possibility aims to provide the thermal energy needed for the community, while keeping the circulating water losses below 20%, the reservoir flow impedance below 1.0 MPa L−1 s−1 and the system thermal drawdown below 1 °C/year. Thus, this study offers a first-order prediction of the performance of engineered geothermal energy systems as off-grid solutions (and constraints on the geomechanical and geological assumptions required) to help in the energy transition of remote northern communities. Additionally, a preliminary evaluation of the levelized cost of energy was undertaken to forecast the economic potential of engineered geothermal energy systems in remote northern regions. Overall, this study may help trigger interest for further geothermal exploration which is fundamental for an accurate evaluation of the deep geothermal energy potential beneath off-grid settlements.
Thus, this study was motivated by the lack of clean energy supply in the majority of the Canadian remote northern communities. The goal was to assess if deep geothermal energy harvested by engineered geothermal energy systems, or enhanced geothermal systems, is a technically and economically viable alternative solution to offset the diesel consumption in such communities. The study was undertaken in the off-grid settlement of Kuujjuaq (Nunavik, Canada) to provide an example for the remaining off-grid communities. The study here described, and the results obtained, although highly speculative due to the lack of deep geothermal exploratory boreholes in the study area, represent an important contribution to understand the potential that deep geothermal energy sources have to offer in the energy transition of diesel-based regions of northern Canada and other arctic regions. Furthermore, this study aims to predict the performance of enhanced geothermal systems in a location with significant geothermal data gaps. This may raise awareness about the potential of geothermal energy in areas considered at first sight unviable and trigger the interest for further geothermal developments. The parameters that require further data collection and possible operational strategies to develop engineered geothermal energy systems for communities on the old Canadian Shield are highlighted in this manuscript.

2. Background Information

2.1. Engineered/Enhanced Geothermal Systems

Enhancing well productivity and the permeability of subsurface rocks has been common practice in the oil and gas industry for many years (e.g., [33]). However, it was only introduced in geothermal energy exploration in 1973. The first hydraulic stimulation experiments in crystalline rocks were done at Fenton Hill (e.g., [34]). The success of these field experiments opened new opportunities to explore geothermal energy sources in areas that were previously considered unviable (e.g., crystalline rocks with very low permeability). In fact, this new concept of recovering the Earth’s heat via a pressurized closed-loop circulation of fluid from the surface through a hydraulically stimulated and confined reservoir several kilometers deep made in crystalline basement rocks marked an important conceptual turning point in the geothermal energy industry.
Since the 1970s, several projects have been started worldwide, applying different stimulation techniques and in geological contexts ranging from crystalline to sedimentary rocks. A review of these projects is given by, for instance, Tester et al. [35], Breede et al. [36], Xie et al. [37], Olasolo et al. [38], Lu [39] and Kumari and Ranjith [40]. This research increased knowledge of the behavior of rock masses and joints subjected to hydraulic stimulation. The current successful commercial projects (e.g., Soultz, Landau and Rittershoffen) were built upon this previously gained know-how. Furthermore, field-scale underground laboratories (e.g., Grimsel, EGS Collab and Utah FORGE) are tackling hydro-thermal-mechanical questions that have remained unresolved in the past. Moreover, although only a few sites are generating geothermal energy power (e.g., Soultz, Landau and Rittershoffen), all of the abandoned, suspended and ongoing projects are still important research facilities, providing a significant scientific database. For example, the Fenton Hill venture is described in great detail by Brown et al. [34] and a summary of the lessons learned can be read in Kelkar et al. [41]. A compilation of the development phases of the Rosemanowes geothermal project, problems faced and unresolved issues are provided by, for example, Parker [42]. Richards et al. [43] discuss the performance and characteristics of the Rosemanowes hydraulically stimulated geothermal reservoir. The authors also discuss the fundamental parameters controlling the impedance, thermal performance and water losses. The contribution of the Soultz project for the scientific community and its development phases are described in detail by, for instance, Genter et al. [44]. The lessons learned from past geothermal projects employing hydraulic stimulation treatments to crystalline rocks were taken into consideration in this study.

2.2. Nunavik’s Geothermal Potential

The geothermal energy potential of Nunavik and northern Québec has been investigated by Majorowicz and Minea [24] and Comeau et al. [45]. These regional-scale studies are based on scarce and sparse data distribution (Figure 1). In fact, only three deep boreholes with heat flow assessment exist in Nunavik (Camp Coulon, Raglan and Asbestos Hill mining sites), another one is in Nunavut (Nielsen Island) and one in Newfoundland and Labrador (Voisey Bay). All of these boreholes lie at distances of 430 to 500 km from Kuujjuaq. The heat flow estimated in these five sites ranges from 22 to 38 mW m−2 (Figure 1). Moreover, Comeau et al. [45] inferred the 1D subsurface temperature distribution among the different geologic provinces of Québec and their results suggest that at 5 km depth, the temperature in the Churchill Province ranges from 49 to 53 °C. However, only two data points exist for the Churchill Province, and these lie at a distance of almost 500 km from Kuujjuaq (Figure 1). Extrapolating such values to Kuujjuaq appears unrealistic, justifying a community-based approach to evaluate the local geothermal potential.

3. Geographical and Geological Setting

The settlement of Kuujjuaq, located north of the 55° parallel, is the administrative capital of the Kativik Regional Government and the largest northern village in the Nunavik region of Québec, Canada (Figure 2a; [46]). The 2016 census indicated 2754 inhabitants [46]. Diorite and paragneiss are the main lithologies, but smaller occurrences of tonalite, gabbro and granite are also observed (Figure 2b). A general description of these units is given in SIGÉOM (Système d’information géominière du Québec) [47].

3.1. Kuujjuaq’s Heating Demand

The community of Kuujjuaq experiences an annual average temperature of about 5.4 negative-degrees and an annual average of 8520 heating degree days below 18 °C (Figure 3; [15]). Although the residential dwellings are built to meet certain regulatory standards of insulation, the harsh climate results in high building heating requirements (Figure 3; [15]). The annual average fuel consumption of a typical residential dwelling in Kuujjuaq has been estimated as about 3100 to 8180 L [14,15]. This represents around 28 to 32 L m−2, with respect to the floor area. There are currently about 973 occupied residential dwellings in Kuujjuaq [33], indicating a total yearly consumption of 3 to 8 million L of oil for space heating. The peak heating load for a residential dwelling is 7 kW [14] to 15 kW, depending on the building heating load and floor area, indicating a peak load for the community of approximately 7–15 MW. The annual heating energy demand has been estimated between 21.6 and 71.3 MWh per dwelling, depending on the floor area [14,15]. Thus, the community’s heating demand is approximately 21 to 69 GWh per year.
The population in Kuujjuaq grew by 16% between 2011 and 2016 [41], representing an annual growth of about 3%. Translating this population increase into heating energy needs, and carrying out projections, suggests that in 30 years, the community’s annual heating demand may be 57 to 188 GWh. These values represent the threshold that the geothermal system designed in this study needs to meet during the 30 years of operation.

3.2. Previous Research Undertaken in Kuujjuaq

The lithologies outcropping in Kuujjuaq were sampled and a detailed description of texture, fabric, main mineral phases and major and trace geochemical elements is given in Miranda et al. [48]. Moreover, Miranda et al. [28,48] present the results of thermal conductivity, volumetric heat capacity, radiogenic heat production, density, porosity and primary rock permeability evaluated for the samples collected. The thermal conductivity and volumetric heat capacity were analyzed considering the samples at dry and water saturated states and for temperatures ranging between 20 and 160 °C [28].
Furthermore, a temperature profile was measured in a groundwater monitoring well located nearby the community and the terrestrial heat flux evaluated following a 1D inverse heat conduction approach, as explained in Miranda et al. [49] (Figure 4a). The heat flux assessment considered several hypotheses for the ground surface temperature history and variable conditions for the thermophysical properties. This approach suggested a heat flux at 10 km depth ranging between 32 and 69 mW m−2. The evaluated subsurface temperature distribution for a depth up to 10 km considered these results and is presented in Miranda et al. [28]. At 4 km depth, the subsurface temperature was estimated to range between 28 °C and 167 °C, with a median of 88 °C, depending on the paleoclimate history and thermophysical properties conditions. The minimum and maximum temperature values were defined in a deterministic manner and may correspond to the least probable scenarios.
Hypotheses for the specific discharge, Peclet number, hydraulic gradient and hydraulic conductivity are also discussed in Miranda et al. [49]. The specific discharge was calculated based on the hydraulic conductivity and hydraulic gradient. The permeability inferred for the paragneiss is 10−19 m2, the hydraulic conductivity was evaluated as 10−13 m s−1, the hydraulic gradient as 10−2, the specific discharge was inferred to be 10−14 m s−1 and the Peclet number 10−6.
Miranda et al. [50] carried out a preliminary evaluation of the fracture network in terms of geometrical and topological properties that was further improved in Miranda et al. [51]. Four main fracture sets were identified but only one, the N-S set, is optimally oriented for slip. The sets E-W and NNW-SSE require higher additional fluid overpressure to be reactivated. Additionally, Miranda et al. [51] present an a priori stress estimation for Kuujjuaq calibrated with published stress data for the Canadian Shield and built upon geological indicators, empirical correlations, analytical models and Monte Carlo simulations to deal with the ensemble of uncertainties (Figure 4b). The orientation of the maximum horizontal stress was inferred N210-220°E. All of these previous results constitute the basis for the numerical simulations carried out in the present study.

4. Methods

4.1. Numerical Simulator

The numerical simulator used in this study is an updated version of Jing et al.’s [52] FRACSIM3D for new joint constitutive laws.
The model requires as input the stimulated rock volume and a given stimulation injection pressure to calculate the 3D shape of the stimulated volume, amounts of shear slip, fracture apertures and the amount of fluid required. This simple approach assumes that the stimulation proceeds as a “shock wave” with low fluid pressure gradients within the already stimulated volume and steep gradients at the outer margins [53]. Following stimulations(s), the steady state flow through the modified fracture system is calculated followed by tracer test stimulations and heat extraction. Multiple well segments, vertical or inclined, can be specified for stimulation and for injection and recovery.
A shear displacement–dilation relationship, in which the shear dilation angle is a reducing function of displacement rather than constant, was formulated and implemented. The rate of opening with displacement associated with this new sliding/opening law is in line with many experimental publications and geological observations, for example, the work of Lee and Cho [54].
Slip commences when the ratio of the shear stress to the normal stress exceeds a threshold. This ratio is expressed as the tangent of the total friction angle, which is a property of both rock material and geometry of the fracture surface, and it can be derived from tilted block experiments. The total friction angle is a function of the normal stress, the total displacement and the accumulated fracture surface damage from past slip movement. This angle is the sum of the basic friction angle of a smooth rock surface and the shear dilation angle, which is the arctangent of the ratio of the amount of fracture opening per small increment of displacement. The shear dilation angle can be derived from the roughness of the natural fracture surface, and the dilation with slip may be derived for various representations of surface roughness and correlation from one side of the fracture to the other. The fracture opening responds to increased effective normal stress by closure and the rate of closure is known as the fracture normal compliance. Theoretical and data driven forms for the closure curve can be derived for specific models of joint roughness and material properties (e.g., [54,55]); for our current purposes, much of the required observational data is unavailable so a simple form is desirable.
The Jing et al. [52] version of FRACSIM3D aims to capture shear dilation and normal compliance in a single equation as:
w = w 0 + U t a n ( ϕ dilation ) 1 + 9 ( σ n σ n _ ref )
where w (m) is the aperture of the compliant fracture and w0 (m) is the initial aperture before induced slip, U (m) is the amount of slip (or shear displacement), σ′n (Pa) is the effective normal stress, σn_ref (Pa) is the reference stress for 90% closure (i.e., closure resistance of the fracture) and ϕdilation (°) is the shear dilation angle at low normal stress.
The shear dilation angle in the original Jing et al. [52] software version is assumed constant, meaning that a fracture tends to open at a constant rate as displacement increases at constant effective normal stress. Rock mechanics experiments (e.g., [54]) have suggested that the shear dilation angle increases over the first few millimeters of displacement as asperities of increasing wavelengths become out of phase, reaching a maximum and then decreasing to lower values after few millimeters of displacement since the longer wavelength asperities have less relative amplitude than those with the shorter wavelength. The assumption of a constant shear dilation angle was substituted by one that is a function of displacement to consider these observations. The relationship is as follows:
  • A linear increase from a starting shear dilation angle value to a maximum value over a certain small shear displacement distance;
  • An exponential decay with displacement thereafter to a low constant value at greater displacement at a user-specified rate.
Microseismic events recorded in the stimulation of crystalline rock may be caused by rapid changes in the shear dilation angle with displacement but rather caused by the sudden failure of macroscopic asperities or jogs. These asperities or jogs, not captured by the geometry of FRACSIM3D, might be related to fracture intersections, where upon the original propagation of one fracture is interrupted when it cuts another. The amount of shear stress needed to break these asperities is likely to be related to the asperity geometry, the number and extent of these asperities, the rock strength at the appropriate scale and the normal stress. To capture the possible effects of these macroasperities, an asperity strength factor was introduced into the code to give the fracture planes some extra resistance to slip. The asperity strength factor is converted to the rupturing shear stress required to overcome friction by:
τ = F asperity σ n t a n ( ϕ basic )
where τ (Pa) is the shear stress, Fasperity is the asperity strength factor and ϕbasic (°) is the basic friction angle.
A brief iterative solution provides the fluid pressure in the fracture required to cause asperity rupture via a reduction in asperity strength through a reduction in effective normal stress. The value of the asperity strength factor is further varied about the mean to create a population of different fracture asperity strengths. The introduction of this factor allows the generation of microseismic events of an appropriate magnitude and numbers within the FRACSIM3D model. Given sufficient stimulation injection pressure, an asperity strength factor of 0.4 and above will generate large simulated microseismic events, while a value of zero will suppress the magnitude of microseismic events.
To calculate post-stimulation flow, the fractures generated are embedded within a regular cubic lattice 3D discretization grid, in which the quantity of fluid flow from block to block is controlled by Darcy’s law with the transmissivity contribution from each fracture governed by the sum of products of the cubes of the fracture apertures and the fracture intersection length with the block face [52]:
Q j = i w h , i 3 l i 12 ω P ,   j = x , y , z
where Qj (m3 s−1) is the flow rate, wh (m) is the fracture flow aperture, li (m) is the fracture intersection length with the block face, ω (kg m−1 s−1) is the fluid dynamic viscosity and ∇P (Pa m−1) is the pressure difference between blocks. Since the aperture of each fracture intersection with a lattice cube face depends on the local fluid pressure, the steady state flow solution is necessarily iterative. A simple over-relaxation scheme was chosen for this computationally intensive part of the simulation to be able to easily interrupt the convergence to update apertures as the solution progresses, to keep the code free of commercial library subroutines and to fit the solution within the confines of the computer memory available at the time the code was written.
In other words, the flow pattern for the discrete fracture network is solved by converting it to an equivalent continuum mesh (with a finer length scale than typical fractures that affect flow) where the steady state flow through fractures is expressed using the “cubic law” (i.e., Reynold’s equation which is the solution of the Navier–Stokes equation for laminar flow between smooth parallel plates). Experiments (e.g., [54,55,56,57,58,59]) have suggested that the mechanical aperture of rough fractures overestimates the flow capacity due to the parallel plate assumption. Therefore, a flow aperture, somewhat less than the calculated mechanical aperture, is implemented. However, a user-definable input variable has been implemented to allow for manual input when desired.
Heat extraction is subsequently calculated under the approximate assumption that a thermal equilibrium is reached instantaneously between each solid element and the water passing through it. Moreover, heat transfer is constrained within the stimulated volume so that no heat transfer occurs at the boundaries of the total model volume. Heat is transferred between elements by both conduction and advection [60]:
ρ c T t = ( λ T ) ρ c fluid u T
where ρc (J m−3 K−1) is volumetric heat capacity, T (°C) is temperature, t (s) is time, λ (W m−1 K−1) is thermal conductivity and u (m s−1) is Darcy velocity.
The thermal energy extracted by the circulating fluid is given by:
e th = ρ c fluid Q recovered ( T recovered T injection )
where eth (W) is the thermal energy. The temperature in the reservoir varies over time as it cools. Further details on the joint constitutive laws, stimulation and steady state flow solutions, water loss approximation and heat extraction are given in Jing et al. [52].

4.2. Model Geometry

A cubic model volume of 4 km3 (with 1.6 km of edge length) discretized into a grid of 200 per 200 per 200 cells is used to carry out the numerical simulations. This grid was selected after carrying out a grid dependency study and the results revealed an influence of, on average, 5%. Slippage takes place in the fractures whose centers fall within the current (gradually refined) estimate of the stimulated volume configuration. Although this volume is fixed and defined by the user, the shape of the reservoir is adjusted progressively as the estimate of the stimulated permeability tensor is refined. The major axis of the stimulated area generally becomes oriented according to the inferred direction of the maximum horizontal principal stress. The boundaries of the model are assumed to be connected to a constant pressure boundary a few kilometers beyond the model boundary and are able to leak off fluid or draw it in during the steady state flow calculation. For a more detailed explanation of the treatment of far field fluid losses/gains, see the discussion in Jing et al. [52]. The model boundaries were kept at a constant temperature since the heat extraction volume happens away from the boundaries. At time zero, the fractures within the model volume are assumed to be filled with the in situ fluid and the pressure distribution is assumed to be hydrostatic. The initial in situ rock mass permeability, from which the initial fracture apertures are calculated, and the state of stress are defined by the user. The far-field, well beyond the model boundaries, is assumed to be at hydrostatic pressure during the stimulation and circulation calculations.
The engineered geothermal energy system was designed as a doublet, with one injector and one producer, with vertical wells. Two configurations for the wells were studied to identify the best location relative to the Lac Pingiajjulik fault (Figure 5) and maximum horizontal principal stress assumed direction (NE-SW). Different possibilities for well spacing and open hole length were considered. The well bore radius was defined as 0.11 m following the example of Soultz’s engineered geothermal energy system (e.g., [44]). The stimulation volume was considered variable. One stimulation was applied in each well and different stimulation and circulation pressures were studied. The operation time was defined for 30 years.
Three working hypotheses for the fracture network were generated and studied due to the current existing uncertainties, as explained in Miranda et al. [51]. These are:
  • Fracture intensity of 0.8 fractures m−1, fracture length and fractal dimension as sampled in the field (Figure 6a);
  • Fracture intensity of 0.8 fractures m−1, fracture length increased by a factor of 10 and fractal dimension as sampled in the field (Figure 6b);
  • Fracture intensity of 0.8 fractures m−1, fracture length as sampled in the field and fractal dimension decreased by a factor of 2 (Figure 6c).
The influence of the Lac Pingiajjulik fault was assessed by adding the fault plane to the model in a deterministic manner. The fracture was assumed to cut approximately through the middle of the model (Figure 5). This structure is described as a thrust fault with late dextral movement (e.g., [61]). Thus, its dip was assumed 45°NE and its dip azimuth N45°E but its exact orientation at depth is unknown. The radius of this fault was assumed to be longer than 5 km and the initial fault aperture at in situ effective stress was considered variable between 0.001 and 0.10 m. These apertures in no way reflect the actual geological fault offset, but rather were simply chosen to generate different flow apertures for the fault structure (low and high). The Lac Pingiajjulik fault, considering the estimated stress field [51], has high effective stresses acting across it at the present day.

4.3. Properties of the Medium

The engineered geothermal energy system was designed for a depth of 4 km since previous analyses suggest a 98% probability that an adequate thermal energy source to meet the community’s heating needs is available at this depth [28]. Preliminary simulations were carried out using a deterministic approach to define the worst-, base- and best-case scenarios in terms of properties of the medium (Table 1). It should be noted that the worst- and best-case scenarios may correspond to the least probable scenarios to occur since all of the best (and worst) values were grouped together.
The magnitude of the principal stresses and in situ fluid pressure were inferred by Miranda et al. [51] (Figure 4b) using empirical correlations, analytical modeling and Monte Carlo simulations calibrated with published literature data for the Canadian Shield. The subsurface temperature distribution was inferred in Miranda et al. [28]. Hypotheses for the in situ permeability are based on the results obtained at the Rosemanowes and Soultz stimulated geothermal projects [62,63]. The thermal conductivity is a function of temperature and pressure and was evaluated by Miranda et al. [28]. The volumetric heat capacity for water-saturated samples has been evaluated by Miranda et al. [28]. The Young’s modulus and Poisson’s ratio were defined based on the findings of Arjang and Herget [64] for the Canadian Shield. Asperity strength factor, initial, peak and ultimate shear dilation angle, peak and residual shear displacement and reference stress for 90% closure are working hypotheses to evaluate their influence on the performance of the hydraulically stimulated geothermal reservoir. It is important to highlight that these values are working hypotheses only and may be subject to wide variation. The basic friction angle was defined according to literature tilt experiments carried out in wet crystalline rock masses (e.g., [65]). The density of both in situ and circulation fluids was calculated as ([66] and references therein):
ρ fluid ( P , T ) = ρ fluid ( T ) + Δ ρ fluid ( P ) + 1000 × T D S
where TDS (kg L−1) is the total dissolved solids and P (Pa) stands for pressure.
The effect of temperature and pressure are given by ([66] and references therein):
ρ fluid ( T ) = 999.84 + 16.95 T + ( 7.99 × 10 3 ) T 2 + ( 4.62 × 10 5 ) T 3 + ( 1.06 × 10 7 ) T 4 + ( 2.81 × 10 10 ) T 5 1 + 0.017 T
Δ ρ fluid ( P ) = ( 4.0625 × 10 7 ) P
The in situ fluid was assumed brine with 100 g L−1 of total dissolved solids [67] and in thermal equilibrium with the geological medium. The dynamic viscosity of the in situ fluid was approximated as [66]:
ω ( T ) = ( 2.414 × 10 5 ) × 10 247.8 T 140
where T (K) is the absolute temperature.
Preliminary calculations were carried out to assess a baseline flow rate that the simulations need to meet. This was done as:
Q = e ρ c fluid ( T reservoir T injection )
where e (W) is the maximum heating energy demand.
The results suggest that flow rates in the range 14 to 45 L s−1 are required for the best-case scenario while for the base-case scenario, the flow rate range necessary is 27 to 89 L s−1. For the worst-case scenario, in turn, flow rates above 200 L s−1 are needed. The flow rates estimated for the best- and base-case scenarios have already been achieved in engineered geothermal energy systems (e.g., Soultz and Rittershoffen; [68]).

4.4. Levelized Cost of Energy

A first-order evaluation of the levelized cost of energy to build and operate an engineered geothermal energy system in Kuujjuaq at a depth of 4 km was carried out as (e.g., [69]):
LCOE = A total e annual = O annual + J annual e annual = O annual + i ( 1 + i ) t ( 1 + i ) t 1 I total e annual
where LCOE ($ MWh−1) is the levelized cost of energy, Oannual (¢ kWh −1) is the annual operation and maintenance cost, eannual (MWh) is the energy provided annually, i (%) is the imputed interest rate, t (year) is the project lifetime and Itotal ($) is the total capital investment. The latter is estimated as (e.g., [70]):
I total = C wells + C stimulation + C plant
where C ($) stands for cost.
Three different well cost models were used. These are (e.g., [71] and references therein):
WellCost   Lite = F s × 10 0.67 + 0.000334 ( z + h ) × N
ThermoGIS = F s × [ 0.2 ( z + l ) 2 + 700 ( z + l ) + 25000 ] × N × 10 6
HSD = 1500 × z × N × 10 6
where Fs is the well cost scaling factor, z (m) is depth, l (m) is the possible extra horizontal length of the well and N is the number of wells. HSD stands for hydrothermal spallation drilling. The parameter Fs was assumed 1 and l was assumed 0 in this work.
The imputed interest rate, stimulation, power plant and other surface facilities and operation and maintenance costs were assumed the same as the values proposed in Sanyal et al. [72]. However, these costs may be underestimated for Kuujjuaq. In fact, due to the remoteness, the infrastructure cost in Nunavik can be 2 to 5 times more expensive than the regular construction/infrastructure costs in southern areas [73]. Therefore, a first analysis was carried out using Sanyal et al. [72] costs and these were posteriorly increased by a factor ranging between 2 and 5 to assess more realistic energy costs for Kuujjuaq (Table 2).
The imputed interest rate was assumed 9% [72]. The annually provided energy was estimated based on the annual average heating (and electricity) consumption inferred for the community. The ensemble of uncertainties was simulated with the Monte Carlo method (Table 3) using the software @RISK [74] and the simulations were carried out for each possible engineered geothermal energy system design computed with FRACSIM3D.

5. Results

The numerical simulations carried out aimed at: (1) studying how the fracture network and properties of the medium influence the initial fracture aperture, fracture shear stiffness and fracture shear displacement, (2) analyze how the stimulation fluid pressure influences the number of sheared fractures and (3) designing an engineered geothermal energy system capable to meet the following targets:
  • Flow rates able to harvest enough thermal energy during all system operation time to meet the community’s heating demand;
  • Water loss maximum 20%;
  • Reservoir flow impedance lower than 1 MPa L−1 s−1;
  • Thermal drawdown lower than 1 °C/year.

5.1. Initial Fracture Aperture

The initial fracture aperture is estimated in FRACSIM3D through a predictive model that considers the measured or estimated average permeability of the undisturbed rock mass, the fracture radius, density, orientation and in situ stresses [52,75,76]. The results reveal that the hydromechanical properties significantly influence the estimated initial aperture of the fractures (Figure 7). The size of the fractures and fractal dimension also play a minor role (Figure 7). The observed variability consequently influences the hydraulic conductivity of the medium and the performance of the hydraulic stimulation procedure. These results were computed assuming the wells in configuration A (Figure 5).

5.2. Fracture Shear Stiffness

The fracture shear stiffness is approximated in FRACSIM3D as the ratio of the rock shear modulus divided by the fracture radius and multiplied by a geometric parameter [52,77]. Longer fractures, small Young’s modulus and high Poisson’s ratio are shown to considerably reduce the fracture shear stiffness of the medium (Figure 8), thus improving the performance of the hydraulic stimulation procedure. These results were computed assuming the wells in configuration A (Figure 5).

5.3. Fracture Shear Displacement

The shear displacement is expressed in FRACSIM3D as the excess shear stress over the required for slip divided by the fracture shear stiffness [52]. Thus, fractures with high shear stiffness tend to experience lower shear displacements and the latter tends to increase as a function of the applied fluid pressure (Figure 9). For the worst-case scenario, slippage occurs at stimulation, or fluid, pressures between 40 MPa and 70 MPa, regardless of the fracture network. For the base-case scenario, slip takes place at stimulation pressures of 20 to 60 MPa while, for the best-case scenario, between the pressures of 5 and 40 MPa. It is convenient to indicate that these results were computed assuming the wells in configuration A (Figure 5).

5.4. Sheared Fractures

The optimally oriented fractures to slip belong to the sets E-W and N-S, followed by the set NNW-SSE and the fractures within the other group (Figure 10). The set NW-SE is not optimally oriented to slip (Figure 10). These observations confirm the Mohr–Coulomb friction and slip tendency analyses carried out by Miranda et al. [51]. Moreover, increasing the stimulation fluid pressure led to an increase on the number of fractures that can be reactivated. However, it is important to highlight that the stress regime plays a major role and large uncertainty still exists as discussed by Miranda et al. [51]. As a note, these results were computed assuming the wells in configuration A (Figure 5) and a friction coefficient of 0.6.

5.5. Engineered Geothermal Energy System Design

The previous results suggest that developing a hydraulically stimulated geothermal reservoir is favored in a mechanically weak (i.e., low magnitude for the principal stresses and fractures with low resistance to deformation and opening) and hydraulically conductive medium. Moreover, they also indicate that smaller fractures tend to slip less, thus, making fluid circulation more difficult through these fractures. Reservoirs of potential interest were not simulated for the worst-case scenario, primarily due to the low reservoir temperature and high flow rates required to meet the demand, and secondly due to the high magnitude inferred for the principal stresses and high resistance of the fractures to slip and opening. Thus, the following sections will be focused on the best- and base-case scenarios, discussing the best simulated design to meet the water loss, reservoir impedance and thermal drawdown targets. The influence of the Lac Pingiajjulik fault and its initial aperture and, thus, hydraulic conductivity, is also considered.

5.5.1. Best-Case Scenario

Reservoirs capable of meeting the aforementioned targets were simulated in this scenario for the three fracture network hypotheses regardless of Lac Pingiajjulik initial fault aperture (Table 4). As previously outlined, these simulations were carried out bearing in mind the estimated required flow rate of 14 to 45 L s−1. The simulations reveal that, in this case study, positioning the wells perpendicular to the fault trace (configuration A; Figure 5) leads to a better performance and sustainability of the system than placing the wells parallel to the fault trace (configuration B; Figure 5). The results also suggest that longer fractures (fracture network 2) tend to improve the performance of the engineered geothermal energy system. Lower stimulation and circulation pressures are required to enhance the transmissivity and circulate the fluid.

5.5.2. Base-Case Scenario

As previously outlined, these simulations were carried out bearing in mind the estimated required flow rate of 27 to 89 L s−1. Reservoirs of potential interest were only generated for fractures longer than sampled in the field (fracture network 2). However, to meet the defined targets, higher stimulation volume and stimulation and circulation pressures were required than for the best scenario previously studied (Table 5).

5.6. Thermal Energy and Potential Heat Output

The thermal energy extracted by each design previously studied is capable of fulfilling the forecast heating energy demand over 30 years of operation (Figure 11). Moreover, the results indicate that configuration A (Figure 5) is more appropriate to meet the upper bound of the projected demand while configuration B (Figure 5) would more likely meet the lower bound of the demand.
Some of the designs largely exceed the heating needs within the first years of the geothermal system operation. Therefore, an additional study was undertaken to assess if the geothermal system can be used for combined heat and power to improve the sustainability of the system. The electricity consumption in Kuujjuaq was evaluated as 18.9 GWh (see [28] for further details) and the yearly growth in energy demand is assumed to be 3% [46]. The results suggest that about 6 to 8 MWth are rejected per each MWe produced when considering the deep geothermal energy sources in Kuujjuaq. If only 50% of the waste heat from power production can be used for space heating, since the remaining is lost from parasitic equipment, an engineered geothermal energy system can be designed to operate in combined heat and power mode during approximately 8–15 years of the project lifetime, meeting both power and heating demand (Figure 12).

5.7. Recovery Factor

The thermally active reservoir volume changes over time as heat is being extracted. The recovery factor was calculated as the ratio between the thermally active volume of rock during the 30 years of operation and the defined stimulated volume (Figure 13). The recovery factor is commonly assumed to be within the theoretical range of 2–20% (e.g., [35,78]). The results obtained agree with these values.

5.8. Levelized Cost of Energy

The three well cost models used in this work suggest that the price of two 4-km wells can be USD 9.8 M (WellCost Lite model; Equation (12a)) and USD 12 M (ThermoGIS and HSD models; Equation (12b) and Equation (12c), respectively). Considering the Sanyal et al. [72] proposed costs, the stimulation of both of the wells ranges between a minimum of USD 1 M and a maximum of USD 2 M, with a mean value of USD 1.5 M. The power plant and other surface facilities estimated for each possible design previously studied (Table 4 and Table 5) is found to vary between an average maximum cost, for configuration A (Figure 5), of USD 54 M and an average minimum of USD 34 M (Figure 14a). For configuration B (Figure 5), the cost varies between an average minimum of USD 12 M and an average maximum of USD 28 M (Figure 14b). Thus, the average minimum capital cost estimated for configuration A (Figure 5) is USD 46 M while the average maximum is USD 66.8 M (Figure 15a). For configuration B (Figure 5), the average minimum is USD 24.8 M and the average maximum is USD 40.8 M (Figure 15b). All these costs are in USD.
However, these costs may be too optimistic since, due to the remoteness, the infrastructure cost in Nunavik can be 2 to 5 times more expensive than the regular construction/infrastructure costs in southern areas. Therefore, a second analysis was carried out increasing the costs by a factor of 2 and 5 (Table 2). The stimulation of both wells is found to vary between USD 2 M and USD 4 M, for a factor of 2, and between USD 5 M and USD 10 M if the costs are increased by 5 times. The power plant and other surface facilities, for a factor of 2, can have a cost varying within the average minimum of USD 24 M to the average maximum of USD 108 M (Table 6) or between USD 60 M and USD 270 M (Table 6) if the cost is 5 times higher than the one proposed by Sanyal et al. [72].
Consequently, the capital cost for each simulated design increases to an average minimum of USD 38 M and an average maximum of USD 122 M if the costs are 2 times higher and to USD 79 M and USD 289 M if the costs are 5 times higher than in Sanyal et al. [72] (Table 7).
Sensitivity analyses indicate that, for configuration A (Figure 5), the power plant and surface facilities cost is the most influential parameter on the capital cost followed by the well cost (Figure 16). The Spearman correlation coefficients range from 0.78 to 0.91 and 0.28 to 0.54, respectively, and 0.02 and 0.13 for the stimulation costs. For configuration B (Figure 5), the power plant and surface facilities and well costs have a similar influence on the capital cost (Figure 16). The Spearman correlation coefficients vary between 0.54 and 0.70 and 0.59 and 0.73, respectively, and between 0.11 and 0.20 for the stimulation cost. Although the values displayed in Figure 15 correspond to the costs proposed by Sanyal et al. [72], the same trend is observed if the costs are 2 or 5 times higher.
Considering the results previously discussed, the levelized cost of energy was evaluated considering the Sanyal et al. [72] costs as the optimistic scenario, the increase by a factor of 2 as the most likely scenario and the increase by a factor of 5 as the pessimistic scenario. The mean value for each of these hypotheses was used in this analysis. The levelized cost of energy considering configuration A (Figure 5) and the best-case scenario (Table 4) for heat production only was evaluated to range between an average minimum of 121 USD MWh−1 and an average maximum of 626 USD MWh−1 (Table 8). These values decreased to 120 USD MWh−1 and 617 USD MWh−1 if the geothermal system is assumed to work in combined heat and power mode during 10 to 15 years (Table 8). The remaining designs can only produce heat. The levelized cost of energy is found to range within an average minimum of 54 USD MWh−1 and an average maximum of 475 USD MWh−1 (Table 8).
Currently, the heating energy cost in Kuujjuaq has been evaluated as 0.19 USD kWh−1 (or 190 USD MWh−1) and electricity as 0.6 USD kWh−1 (or 600 USD MWh−1; [11,79]). Thus, a probabilistic study was undertaken to assess if the deep geothermal energy source, for both heat production and electricity generation, is cost-competitive compared to current oil furnaces for space heating and diesel power plants for electricity. The results for heat production reveal that configuration A for the best-case scenario have a probability of 8 to 18% of providing heating energy at a lower energy cost than the business-as-usual scenario (Figure 17a). A probability of 25 to 33% was estimated for configuration A for the base-case scenario (Figure 17b), while configuration B reveals probabilities of 49 to 91% (Figure 17c). Per contra, the probability of combined heat and power to provide electricity at a lower cost than the business-as-usual scenario is higher than 99% (Figure 17d). This analysis was undertaken by combining the results obtained for the optimistic, most likely and pessimistic scenarios. Considering the most likely scenario only, the engineered geothermal energy system has 13 to 100% probability of providing heating energy at a lower cost than the oil furnaces and 100% probability of providing electricity at a lower cost than the diesel power plants.
A detailed analysis was undertaken to assess how the probability of providing affordable heating energy with an engineered geothermal energy system can be increased. The results revealed that the energy consumption and capital cost play the major roles, with Spearman coefficients of, on average, −0.70 and 0.60, respectively. The operation and maintenance costs, per contra, reveal a weak correlation (on average 0.02) with the levelized cost of energy. This trend is observed for all of the simulated scenarios (Figure 18). In fact, for the best- and base-case scenarios assuming configuration A, if the energy consumption is near or greater than its maximum value and the capital cost near or below its minimum value, then the geothermal systems simulated will have a greater probability of providing heating energy at a lower cost than the business-as-usual scenario (Figure 19). However, for configuration B at the best-case scenario, the energy consumption only needs to be greater than its 30th to 50th percentile and the capital cost below, its 50th percentile (Figure 20). For the base-case scenario, in turn, the energy consumption can be near its minimum value and the capital cost near its maximum value (Figure 20).

6. Discussion

The current energetic framework of the 239 Canadian off-grid communities relying solely on diesel entails: (1) high costs to buy, transport and store diesel, (2) carbon emissions that contribute to climate change, (3) potential spills and leakages that damage the local environment and (4) low energy security that constrains community development (e.g., [80]); thus, opposite to the three main axes that constitute the concept of sustainable energy markets: economic affordability, environmental compatibility and energy security [69,81]. A sustainable energy market intends to maintain and improve living standards at an affordable cost by replacing the consumption of environmentally harmful sources of energy by environmentally friendly alternatives. Therefore, introducing renewable off-grid technologies into the communities can break this energy poverty cycle that has held back their socio-economic progress. Deep geothermal energy can be a viable alternative to fossil fuels in remote northern communities, if the following criteria are met [82]:
  • It is sufficiently abundant to meet a significant percentage of the market demand;
  • It can be obtained at a cost competitive with existing energy sources.
A previous deep geothermal energy source assessment suggested that the thermal energy stored beneath the community of Kuujjuaq at depths greater than 4 km was capable of fulfilling an estimated community heating energy demand [28]. Thus, the current study was carried out with the goal of designing an engineered geothermal energy system to harvest this energy and evaluating its economic potential compared to the current use of fossil fuels. However, the results from this study are subjected to high levels of uncertainty due to the current poor knowledge of both the geological structures and thermo-hydro-mechanical properties of rock at depth in the area. Nevertheless, this work is a contribution to assess if deep geothermal energy sources can supplant or displace reliance on fossil fuels. Moreover, although this work was undertaken in the community of Kuujjuaq, the methodological approach followed can be extended to other remote northern settlements facing important energy issues and with similar geothermal exploration data gaps.
The first research questions of this study were: Will the hydraulic stimulation technique applied in crystalline basement rocks develop a well-connected flowing system in Kuujjuaq? How can this be done? What further local geological and thermo-hydro-mechanical data is required for more accurate predictions? A shear-dilation-based model (FRACSIM3D; [52]), updated for new joint constitutive laws, was used in this work to help on the design of an engineered geothermal energy system for each uncertain geological (fracture network) and thermo-hydro-mechanical (properties of the medium) scenario. Although McClure and Horne [83] propose different hydraulic stimulation mechanisms (pure opening mode, pure shear stimulation, primary fracturing with shear stimulation leakoff and mixed-mechanism stimulation), shear displacement is still the most widely accepted mechanism to explain permanently enhanced transmissivity of natural fracture networks (e.g., [44,84,85]). Jacking (i.e., fracture-normal dilation) can also occur near the open-hole section, but its effect declines over time with the depressurization (e.g., [84,85]), but may increase over time with thermoelastic effects following cooling. Nevertheless, field observations at the Fenton Hill test site have supported mixed-mechanism stimulation, involving the propagation of hydraulic splay fractures due to stress changes induced as fractures opened and failed to shear [86]. Thus, the present state of knowledge highlights the importance of continued field tests to better characterize the subsurface and understand stimulation mechanisms.
Studies have shown that considering the effect of thermoelasticity tends to improve the flow rates, decrease water loss and reduce the hydraulic impedance (e.g., [87]). In this study, although a thermoelasticity module has already been developed for an earlier version of FRACSIM3D [87], the authors focused only on the hydro-mechanical and thermo-hydro coupling processes. Nevertheless, future improvements could include this module and assess the effects of considering only poroelasticity and considering both thermo and poroelasticity effects. The simulations reveal that the development of hydraulically stimulated geothermal reservoirs is favored in a high temperature, hydraulically conductive and mechanically weak (in terms of low magnitude of the principal stresses and low resistance of fractures to deformation and opening) subsurface, i.e., the best-case scenario assumed in this study. If the reservoir temperature and subsurface hydraulic conductivity are decreased and the magnitude of the principal stresses and the fractures resistance to deformation and opening is increased, then developing engineered geothermal energy systems is more difficult. In fact, the base-case scenario assumed in this study only generated reservoirs of potential interest assuming fracture longer than those sampled in the field. Moreover, if the worst-case scenario is the prevailing scenario in the study area, then engineered geothermal energy systems are not a viable alternative off-grid technology to offset fossil fuels’ consumption in the community of Kuujjuaq. Nevertheless, the work conducted suggests additional geothermal exploration to improve these predictions is warranted. Exploratory boreholes deeper than 300 m are necessary to carry out more accurate estimations of the subsurface temperature, stress field regime and fracture network. Hydraulic stimulation field tests are also required to obtain more accurate information on the hydro-mechanical behavior of the subsurface. Another important aspect to accurately assess is the characterization of the fault hydraulic properties that may play a major role for the development of engineered geothermal energy systems in Kuujjuaq. Although numerical models are fundamental to support the design and operational planning of engineered geothermal energy systems, reservoir simulation is a complex endeavor that requires careful calibration [88]. Hydraulic stimulation is accompanied by induced seismicity and the cloud of events is an important monitoring tool to understand the direction in which the reservoir can develop (e.g., [34,63,89]). Moreover, fracture seismic observations recorded before, during and after hydraulic stimulation enable the building of time-lapse images mapping the fluid flow in the rocks, and thus estimate the reservoir connectivity changes over time (e.g., [90]). The seismic method used to analyze fracture behavior during stimulation is based on the recording and analyzing of passive seismic data. Drilling the wells without previous knowledge of the reservoir growth direction, jointing pattern and in situ stress field can severely compromise a stimulated geothermal project (e.g., [34,89]). A lesson learned from Rosemanowes’s stimulated geothermal project is that the wells should be drilled with their azimuths parallel to the minimum in situ horizontal stress direction (e.g., [43]). In the current work, the wells were assumed to be vertical, but further simulations can be envisioned by changing their dip, dip directions and relative azimuths. However, the uncertainty associated with the direction of the stress field must be resolved at an early stage. Microseismic monitoring can also be used to infer the fracture shear displacement (e.g., [63]). Additionally, an artificial neural network–genetic algorithm-based displacement back analysis has been proposed by Zhang et al. [91] to estimate fracture stiffness, in situ stresses and elastic parameters from borehole displacements during drilling. Laboratory tests are also an important component to enhance the understanding of joint behavior (e.g., [92,93,94,95,96,97,98,99]). However, care must be taken due to scaling effects. Nevertheless, the laboratory tests would provide first-order calibrations for the fractures opening parameters. Laboratory tests were not in the scope of this study but may be beneficial to refine the range of applicable input parameters. Moreover, further simulations can be carried out by changing the configuration from a doublet to a triplet, for instance. Two injectors and one producer may help to develop hydraulically stimulated geothermal reservoirs for the base-case scenario with shorter fractures. This can, however, impact the capital cost attributable to wells.
Chemical mineral dissolution by alkaline or acidic additives was used at Soultz to clean joints, faults and pore volume and, thus, helping to enhance the natural permeability of these structures [44,100]. The application of chemical stimulation may be an additional option for further studies in Kuujjuaq if future geothermal exploration reveals the presence of fractures infilling material that can react to this procedure. If this is observed, then a 3D water/rock chemical interaction module has been developed for FRACSIM3D [101] that can be used to simulate both hydraulic and chemical stimulations and evaluate the system performance response. Moreover, laboratory experiments with lithologies analogous to the potential reservoir rock in Kuujjuaq can provide further insights as to whether chemical stimulation may be effective (e.g., [102]). Thermal stimulation can also be applied to increase the near-wellbore productivity (e.g., [103,104]).
The numerical simulations carried out considered water as the working fluid. Further improvements can be carried out to assess the change in performance if supercritical CO2 is used instead. This concept was proposed by Brown [105] and the advantages of operating engineered geothermal energy systems with supercritical CO2 are twofold: (1) CO2 has certain thermophysical and chemical properties that make it an attractive heat transfer medium and (2) such systems can promote geological storage of CO2 as an ancillary benefit. Numerical simulations have been carried out by several authors, and their results suggest a significant performance improvement compared to water [106,107,108]. Pruess’s [106] numerical simulations reveal thermal extraction rates approximately 50% larger for CO2 compared to water. This suggests that flow rates 50% lower than the ones currently estimated for Kuujjuaq may be enough to extract the same amount of energy. In 2012, the 14 diesel power plants in Nunavik contributed with 65,000 tonnes/year of CO2 equivalent [109]. This equals to a flow rate of approximately 2.1 L s−1. Assuming that the required flow rate to meet Kuujjuaq’s minimum heating demand threshold is now 7 L s−1 instead of the 14 L s−1 estimated using water as the working fluid, then 3.5 years of capture and storage of CO2 are needed to reach that flow. This suggests that although CO2 can be captured from the thermal plants, used in engineered geothermal energy systems (e.g., [110]) and its thermal performance was proven superior to water (e.g., [106]), this solution is not technically or economically viable for northern communities. For example, the amount of CO2 produced by a single diesel power plant is insufficient to meet Kuujjuaq’s demand and shipping additional CO2 from other communities and southern areas can increase the costs of the geothermal project and may bring additional environmental impacts. A review of large-scale CO2 shipping has been undertaken by Baroudi et al. [111] and the transport costs range between USD 10 and USD 167 per tonne CO2, depending on the travel distance and transport capacity. An alternative working fluid, supercritical N2O, was proposed by Olasolo et al. [112] due to its thermodynamics properties but may still not be a viable solution for remote northern communities.
The second research question was: Are the deep geothermal energy sources harvested by engineered geothermal energy systems in Kuujjuaq cost-competitive compared to fossil fuels? The levelized cost of energy was estimated based on literature values for the US proposed by Sanyal et al. [72] and increased by factors of 2 and 5 to be more in the range expected in remote northern regions. Empirical well cost models were also used (e.g., [71] and references therein). Two 4-km-deep wells were inferred to cost from USD 9.8 to 12.1 M. However, Minnick et al. [25] refer that in Nunavut (northern Canada) a full-size 4-km production well can cost approximately USD 12 M. This corresponds to twice the estimated price in this work since Minnick et al.’s [25] estimate includes expenses associated to Nunavut’s challenging environment. Such a difference leads to an increase of 14 to 31% of the capital cost, which in turn increases the levelized cost of energy. The levelized cost of energy, considering the literature well cost models and assuming heat production only, ranges between 83 USD MWh−1 and 265 USD MWh−1 for the likely scenario. Doubling the well costs to be in line with anticipated Nunavik costs leads to values varying between 108 USD MWh−1 and 345 USD MWh−1.
Nevertheless, the first-order evaluation carried out in this study suggests that engineered geothermal energy systems may have commercial interest. The probability of geothermal heat production being less expensive than the business-as-usual scenario (oil furnaces) ranges from a minimum of 8% to a maximum of 91%. This probability can be increased if the energy consumption of the community is near or above its maximum estimated value and if the capital cost is decreased to its minimum value.
Furthermore, if the best-case scenario prevails in the study area, then combined heat and power may be a viable option to not only increase the sustainability of the geothermal system, but also to decrease the levelized cost of energy and, thus, increase the economic potential. In fact, combined heat and power have more than 99% of probability to provide electricity at a lower cost than the current diesel power plants. It is important to highlight that the geothermal system was assumed to be working 24 h per 7 days during the 30 years of project lifetime. The most likely levelized cost of energy, for heat production only, in this study, was found to range between 83 USD MWh−1 and 265 USD MWh−1. For combined heat and power, the most likely levelized cost of energy ranges between 218 USD MWh−1 and 262 USD MWh−1. These values are within the range referred by Tester et al. [35] and Augustine [113], 100 to 1000 $USD MWh−1 and 140 to 310 USD MWh−1, respectively. However, Tester et al. [35] also mentioned that with mature and cheaper technology, the levelized cost of energy can reach values as low as 36 to 92 USD MWh−1. Note that the levelized cost of energy described by Tester et al. [35] and Augustine [113] are for electricity generation only and not heat production or combined heat and power as studied in Kuujjuaq.
A previous economic potential assessment of geothermal energy sources in northern Canada undertaken by Majorowicz and Grasby [23] suggests a total cost of the project of USD 26 M. These estimations were carried out for a depth of 3 to 5 km, doublet spaced 550–700 m, reservoir temperature of 120 °C and flow rate of 30 L s−1. The remaining estimations presented by these authors were done for electricity generation only, thus difficult to compare with the results obtained in the current study. Nevertheless, these authors evaluated a cost ranging from 0.50 to 0.84 USD kWh−1. Additionally, Richard [114] carried out an assessment of the viability of engineered geothermal energy systems for electricity generation in Québec. The base-case estimated production costs were found to vary between 1.25 USD kWh−1 at 4 km and for a reservoir temperature of 100 °C and 0.18 USD kWh−1 at 10 km depth and 250 °C using a reduced drilling cost. In this work, an engineered geothermal energy system in Kuujjuaq for direct use applications indicates a most likely cost of 0.08 to 0.27 USD kWh−1, while combined heat and power is 0.22 to 0.26 USD kWh−1. Thus, deep geothermal energy sources harvested by engineered geothermal energy systems may be cost-competitive compared to fossil fuels, and further gathering of information is worthwhile to improve the quality of these predictions.

6.1. Comparison with Other Deep Geothermal Energy Projects

Engineered geothermal energy systems, or enhanced geothermal systems (EGS), are always site-specific in terms of both the subsurface environment and the local energy demand and price sensitivity. Thus, the comparison of the results obtained in this study with other deep geothermal energy projects developed in different geologic–economic contexts is challenging. The thermo-hydro-mechanical characteristics of the subsurface, which highly influence the design and development of EGS, vary from place to place. To date, the technological and economic feasibility of EGS in the Canadian Shield has not yet been fully assessed. Due to a lack of field experiments, history matching to calibrate the numerical simulations and restricting the range of assumptions made in this study is not possible. Nevertheless, in order to design an EGS for Kuujjuaq, the approach was strongly influenced by the results obtained and the lessons learned from the hydraulically stimulated geothermal projects of Fenton Hill, Rosemanowes, Soultz and Rittershoffen [34,41,42,43,44], many of which are incorporated into the design of the FRACSIM3D model. For example, the acceptable targets chosen for water loss, thermal impedance and thermal drawdown were based on targets adopted for those geothermal projects (e.g., [34,41,42,43,44]) even though the economics at Kuujjuaq differ substantially. The flow rates applied to the numerical models also considered the values targeted and observed in the commercial projects of Soultz and Rittershoffen [68]. The effective stimulation pressure required to induce shear slip in the pre-existing natural fractures is also found to be comparable to the values found at those four geothermal sites [37]. The potential heat and power output from an EGS developed in Kuujjuaq is also similar to the power generated and heat produced at Soultz and Rittershoffen [68]. The Soultz-sous-Forêts geothermal power plant is producing electricity using an organic Rankine cycle and has an installed gross capacity of 1.7 MWe. The results of this study suggest a potential electricity production of 2.9 to 3.5 MWe provided that the “best case” in situ temperatures are found. The Rittershoffen geothermal plant is providing superheated water for industrial needs and has an installed capacity of 24 MWth. The results in Kuujjuaq suggest a potential heat production between 5.7 and 29 MWth. Finally, the capital cost and levelized cost of energy are also within values published for northern Canada and southern Québec [23,114].
A geothermal cascade system can be an option for the deep geothermal energy source in Kuujjuaq as it has been proposed in southern Poland [115]. In southern Poland, a deep geothermal well with a geothermal water temperature of 82 °C and a maximum flow rate of 51.22 kg s−1 has been considered to assess the potential for electricity generation. The thermodynamic calculations carried out by Kaczmarczyk et al. [115] indicate that the gross capacity in the most optimistic variant will not exceed 250 kW for the organic Rankine cycle and 440 kW for the Kalina cycle, and that the gross electricity generation will not exceed 1.9 GWh/year for the organic Rankine cycle and 3.5 GWh/year for the Kalina cycle. If the same order of values is found in Kuujjuaq, electricity generation from geothermal waters of 82 °C and flow rates of 51 kg s−1 does not seem to be a viable option to replace diesel since the annual average electricity demand in the community is almost 19 GWh. Nevertheless, the use of geothermal energy, regarding the applications, has several positive aspects, both environmental and social, that can significantly improve the living conditions at Kuujjuaq and, indeed, in any region as highlighted by, for instance, Operacz and Chowaniec [116] and Sowizdzal et al. [117].

6.2. Comparison with Other Renewable Energy Alternatives

Another geothermal solution to reduce diesel consumption in remote northern regions is ground-coupled heat pumps [15] and borehole thermal energy storage [118]. Within the components of a shallow geothermal system exploited by ground-source heat pumps, the ground heat exchanger is a key component, playing a major role in achieving a high coefficient of performance [119]. Important parameters for the ground heat exchangers include the geometry, the pipe material, the backfill material, the working fluid, the depth of the ground heat exchanger, the heat transfer coefficient, the outlet temperature, the thermal resistance and the pressure drop [119,120,121]. Although ground-source heat pumps can be viable to supply heat to residential dwellings [15], the extracted energy with such shallow systems (100 to 300 m) is generally insufficient to fulfill the heating demand of a single residential dwelling, relying on an energy input to cover the remaining load.
Biomass, wind, solar, biofuels and hydro have also been studied as alternatives to offset diesel consumption in remote Canadian communities [6,7,10,12,14,17]. Thompson and Duggirala [7] and Yan et al. [14] indicate biomass as the most favorable and competitive renewable energy technology compared to natural gas, gasification of domestic waste, wind and solar. Furthermore, Stephen et al. [10] carried out a study to determine the techno-economic feasibility of biomass utilization for space heating and concluded that biomass has the potential to reduce heat costs, reduce the cost of electricity subsidizations for electrical utilities, reduce greenhouse gas emissions and increase energy independence. However, biomass resources need to be transported and stored similarly to diesel, which is disadvantageous compared with the development of a local source of energy. A fully hybrid wind-solar-battery-diesel system was selected by a multi-objective genetic algorithm as a viable solution for northern communities [17]. This system suggests a reduction of 50% of the levelized cost of energy compared to a diesel-only scenario and it may help to displace 675 MWh/year of energy from diesel. However, wind and solar remain weather-dependent and, therefore, their supply of energy is intermittent [122]. McFarlan [12] conducted a techno-economic assessment of replacing diesel for electricity generation with clean biofuels (methanol and dimethyl ether). The results revealed an increase of the cost compared to the diesel scenario. However, McFarlan [12] argued that although clean biofuels are more expensive than diesel there is potential socio-economic benefits from switching to these energy sources. Biofuels are still in a nascent stage and more research and development are needed. Micro-hydropower systems in off-grid communities have been identified as a favorable solution to displace the use of diesel and reduce greenhouse gas emissions [6]. Most of these systems are run-of-river without a need for a dam or reservoir. The best geographical locations would be steep rivers, streams, creeks or springs flowing year-round [6]. Although water resources may not be a current issue throughout Canada, climate change may impact existing and proposed hydropower projects, especially in the boreal, subarctic, arctic unique and complex environments [123].

7. Conclusions

Geothermal energy off-grid technologies are a potential solution for improving the energetic framework of the 239 Canadian remote northern communities that rely solely on diesel for electricity and space heating. Although subject to high uncertainty, the results of this work suggest that engineered geothermal energy systems are technically and economically viable in the community of Kuujjuaq (Nunavik, Canada) to supplant the reliance on fossil fuels. A “what-if” approach was followed in this study to deal with the poor subsurface knowledge. The engineered geothermal energy systems designed for each highly uncertain scenario needed to provide enough thermal energy for direct use applications during 30 years of operation and to be within certain defined performance parameters limits. These were: water loss lower than 20%, reservoir flow impedance lower than 1 MPa L−1 s−1 and thermal drawdown lower than 1 °C/year. The numerical simulations revealed that developing hydraulically stimulated geothermal reservoirs is favored in a high temperature, hydraulically conductive and mechanically weak subsurface (in terms of low magnitude of the principal stresses, high differential stress and low resistance of fractures to slip and opening). Decreasing the reservoir temperature and hydraulic conductivity and increasing the magnitude of the principal stresses and the fractures resistance to slip and opening decreases the performance of the system. Moreover, where an NW-SE striking fault is present, placing the wells relative positions parallel to the inferred direction of the maximum principal stress (NE-SW) revealed better performance than if the wells are located parallel to the fault plane. This may also be caused by an overlap of the stimulation volumes. In fact, the simulation results revealed to be relatively insensitive to the assumed initial fault apertures, suggesting that the fault itself plays little role in fluid flow circulation. However, if the fault segment runs quasi-parallel to the maximum principal stress (and so being capable of further stimulations) may make a significant difference to the results. Furthermore, longer fractures tend to improve the performance of the system. Smaller fractures have higher shear stiffness and tend to slip less, making fluid circulation more difficult. Additionally, the best-case scenario suggests that combined heat and power is possible during the first 10 to 15 years of the geothermal system operation if the in situ temperature assumptions are met.
A first-order evaluation of the capital cost and levelized cost of energy was carried out using the Monte Carlo method. The results for capital cost is in the range USD 25 to 67 M$ for the optimistic scenario, USD 38 to 122 M for the likely scenario and USD 79 to 289 M for the pessimistic scenario. The global sensitivity analysis based on the correlation coefficient revealed that the power plant and surface facilities cost is the most influential parameter on the capital cost, followed by the well cost and stimulation cost. The levelized cost of energy, assuming heat production only, was estimated to range within 54 and 145 USD MWh−1 for the optimistic scenario, between 83 and 265 USD MWh−1 for the likely scenario and between 170 to 626 USD MWh−1 for the pessimistic scenario. Combined heat and power revealed a levelized cost of energy varying between 120 and 143 USD MWh−1 for the optimistic scenario, 218 and 262 USD MWh−1 for the likely scenario and 510 and 617 USD MWh−1 for the pessimistic scenario. The probabilistic analysis carried out indicates that engineered geothermal energy systems have an 8 to 91% probability of providing heating energy at a lower cost than the current oil furnaces and more than a 99% chance of providing electricity at a lower cost than the diesel power plants currently in place. Given that geothermal energy is a local source available for remote community, further geothermal exploration is recommended and indispensable to decrease the existing uncertainties and support decisions to develop this energy alternative. Hence, helping remote northern communities moving towards a greener and sustainable energetic future.
The work described in this study highlights several uncertain geological and thermo-hydro-mechanical parameters that need further gathering of information to obtain more accurate estimates of the techno-economic potential of EGS in Kuujjuaq. Exploratory boreholes deeper than 1 km are required for accurate heat flux and subsurface temperature assessments (e.g., [34,42,43,44,124,125]). These boreholes are also needed to carry out hydraulic tests (e.g., [34,42,43,44,124]) and stress measurements (e.g., [34,42,43,44,124]). Borehole televiewer can also be useful to have a better characterization of the underground fracture network (e.g., [34,42,43,44,124]). The effect of temperature (or thermal lift) on the hydraulic properties (e.g., [126]) could be additionally evaluated. Thus, the next step that is justified by this research is the drilling of exploration boreholes to obtain more accurate data for the numerical simulations that can be improved with history matching.

Author Contributions

Conceptualization, M.M., J.R., J.W.-R. and C.D.; methodology, M.M. and J.W.-R.; software, J.W.-R.; validation, M.M.; formal analysis, M.M.; investigation, M.M.; resources, M.M., J.R., J.W.-R. and C.D.; data curation, M.M.; writing—original draft preparation, M.M., J.R., J.W.-R. and C.D.; writing—review and editing, M.M., J.R., J.W.-R. and C.D.; visualization, M.M.; supervision, J.R. and C.D.; project administration, J.R.; funding acquisition, J.R. All authors have read and agreed to the published version of the manuscript.

Funding

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.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All the relevant dataset are presented in this study.

Acknowledgments

The authors would like to acknowledge Félix-Antoine Comeau, Inès Kanzari, Jean-François Dutil, Sérgio Seco and Stefan Premont for the support during the analyses of the thermophysical properties. Acknowledgments are extended to Cynthia Brind’Amour-Cote for the support during the field campaign. A special thanks is given to Fiona Chapman for the English language review. The authors are also grateful to the three anonymous reviewers whose comments and suggestions helped to make the manuscript clearer.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Notation

SymbolDefinitionUnit
AAream2
CCost$
eThermal energyW; MWh
FFactor---
HHydraulic impedanceMPa L−1 s−1
iImputed interest rate%
ITotal capital investment$
I*Injection well---
lLengthm
LCOELevelized cost of energy$ MWh−1
NNumber---
OAnnual operation and maintenance cost¢ kWh −1
PPressurePa
PPressure gradientPa m−1
QFlow ratem3 s−1
qHeat fluxmW m−2
R*Recovery well---
TTemperature°C; K
tTimes; year
TDSTotal dissolved solidskg L−1
UShear displacementm
uDarcy velocitym s−1
VVolumem3; km3
WWater%
wAperturem
zDepthm
Greek letters
λThermal conductivityW m−1 K−1
ρcVolumetric heat capacityJ m−3 K−1
σPrincipal stressPa
σ’nEffective normal stressPa
n_refReference stress for 90% closurePa
τShear stressPa
ϕAngle°
ωDynamic viscositykg m−1 s−1
Subscript
0Initial
circCirculation
HMaximum horizontal principal stress
hMinimum horizontal principal stress
recRecovered
sScaling
stimStimulation
thThermal
VVertical principal stress
Abbreviation
HSDHydrothermal spallation drilling

References

  1. Sustainable Development Goals. Available online: https://www.un.org/sustainabledevelopment/energy/ (accessed on 25 January 2021).
  2. Off-Grid Renewable Energy Solutions to Expand Electricity Access: An Opportunity Not to Be Missed. Available online: https://www.irena.org/-/media/Files/IRENA/Agency/Publication/2019/Jan/IRENA_Off-grid_RE_Access_2019.pdf (accessed on 22 February 2021).
  3. Market Snapshot: Overcoming the Challenges of Powering Canada’s Off-Grid Communities. Available online: https://www.cer-rec.gc.ca/en/data-analysis/energy-markets/market-snapshots/2018/market-snapshot-overcoming-challenges-powering-canadas-off-grid-communities.html (accessed on 22 February 2021).
  4. Energy Access—The Canadian Context. Available online: https://wise.uwaterloo.ca/documents/events/marketing-materials/ae4h/ae4h-blogs/energy_access_canadian_context_infographic_spread_and_section_openac (accessed on 18 November 2019).
  5. Khan, F.I.; Hawbolt, K.; Iqbal, M.T. Life cycle analysis of wind-fuel cell integrated system. Renew. Energy 2005, 30, 157–177. [Google Scholar] [CrossRef]
  6. Ranjitkar, G.; Huang, J.; Tung, T. Application of micro-hydropower technology for remote regions. In Proceedings of the IEEE EIC Climate Change Conference, Ottawa, ON, Canada, 10–12 May 2006; IEEE: Ottawa, ON, Canada, 2006. [Google Scholar]
  7. Thompson, S.; Duggirala, B. The feasibility of renewable energies at an off-grid community in Canada. Renew. Sustain. Energy Rev. 2009, 13, 2740–2745. [Google Scholar] [CrossRef]
  8. Ibrahim, H.; Younès, R.; Ilinca, A.; Ramdenee, D.; Dimitrova, M.; Perron, J.; Adegnon, M.; Boulay, D.; Arbez, C. Potential of a hybrid wind-diesel-compressed air system for Nordic remote Canadian areas. Energy Procedia 2011, 6, 795–804. [Google Scholar] [CrossRef] [Green Version]
  9. Guo, L.; Yu, Z.; Wang, C.; Li, F.; Schiettekatte, J.; Deslauriers, J.-C.; Bai, L. Optimal design of battery energy storage system for a wind-diesel off-grid power system in a remote Canadian community. IET Gener. Transm. Distrib. 2016, 10, 608–616. [Google Scholar] [CrossRef]
  10. Stephen, J.D.; Mabee, W.E.; Pribowo, A.; Pledger, S.; Hart, R.; Tallio, S.; Bull, G.Q. Biomass for residential and commercial heating in a remote Canadian aboriginal community. Renew. Energy 2016, 86, 563–575. [Google Scholar] [CrossRef] [Green Version]
  11. Belzile, P.; Comeau, F.-A.; Raymond, J.; Lamarche, L.; Carreau, M. Arctic climate horizontal ground-coupled heat pump. GRC Trans. 2017, 41, 1958–1978. [Google Scholar]
  12. McFarlan, A. Techno-economic assessment of pathways for electricity generation in northern remote communities in Canada using methanol and dimethyl ether to replace diesel. Renew. Sustain. Energy Rev. 2018, 90, 863–876. [Google Scholar] [CrossRef]
  13. Kanzari, I. Évaluation du Potentiel des Pompes à Chaleur Géothermique Pour la Communauté Nordique de Kuujjuaq. Master’s Thesis, Institut National de la Recherche Scientifique (INRS), Québec, QC, Canada, 2019. [Google Scholar]
  14. Yan, C.; Rousse, D.; Glaus, M. Multi-criteria decision analysis ranking alternative heating systems for remote communities in Nunavik. J. Clean. Prod. 2019, 208, 1488–1497. [Google Scholar] [CrossRef]
  15. Gunawan, E.; Giordano, N.; Jensson, P.; Newson, J.; Raymond, J. Alternative heating systems for northern remote communities: Techno-economic analysis of ground-coupled heat pumps in Kuujjuaq, Nunavik, Canada. Renew. Energy 2020, 147, 1540–1553. [Google Scholar] [CrossRef]
  16. Mahbaz, S.B.; Dehghani-Sanij, A.R.; Dusseault, M.B.; Nathawani, J.S. Enhanced and integrated geothermal systems for sustainable development of Canada’s northern communities. Sustain. Energy Technol. Assess. 2020, 37, 100565. [Google Scholar] [CrossRef]
  17. Quitoras, M.R.; Campana, P.E.; Crawford, C. Exploring electricity generation alternatives for Canadian Arctic communities using a multi-objective genetic algorithm approach. Energy Convers. Manag. 2020, 2010, 112471. [Google Scholar] [CrossRef]
  18. Majorowicz, J.; Grasby, S.E. Heat flow, depth-temperature variations and stored thermal energy for enhanced geothermal systems in Canada. J. Geophys. Eng. 2010, 7, 232–241. [Google Scholar] [CrossRef] [Green Version]
  19. Majorowicz, J.; Grasby, S.E. High potential regions for enhanced geothermal systems in Canada. Nat. Resour. Res. 2010, 19, 177–188. [Google Scholar] [CrossRef]
  20. Kunkel, T.; Ghomshei, M.; Ellis, R. Geothermal energy as an indigenous alternative energy source in British Columbia. J. Ecosyst. Manag. 2012, 13, 1–16. [Google Scholar]
  21. Grasby, S.E.; Majorowicz, J.; McCune, G. Geothermal Energy Potential for Northern Communities; Report No.: Open File 7350; Geothermal Survey of Canada: Ottawa, ON, Canada, 2013; 50p. [Google Scholar]
  22. Walsh, W. Geothermal resource assessment of the Clarke Lake Gas Field, Fort Nelson, British Columbia. Bull. Can. Pet. Geol. 2013, 61, 241–251. [Google Scholar] [CrossRef]
  23. Majorowicz, J.; Grasby, S.E. Geothermal energy for northern Canada: Is it economical? Nat. Resour. Res. 2014, 23, 159–173. [Google Scholar] [CrossRef]
  24. Majorowicz, J.; Minea, V. Shallow and deep geothermal energy potential in low heat flow/cold climate environment: Northern Quebec, Canada, case study. Environ. Earth Sci. 2015, 74, 5233–5244. [Google Scholar] [CrossRef]
  25. Minnick, M.; Shewfelt, D.; Hickson, C.; Majorowicz, J.; Rowe, T. Nunavut Geothermal Feasibility Study; Report No.: Topical Report RSI-2828; Qulliq Energy Corporation: Iqaluit, NT, Canada, 2018; 42p. [Google Scholar]
  26. Majorowicz, J.A.; Grasby, S.E. Heat transition for major communities supported by geothermal energy development of the Alberta Basin, Canada. Geothermics 2020, 88, 101883. [Google Scholar] [CrossRef]
  27. Majorowicz, J.; Grasby, S.E. Deep geothermal heating potential for the communities of the Western Canadian Sedimentary Basin. Energies 2021, 14, 706. [Google Scholar] [CrossRef]
  28. Miranda, M.M.; Raymond, J.; Dezayes, C. Uncertainty and risk evaluation of deep geothermal energy source for heat production and electricity generation in remote northern regions. Energies 2020, 13, 4221. [Google Scholar] [CrossRef]
  29. Grasby, S.E.; Allen, D.M.; Bell, S.; Chen, Z.; Ferguson, G.; Jessop, A.; Kelman, M.; Ko, M.; Majorowicz, J.; Moore, M.; et al. Geothermal Energy Resource Potential of Canada; Report No.: Open File 6914; Geological Survey of Canada: Ottawa, ON, Canada, 2012; 322p. [Google Scholar]
  30. Hofmann, H.; Weides, S.; Babadagli, T.; Zimmermann, G.; Moeck, I.; Majorowicz, J.; Unsworth, M. Potential for enhanced geothermal systems in Alberta, Canada. Energy 2014, 69, 578–591. [Google Scholar] [CrossRef]
  31. Hofmann, H.; Babadagli, T.; Yoon, J.S.; Blöcher, G.; Zimmermann, G. A hybrid discrete/finite element modeling study of complex hydraulic fracture development for enhanced geothermal systems (EGS) in granitic basements. Geothermics 2016, 64, 362–381. [Google Scholar] [CrossRef]
  32. Kazemi, A.R.; Mahbaz, S.B.; Dehghani-Sanij, A.R.; Dusseault, M.B.; Fraser, R. Performance evaluation of an enhanced geothermal system in the Western Canada Sedimentary Basin. Renew. Sustain. Energy Rev. 2019, 113, 109278. [Google Scholar] [CrossRef]
  33. Gidley, J.L.; Holditch, S.A.; Nierode, D.A.; Veatch, R.W., Jr. Recent Advances in Hydraulic Fracturing; Society of Petroleum Engineers: Richardson, TX, USA, 1989; 461p. [Google Scholar]
  34. Brown, D.W.; Duchane, D.V.; Heiken, G.; Hriscu, V.T. Mining the Earth’s Heat: Hot Dry Rock Geothermal Energy; Springer: Berlin, Germany, 2012; 669p. [Google Scholar]
  35. Tester, J.W.; Anderson, B.J.; Batchelor, A.S.; Blackwell, D.D.; DiPippo, R.; Drake, E.M.; Garnish, J.; Livesay, B.; Moore, M.C.; Nichols, K.; et al. The Future of Geothermal Energy: Impact of Enhanced Geothermal Systems (EGS) on the United States in the 21st Century; MIT: Idaho Falls, ID, USA; 372p.
  36. Breede, K.; Dzebisashvili, K.; Liu, X.; Falcone, G. A systematic review of enhanced (or engineered) geothermal systems: Past, present and future. Geotherm. Energy 2013, 1, 4. [Google Scholar] [CrossRef] [Green Version]
  37. Xie, L.; Min, K.-B.; Song, Y. Observations of hydraulic stimulations in seven enhanced geothermal system projects. Renew. Energy 2015, 79, 56–65. [Google Scholar] [CrossRef]
  38. Olasolo, P.; Juárez, M.C.; Olasolo, J.; Morales, M.P.; Valdani, D. Enhanced geothermal systems (EGS): A review. Renew. Sustain. Energy Rev. 2016, 56, 133–144. [Google Scholar] [CrossRef]
  39. Lu, S.-M. A global review of enhanced geothermal system (EGS). Renew. Sustain. Energy Rev. 2018, 81, 2902–2921. [Google Scholar] [CrossRef]
  40. Kumari, W.G.P.; Ranjith, P.G. Sustainable development of enhanced geothermal systems based on geotechnical research—A review. Earth-Sci. Rev. 2019, 199, 102955. [Google Scholar] [CrossRef]
  41. Kelkar, S.; Gabriel, G.W.; Rehfeldt, K. Lessons learned from the pioneering hot dry rock project at Fenton Hill, USA. Geothermics 2016, 63, 5–14. [Google Scholar] [CrossRef] [Green Version]
  42. Parker, R. The rosemanowes HDR project 1983-1991. Geothermics 1999, 28, 603–615. [Google Scholar] [CrossRef]
  43. 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.; Willis-Richards, J. The performance and characteristics of the experimental hot dry rock geothermal reservoir at Rosemanowes, Cornwall (1985–1988). Geothermics 1994, 23, 73–109. [Google Scholar] [CrossRef]
  44. Genter, A.; Evans, K.; Cuenot, N.; Fritsch, D.; Sanjuan, B. Contribution of the exploration of deep crystalline fractured reservoir of Soultz to the knowledge of enhanced geothermal systems (EGS). Comptes Rendus Geosci. 2010, 342, 502–516. [Google Scholar] [CrossRef]
  45. Comeau, F.-A.; Raymond, J.; Malo, M.; Dezayes, C.; Carreau, M. Geothermal potential of northern Quebec: A regional assessment. GRC Trans. 2017, 41, 1076–1094. [Google Scholar]
  46. Census Profile, 2016 Census. Available online: https://www12.statcan.gc.ca/census-recensement/2016/dp-pd/prof/details/page.cfm?Lang=E&Geo1=POPC&Code1=1466&Geo2=PR&Code2=24&SearchText=Kuujjuaq&SearchType=Begins&SearchPR=01&B1=All&GeoLevel=PR&GeoCode=1466&TABID=1&type=0 (accessed on 6 October 2020).
  47. Carte Interactive. Available online: https://sigeom.mines.gouv.qc.ca/signet/classes/I1108_afchCarteIntr (accessed on 5 May 2019).
  48. Miranda, M.M.; Giordano, N.; Raymond, J.; Pereira, A.J.S.C.; Dezayes, C. Thermophysical properties of surficial rocks: A tool to characterize geothermal resources of remote northern regions. Geotherm. Energy 2020, 8, 4. [Google Scholar] [CrossRef]
  49. Miranda, M.M.; Marquez, M.I.V.; Raymond, J.; Dezayes, C. A numerical approach to infer terrestrial heat flux from shallow temperature profiles in remote northern regions. Geothermics 2021, 93, 102064. [Google Scholar] [CrossRef]
  50. Miranda, M.M.; Dezayes, C.; Giordano, N.; Kanzari, I.; Raymond, J.; Carvalho, J. Fracture network characterization as input for geothermal energy research: Preliminary data from Kuujjuaq, northern Quebec, Canada. In Proceedings of the 43rd Workshop on Geothermal Reservoir Engineering, Stanford, CA, USA, 12–14 February 2018. [Google Scholar]
  51. Miranda, M.M. Assessment of the Deep Geothermal Energy Source Potential in Remote Northern Regions: A Study Undertaken in the Subarctic Off-Grid Community of Kuujjuaq, Nunavik, Canada. Ph.D. Thesis, Institut National de la Recherche Scientifique (INRS), Québec, QC, Canada, 2021. [Google Scholar]
  52. Jing, Z.; Willis-Richards, J.; Watanabe, K.; Hashida, T. A three-dimensional stochastic rock mechanics model of engineered geothermal systems in fractured crystalline rock. J. Geophys. Res. 2000, 150, 23663–23679. [Google Scholar] [CrossRef]
  53. Murphy, H.; Dash, Z. The shocking behavior of fluid flow in deformable joints. GRC Transactions 1985, 9, 115–118. [Google Scholar]
  54. Lee, H.S.; Cho, T.F. Hydraulic characteristics of rough fractures in linear flow under normal and shear load. Rock Mech. Rock Eng. 2002, 35, 299–318. [Google Scholar] [CrossRef]
  55. Goodman, R.E. Methods of Geological Engineering in Discontinuous Rocks; West Publishing Company: New York, NY, USA, 1976; 484p. [Google Scholar]
  56. Bandis, S.C.; Lumsden, A.C.; Barton, N.R. Fundamentals of rock joint deformation. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1983, 20, 249–268. [Google Scholar] [CrossRef]
  57. Barton, N.; Quadros, E.F. Joint aperture and roughness in the prediction of flow and groutability of rock masses. Int. J. Rock Mech. Min. Sci. 1997, 34, 252.e1–252.e14. [Google Scholar] [CrossRef]
  58. Keller, A.A. High resolution cat imaging of fractures in consolidated materials. Int. J. Rock Mech. Min. Sci. 1997, 34, 155.e1–155.e16. [Google Scholar] [CrossRef]
  59. Keller, A.A. High resolutions, non-destructive measurement and characterization of fracture apertures. Int. J. Rock Mech. Min. Sci. 1998, 35, 1037–1050. [Google Scholar] [CrossRef]
  60. Fomin, S.; Jing, Z.; Hashida, T. The effect of thermal, chemical, hydrological, and mechanical factors on water/rock interaction in HDR geothermal systems. In Engineering Properties of Rocks; Zhang, L., Ed.; Elsevier: Berlin, Germany, 2004; Volume 2, pp. 649–654. [Google Scholar]
  61. Simard, M.; Lafrance, I.; Hammouche, H.; Legouix, C. 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, QC, Canada, 2013; 62p. [Google Scholar]
  62. Lanyon, G.W.; Batchelor, A.S.; Ledingham, P. Results from a discrete fracture network model of a hot dry rock system. In Proceedings of the 18th Workshop on Geothermal Reservoir Engineering, Stanford, CA, US, 26–28 January 1993. [Google Scholar]
  63. Evans, K.F.; Moriya, H.; Niitsuma, H.; Jones, R.H.; Phillips, W.S.; Genter, A.; Sausse, J.; Jung, R.; Baria, R. Microseismicity and permeability enhancement of hydrogeologic structures during massive fluid injections into granite at 3 km depth at the Soultz HDR site. Geophys. J. Int. 2005, 160, 388–412. [Google Scholar] [CrossRef]
  64. Arjang, B.; Herget, G. In situ ground stresses in the Canadian hardrock mines: An update. Int. J. Rock Mech. Min. Sci. 1997, 34, 15. [Google Scholar] [CrossRef]
  65. Alejano, L.R.; González, J.; Muralha, J. Comparison of different techniques of tilt testing and basic friction angle variability assessment. Rock Mech. Rock Eng. 2012, 45, 1023–1035. [Google Scholar] [CrossRef]
  66. Achtziger-Zupancic, P.; Loew, S.; Mariéthoz, G. A new global database to improve predictions of permeability distribution in crystalline rocks at site scale. J. Geophys. Res. Solid Earth 2017, 122, 3513–3539. [Google Scholar] [CrossRef]
  67. Bottomley, D.J. A Review of Theories on the Origins of Saline Waters and Brines in the Canadian Precambrian Shield; Atomic Energy Control Board: Ottawa, ON, Canada, 1996; 32p. [Google Scholar]
  68. Mouchot, J.; Genter, A.; Cuenot, N.; Scheiber, J.; Seibel, O.; Bosia, C.; Ravier, G. First year of operation from EGS geothermal plants in Alsace, France: Scaling issues. In Proceedings of the 43rd Workshop on Geothermal Reservoir Engineering, Stanford, CA, USA, 12–14 February 2018. [Google Scholar]
  69. Frick, S.; van Wees, J.D.; Kaltschmitt, M.; Schröder, G. Economic performance and environmental assessment. In Geothermal Energy Systems: Exploration, Development, and Utilization; Huenges, E., Ed.; Wiley: Weinheim, Germany, 2010; pp. 373–421. [Google Scholar]
  70. Beckers, K.F.; McCabe, K. GEOPHIRES V2.0: Updated geothermal techno-economic simulation tool. Geotherm. Energy 2019, 7, 5. [Google Scholar] [CrossRef]
  71. Limberger, J.; Calcagno, P.; Manzella, A.; Trumpy, E.; Boxem, T.; Pluymaekers, M.P.D.; van Wees, J.-D. Assessing the prospective resource base for enhanced geothermal systems in Europe. Geotherm. Energy Sci. 2014, 2, 55–71. [Google Scholar] [CrossRef] [Green Version]
  72. Sanyal, S.K.; Morrow, J.W.; Butler, S.J.; Robertson-Tait, A. Cost of electricity from enhanced geothermal systems. In Proceedings of the 32nd Workshop on Geothermal Reservoir Engineering, Stanford, CA, USA, 22–24 January 2007. [Google Scholar]
  73. Lance, J. Personal Communication; FCNQ: Quebec, QC, Canada, 2021. [Google Scholar]
  74. @RISK. Available online: https://www.palisade.com/risk/default.asp (accessed on 18 June 2020).
  75. Snow, D.T. Anisotropic permeability of fracture media. Water Resour. Res. 1969, 5, 1273–1289. [Google Scholar] [CrossRef]
  76. Willis-Richards, J.; Watanabe, K.; Takahashi, H. Progress toward a stochastic rock mechanics model of engineered geothermal systems. J. Geophys. Res. 1996, 101, 17481–17496. [Google Scholar] [CrossRef]
  77. Eshelby, J.D. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. Lond. A 1957, 241, 376–396. [Google Scholar]
  78. Beardsmore, G.R.; Rybach, L.; Blackwell, D.; Baron, C. A protocol for estimating and mapping global EGS potential. GRC Trans. 2010, 34, 301–312. [Google Scholar]
  79. Giordano, N.; Kanzari, I.; Miranda, M.M.; Dezayes, C.; Raymond, J. Underground thermal energy storage in subarctic climates: A feasibility study conducted in Kuujjuaq (QC, Canada). In Proceedings of the IGSHPA, Stockholm, Sweden, 18–20 September 2018. [Google Scholar]
  80. Karanasios, K.; Parker, P. Tracking the transition to renewable electricity in remote indigenous communities in Canada. Energy Policy 2018, 118, 169–181. [Google Scholar] [CrossRef]
  81. Zweifel, P.; Praktiknjo, A.; Erdmann, G. Energy Economics: Theory and Applications; Springer: Berlin, Germany, 2017; 333p. [Google Scholar]
  82. Glassley, W.E. Geothermal Energy: Renewable Energy and the Environment; CRC Press: Boca Raton, FL, USA, 2010; 285p. [Google Scholar]
  83. McClure, M.W.; Horne, R.N. An investigation of stimulation mechanisms in Enhanced Geothermal Systems. Int. J. Rock Mech. Min. Sci. 2014, 72, 242–260. [Google Scholar] [CrossRef] [Green Version]
  84. Evans, K.; Cornet, F.; Hayashi, K.; Hashida, T.; Ito, T.; Matsuki, K.; Wallroth, T. Stress and rock mechanics issues of relevance to HDR/HWR engineered geothermal systems: A review of developments during the past 15 years. Geothermics 1999, 28, 455–474. [Google Scholar] [CrossRef]
  85. Jalali, M.; Gischig, V.; Doetsch, J.; Näf, R.; Krietsch, H.; Klepikova, M.; Amann, F.; Giardini, D. Transmissivity changes and microseismicity induced by small-scale hydraulic fracturing tests in crystalline rock. Geophys. Res. Lett. 2018, 45, 2265–2273. [Google Scholar] [CrossRef] [Green Version]
  86. Norbeck, J.H.; McClure, M.W.; Horne, R.N. Field observations at the Fenton Hill enhanced geothermal system test site support mixed-mechanism stimulation. Geothermics 2018, 74, 135–149. [Google Scholar] [CrossRef]
  87. Jing, Y.; Jing, Z.; Willis-Richards, J.; Hashida, T. A simple 3-D thermoelastic model for assessment of the long-term performance of the Hijiori deep geothermal reservoir. J. Volcanol. Geotherm. Res. 2014, 269, 14–22. [Google Scholar] [CrossRef]
  88. Gilman, J.R.; Ozgen, C. Reservoir Simulation: History Matching and Forecasting; Society of Petroleum Engineers: Richardson, TX, USA, 2013; 129p. [Google Scholar]
  89. Pine, R.J.; Batchelor, A.S. Downward migration of shearing in jointed rock during hydraulic injections. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1984, 21, 249–263. [Google Scholar] [CrossRef]
  90. Sicking, C.; Malin, P. Fracture seismic: Mapping subsurface connectivity. Geosciences 2019, 9, 508. [Google Scholar] [CrossRef] [Green Version]
  91. Zhang, S.; Yin, S.; Yuan, Y. Estimation of fracture stiffness, in situ stresses, and elastic parameters of naturally fractured geothermal reservoirs. Int. J. Geomech. 2015, 15, 04014033. [Google Scholar] [CrossRef]
  92. Chen, Z.; Narayan, S.P.; Yang, Z.; Rahman, S.S. An experimental investigation of hydraulic behaviour of fractures and joints in granitic rock. Int. J. Rock Mech. Min. Sci. 2000, 37, 1061–1071. [Google Scholar] [CrossRef]
  93. Barton, N. Shear strength criteria for rock, rock joints, rockfill and rock masses: Problems and some solutions. J. Rock Mech. Geotech. Eng. 2013, 5, 249–261. [Google Scholar] [CrossRef] [Green Version]
  94. Frash, L.P. Laboratory-Scale Study of Hydraulic Fracturing in Heterogeneous Media for Enhanced Geothermal Systems and General Well Stimulation. Ph.D. Thesis, Colorado School of Mines, Golden, CO, USA, 2014. [Google Scholar]
  95. Mao, R.; Feng, Z.; Liu, Z.; Zhao, Y. Laboratory hydraulic fracturing test on large-scale pre-cracked granite specimens. J. Nat. Gas Sci. Eng. 2017, 44, 278–286. [Google Scholar] [CrossRef]
  96. Deb, P.; Düber, S.; Carducci, C.G.C.; Clauser, C. Laboratory-scale hydraulic fracturing dataset for benchmarking of enhanced geothermal system simulation tools. Sci. Data 2020, 7, 220. [Google Scholar] [CrossRef]
  97. Hu, L.; Ghassemi, A.; Pritchett, J.; Garg, S. Characterization of laboratory-scale hydraulic fracturing for EGS. Geothermics 2020, 83, 101706. [Google Scholar] [CrossRef]
  98. Deb, P.; Salimzadeh, S.; Vogler, D.; Düber, S.; Clauser, C.; Settgast, R.R. Verification of coupled hydraulic fracturing simulators using laboratory-scale experiments. Rock Mech. Rock Eng. 2021, 54, 2881–2902. [Google Scholar] [CrossRef]
  99. Weydt, L.M.; Ramírez-Guzmán, Á.A.; Pola, A.; Lepillier, B.; Kummerow, J.; Mandrone, G.; Comina, C.; Deb, P.; Norini, G.; Gonzalez-Partida, E.; et al. Petrophysical and mechanical rock property database of the Los Humeros and Acoculco geothermal fields (Mexico). Earth Syst. Sci. Data 2021, 13, 571–598. [Google Scholar] [CrossRef]
  100. Koelbel, T.; Genter, A. Enhanced geothermal systems: The soultz-sous-forêts project. In Towards 100% Renewable Energy: Techniques, Costs and Regional Case-Studies; Uyar, T.S., Ed.; Springer: Cham, Switzerland, 2017; pp. 243–248. [Google Scholar]
  101. Jing, Z.; Watanabe, K.; Willis-Richards, J.; Hashida, T. A 3-D water/rock chemical interaction model for prediction of HDR/HWR geothermal reservoir performance. Geothermics 2002, 31, 1–28. [Google Scholar] [CrossRef]
  102. Farquharson, J.I.; Kushnir, A.R.L.; Wild, B.; Baud, P. Physical property evolution of granite during experimental chemical stimulations. Geotherm. Energy 2020, 8, 14. [Google Scholar] [CrossRef]
  103. Saadat, A.; Frick, S.; Kranz, S.; Regenspurg, S. Energetic use of EGS reservoirs. In Geothermal Energy Systems: Exploration, Development, and Utilization; Huenges, E., Ed.; Wiley: Weinheim, Germany, 2010; pp. 303–372. [Google Scholar]
  104. Grant, M.A.; Bixley, P.F. Geothermal Reservoir Engineering; Academic Press: Oxford, UK, 2011; 379p. [Google Scholar]
  105. Brown, D. A Hot Dry Rock geothermal energy concept utilizing supercritical CO2 instead of water. In Proceedings of the 25th Workshop on Geothermal Reservoir Engineering, Stanford, CA, USA, 24–26 January 2000. [Google Scholar]
  106. Pruess, K. Enhanced geothermal systems (EGS) using CO2 as working fluid—A novel approach for generating renewable energy with simultaneous sequestration of carbon. Geothermics 2006, 35, 351–367. [Google Scholar] [CrossRef] [Green Version]
  107. Pruess, K. On production behavior of enhanced geothermal systems with CO2 as working fluid. Energy Convers. Manag. 2008, 49, 1446–1454. [Google Scholar] [CrossRef] [Green Version]
  108. Guo, T.; Gong, F.; Wang, X.; Lin, Q.; Qu, Z.; Zhang, W. Performance of enhanced geothermal system (EGS) in fractured geothermal reservoirs with CO2 as working fluid. Appl. Therm. Eng. 2019, 152, 215–230. [Google Scholar] [CrossRef]
  109. Karanasios, K.; Parker, P. Recent developments in renewable energy in remote aboriginal communities, Quebec, Canada. Pap. Can. Econ. Dev. 2016, 16, 98–108. [Google Scholar]
  110. Bonalumi, D. Enhanced geothermal systems with captured CO2. Energy Procedia 2018, 148, 744–750. [Google Scholar] [CrossRef]
  111. Baroudi, H.A.; Awoyomi, A.; Patchigolla, K.; Jonnalagadda, K.; Anthony, E.J. A review of large-scale CO2 shipping and marine emissions management for carbon capture, utilization and storage. Appl. Energy 2021, 287, 116510. [Google Scholar] [CrossRef]
  112. Olasolo, P.; Juárez, M.C.; Morales, M.P.; Olasolo, A.; Agius, M.R. Analysis of working fluids applicable in Enhanced Geothermal Systems: Nitrous oxide as an alternative working fluid. Energy 2018, 157, 150–161. [Google Scholar] [CrossRef]
  113. Augustine, C.R. Updated US Geothermal Supply Characterization and Representation for Market Penetration Model Input; Report No.: NREL/TP-6A20-47459; NREL: Golden, CO, USA, 2011; 103p. [Google Scholar]
  114. Richard, M.-A. Production D’électricité avec des Systèmes Géothermiques Stimulés au Québec: Analyse des Résultats d’un Outil de Simulation; Report No.: IREQ-2016-0001; Hydro Québec Institut de Recherche: Varennes, QC, Canada, 2016; 164p. [Google Scholar]
  115. Kaczmarczyk, M.; Tomaszewska, B.; Operacz, A. Sustainable utilization of low entalphy geothermal resources to electricity generation through a cascade system. Energies 2020, 13, 2495. [Google Scholar] [CrossRef]
  116. Operacz, A.; Chowaniec, J. Perspectives of geothermal water use in the Podhale Basin according to geothermal step distribution. Geol. Geophys. Environ. 2018, 44, 379–389. [Google Scholar] [CrossRef] [Green Version]
  117. Sowizdzal, A.; Chmielowska, A.; Tomaszewska, B.; Operacz, A.; Chowaniec, J. Could geothermal water and energy use improve living conditions? Environmental effects from Polans. Arch. Environ. Prot. 2019, 45, 109–118. [Google Scholar]
  118. Giordano, N.; Raymond, J. Alternative and sustainable heat production for drinking water needs in a subarctic climate (Nunavik, Canada): Borehole thermal energy storage to reduce fossil fuel dependency in off-grid communities. Appl. Energy 2019, 252, 113463. [Google Scholar] [CrossRef]
  119. Javadi, H.; Ajarostaghi, S.S.M.; Rosen, M.A.; Pourfallah, M. Performance of ground heat exchangers: A comprehensive review of recent advances. Energy 2019, 178, 207–233. [Google Scholar] [CrossRef]
  120. Javadi, H.; Ajarostaghi, S.S.M.; Rosen, M.A.; Pourfallah, M. A comprehensive review of backfill materials and their effect on ground heat exchanger performance. Sustainability 2018, 10, 4486. [Google Scholar] [CrossRef] [Green Version]
  121. Javadi, H.; Ajarostaghi, S.S.M.; Pourfallah, M.; Zaboli, M. Performance analysis of helical ground heat exchangers with different configurations. Appl. Therm. Eng. 2019, 154, 24–36. [Google Scholar] [CrossRef]
  122. Bremen, L.V. Large-scale variability of weather dependent renewable energy sources. In Management of Weather and Climate Risk in the Energy Industry; Troccoli, A., Ed.; Springer: Dordrecht, The Netherlands, 2010; pp. 186–206. [Google Scholar]
  123. Cherry, J.E.; Knapp, C.; Trainor, S.; Ray, A.J.; Tedesche, M.; Walker, S. Planning for climate change impacts on hydropower in the Far North. Hydrol. Earth Syst. Sci. 2017, 21, 133–151. [Google Scholar] [CrossRef] [Green Version]
  124. Reinecker, J.; Gutmanis, J.; Foxford, A.; Cotton, L.; Dalby, C.; Law, R. Geothermal exploration and reservoir modelling of the United Downs deep geothermal project, Cornwall (UK). Geothermics 2021, 97, 102226. [Google Scholar] [CrossRef]
  125. Somma, R.; Blessent, D.; Raymond, J.; Constance, M.; Cotton, L.; Natale, G.; Fedele, A.; Jurado, M.J.; Marcia, K.; Miranda, M.M.; et al. Review of recent drilling projects in unconventional geothermal resources at Campi Flegrei Caldera, Cornubian Batholith, and Williston Sedimentary Basin. Energies 2021, 14, 3306. [Google Scholar] [CrossRef]
  126. Operacz, A.; Bielec, B.; Tomaszewska, B.; Kaczmarczyk, M. Physicochemical composition variability and hydraulic conditions in a geothermal borehole—The latest study in Podhale Basin, Poland. Energies 2020, 13, 3882. [Google Scholar] [CrossRef]
Figure 1. (a) Distribution of heat flow data available and (b) geothermal heat flow density for northern Québec. In (a) pink—Superior Province, green—Churchill Province, orange—Grenville Province, yellow—Appalachian Province, dark blue—Hudson Bay Platform and light blue—St. Lawrence Lowlands. Redrawn from Majorowicz and Minea [24] and Comeau et al. [45].
Figure 1. (a) Distribution of heat flow data available and (b) geothermal heat flow density for northern Québec. In (a) pink—Superior Province, green—Churchill Province, orange—Grenville Province, yellow—Appalachian Province, dark blue—Hudson Bay Platform and light blue—St. Lawrence Lowlands. Redrawn from Majorowicz and Minea [24] and Comeau et al. [45].
Water 13 03526 g001
Figure 2. (a) Geographical location of Kuujjuaq and the remaining communities in Nunavik; (b) geological setting of the study area. LP—Lac Pingiajjulik fault, LG—Lac Gabriel fault. Adapted from SIGÉOM [47] and Miranda et al. [48].
Figure 2. (a) Geographical location of Kuujjuaq and the remaining communities in Nunavik; (b) geological setting of the study area. LP—Lac Pingiajjulik fault, LG—Lac Gabriel fault. Adapted from SIGÉOM [47] and Miranda et al. [48].
Water 13 03526 g002
Figure 3. Average daily temperature and heating load profile of a typical residential dwelling in Kuujjuaq, redrawn from Gunawan et al. [15].
Figure 3. Average daily temperature and heating load profile of a typical residential dwelling in Kuujjuaq, redrawn from Gunawan et al. [15].
Water 13 03526 g003
Figure 4. (a) Measured (solid line) and simulated (dashed line) temperature profile near the community of Kuujjuaq. The simulated profile corresponds to a heat flux of 41.6 mW m−2 at 10 km depth (reprinted with permission from Miranda et al. [49], 2021 Elsevier Ltd.); (b) magnitude and orientation of the principal stresses. The error bar represents the range of values inferred. σH—maximum horizontal principal stress, σh—minimum horizontal principal stress, σV—vertical principal stress, Ppore—in situ fluid pressure.
Figure 4. (a) Measured (solid line) and simulated (dashed line) temperature profile near the community of Kuujjuaq. The simulated profile corresponds to a heat flux of 41.6 mW m−2 at 10 km depth (reprinted with permission from Miranda et al. [49], 2021 Elsevier Ltd.); (b) magnitude and orientation of the principal stresses. The error bar represents the range of values inferred. σH—maximum horizontal principal stress, σh—minimum horizontal principal stress, σV—vertical principal stress, Ppore—in situ fluid pressure.
Water 13 03526 g004
Figure 5. Top view slice cutting the center of the model and illustrating the different wells configuration studied. In (A) the wells are placed perpendicular to the fault trace and parallel to the direction of the maximum horizontal principal stress. In (B) the wells are placed parallel to the fault trace and minimum horizontal principal stress. Dots—vertical wells, dashed line—Lac Pingiajjulik fault trace, ellipsoids—stimulation volume per well tried (larger ellipsoid, V = 0.4 km3; smaller ellipsoid, V = 0.2 km3; the sum of the stimulated volumes in each well were defined to correspond to 20% and 10%, respectively, of the total model volume).
Figure 5. Top view slice cutting the center of the model and illustrating the different wells configuration studied. In (A) the wells are placed perpendicular to the fault trace and parallel to the direction of the maximum horizontal principal stress. In (B) the wells are placed parallel to the fault trace and minimum horizontal principal stress. Dots—vertical wells, dashed line—Lac Pingiajjulik fault trace, ellipsoids—stimulation volume per well tried (larger ellipsoid, V = 0.4 km3; smaller ellipsoid, V = 0.2 km3; the sum of the stimulated volumes in each well were defined to correspond to 20% and 10%, respectively, of the total model volume).
Water 13 03526 g005
Figure 6. Top view of an 80 per 80 m section of the total generated fracture network model: (a) base-case fracture network; (b) fracture network with fracture lengths increased by a factor of 10; (c) fracture network with fractal dimension decreased by a factor of 2; (d) Wulff projection stereo plot of field fracture data. Adapted from Miranda et al. [51].
Figure 6. Top view of an 80 per 80 m section of the total generated fracture network model: (a) base-case fracture network; (b) fracture network with fracture lengths increased by a factor of 10; (c) fracture network with fractal dimension decreased by a factor of 2; (d) Wulff projection stereo plot of field fracture data. Adapted from Miranda et al. [51].
Water 13 03526 g006
Figure 7. Boxplot of initial fracture aperture with whiskers from minimum to maximum. The reader is referred to Section 4.2 for further details on the fracture network and to Table 1 for further information on the properties of the medium.
Figure 7. Boxplot of initial fracture aperture with whiskers from minimum to maximum. The reader is referred to Section 4.2 for further details on the fracture network and to Table 1 for further information on the properties of the medium.
Water 13 03526 g007
Figure 8. Boxplot of fracture shear stiffness with whiskers from minimum to maximum. The reader is referred to Section 4.2 for further details on the fracture network and to Table 1 for further information on the properties of the medium.
Figure 8. Boxplot of fracture shear stiffness with whiskers from minimum to maximum. The reader is referred to Section 4.2 for further details on the fracture network and to Table 1 for further information on the properties of the medium.
Water 13 03526 g008
Figure 9. Boxplot of fracture shear displacement with whiskers from minimum to maximum: (a) stimulation pressure of 20 MPa; (b) stimulation pressure of 30 MPa; (c) stimulation pressure of 40 MPa; (d) stimulation pressure of 50 MPa. The reader is referred to Section 4.2 for further details on the fracture network and to Table 1 for further information on the properties of the medium.
Figure 9. Boxplot of fracture shear displacement with whiskers from minimum to maximum: (a) stimulation pressure of 20 MPa; (b) stimulation pressure of 30 MPa; (c) stimulation pressure of 40 MPa; (d) stimulation pressure of 50 MPa. The reader is referred to Section 4.2 for further details on the fracture network and to Table 1 for further information on the properties of the medium.
Water 13 03526 g009
Figure 10. Sheared fractures grouped by sets for different stimulation fluid pressure. The reader is referred to Section 4.2 for further details on the fracture network and to Table 1 for further information on the properties of the medium. Relative frequency, in this context, corresponds to the percentage of fractures that slipped in each set from the total of fractures that slip.
Figure 10. Sheared fractures grouped by sets for different stimulation fluid pressure. The reader is referred to Section 4.2 for further details on the fracture network and to Table 1 for further information on the properties of the medium. Relative frequency, in this context, corresponds to the percentage of fractures that slipped in each set from the total of fractures that slip.
Water 13 03526 g010
Figure 11. Heating energy produced (color lines) and projected demand (grey polygon). Lower bound of grey polygon—annual heating energy demand estimated based on Yan et al. [14], upper bound of grey polygon—annual heating energy demand estimated based on Gunawan et al. [15]. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Figure 11. Heating energy produced (color lines) and projected demand (grey polygon). Lower bound of grey polygon—annual heating energy demand estimated based on Yan et al. [14], upper bound of grey polygon—annual heating energy demand estimated based on Gunawan et al. [15]. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Water 13 03526 g011
Figure 12. (a) Power and (b) heating energy produced from the waste heat (color lines) and projected demand (black line and grey polygon). Lower bound of grey polygon—annual heating energy demand estimated based on Yan et al. [14], upper bound of grey polygon—annual heating energy demand estimated based on Gunawan et al. [15]. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Figure 12. (a) Power and (b) heating energy produced from the waste heat (color lines) and projected demand (black line and grey polygon). Lower bound of grey polygon—annual heating energy demand estimated based on Yan et al. [14], upper bound of grey polygon—annual heating energy demand estimated based on Gunawan et al. [15]. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Water 13 03526 g012
Figure 13. Recovery factor over time. Grey polygon—theoretical values for the recovery factor, lower bound = 2%, upper bound = 20% (e.g., [35,78]). The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Figure 13. Recovery factor over time. Grey polygon—theoretical values for the recovery factor, lower bound = 2%, upper bound = 20% (e.g., [35,78]). The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Water 13 03526 g013
Figure 14. Power plant and other surface facilities cost for each simulated design considering Sanyal et al. [72] costs: (a) configuration A (Figure 5); (b) configuration B (Figure 5). The cost is in USD. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Figure 14. Power plant and other surface facilities cost for each simulated design considering Sanyal et al. [72] costs: (a) configuration A (Figure 5); (b) configuration B (Figure 5). The cost is in USD. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Water 13 03526 g014
Figure 15. Capital cost for each simulated design considering Sanyal et al. [72] costs: (a) configuration A (Figure 5); (b) configuration B (Figure 5). The cost is in USD. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Figure 15. Capital cost for each simulated design considering Sanyal et al. [72] costs: (a) configuration A (Figure 5); (b) configuration B (Figure 5). The cost is in USD. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Water 13 03526 g015
Figure 16. Input parameters ranked according to their influence on the capital cost. Baseline—overall simulated mean value, solid color—positive impact on the output, transparency—negative impact on the output. The cost is in USD. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Figure 16. Input parameters ranked according to their influence on the capital cost. Baseline—overall simulated mean value, solid color—positive impact on the output, transparency—negative impact on the output. The cost is in USD. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Water 13 03526 g016
Figure 17. Levelized cost of energy and probability of providing energy at lower cost than the business-as-usual (BAU) scenario: (a) configuration A and best-case scenario—heat only; (b) configuration A and base-case scenario—heat only; (c) configuration B—heat only; (d) configuration A and best-case scenario—combined heat and power. The cost is in USD. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Figure 17. Levelized cost of energy and probability of providing energy at lower cost than the business-as-usual (BAU) scenario: (a) configuration A and best-case scenario—heat only; (b) configuration A and base-case scenario—heat only; (c) configuration B—heat only; (d) configuration A and best-case scenario—combined heat and power. The cost is in USD. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Water 13 03526 g017
Figure 18. Input parameters ranked according to their influence on the levelized cost of energy. Baseline—overall simulated mean value, solid color—positive impact on the output, transparency—negative impact on the output. The cost is in USD. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Figure 18. Input parameters ranked according to their influence on the levelized cost of energy. Baseline—overall simulated mean value, solid color—positive impact on the output, transparency—negative impact on the output. The cost is in USD. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Water 13 03526 g018
Figure 19. Levelized cost of energy as a function of the uncertain parameters’ percentile for configuration A. Dashed line—heating energy cost with oil furnaces (see text for further details). The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Figure 19. Levelized cost of energy as a function of the uncertain parameters’ percentile for configuration A. Dashed line—heating energy cost with oil furnaces (see text for further details). The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Water 13 03526 g019
Figure 20. Levelized cost of energy as a function of the uncertain parameters’ percentile for configuration B. Dashed line—heating energy cost with oil furnaces (see text for further details). The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Figure 20. Levelized cost of energy as a function of the uncertain parameters’ percentile for configuration B. Dashed line—heating energy cost with oil furnaces (see text for further details). The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Water 13 03526 g020
Table 1. Sensitivity analysis for the properties of the medium.
Table 1. Sensitivity analysis for the properties of the medium.
ParameterSymbolUnitWorst-Case
Scenario
Base-Case
Scenario
Best-Case
Scenario
References
Maximum horizontal principal stressσHMPa213181152[51]
Direction maximum horizontal principal stressσH°N215°E[51]
Minimum horizontal principal stressσhMPa13810684[51]
Direction minimum horizontal principal stressσh°N305°E[51]
Vertical principal stressσVMPa12110897[51]
In situ fluid pressurePporeMPa494339[51]
Reservoir temperatureTreservoir°C3388167[28]
Permeabilityκm210−1810−1710−16[62,63]
Geological materials
Thermal conductivityλrockW m−1 K−12.42.01.5[28]
Volumetric heat capacityρcrockMJ m−3 K−12.12.42.7[28]
Young’s modulusEGPa10071.543[64]
Poisson’s ratioν---0.160.230.30[64]
Asperity strength factorFasperity---0.60.50.4---
Basic friction angleϕbasic°2624.523[65]
Initial shear dilation angleϕdilation, 0°02.55---
Peak shear dilation angleϕdilation, peak°5.01020---
Ultimate shear dilation angleϕdilation, ultimate°2.55.010---
Peak shear displacementUpeakmm2.55.010---
Residual shear displacementUresidualmm1.252.55---
Reference stress for 90% closureσn, refMPa405060---
In situ fluid
Densityρfluidkg m−3101210801112---
Dynamic viscosityωkg m−1 s−1 7.48   ×   10−4 3.19   ×   10−4 1.62   ×   10−4---
Circulation fluid
Re-injection temperatureTinjection°C303050---
Densityρfluidkg m−3993993985---
Specific heatcfluidJ kg−1 K−1418041804180---
Table 2. Sensitivity analysis for the costs.
Table 2. Sensitivity analysis for the costs.
CostSanyal et al. [72]
(Optimistic Scenario)
Factor of 2
(Likely Scenario)
Factor of 5
(Pessimistic Scenario)
Stimulation per well
(M$)
Minimum0.51.02.5
Mean0.751.53.75
Maximum1.02.05.0
Power plant and other surface facilities
($ kWh−1)
Minimum180036009000
Mean2000400010,000
Maximum2200440011,000
Annual operation and maintenance
(¢ kWh −1)
Minimum2.04.010.0
Maximum3.57.017.5
The costs are in USD.
Table 3. Levelized cost of energy and uncertainty.
Table 3. Levelized cost of energy and uncertainty.
Parameter CodeParameter DescriptionUnitVariable TypeDistribution
CwellsWells cost$DiscreteUniform
CstimulationStimulation cost$ContinuousTriang (min, mean, max)
CplantPower plant and other surface facilities cost$ kW−1ContinuousTriang (min, mean, max)
OannualAnnual operations and maintenance cost¢ kWh−1ContinuousUniform (min, max)
eannualAnnually consumed energyMWhContinuousUniform (min, max)
iImputed interest rate%ContinuousSingle value
tProject lifetimeyearContinuousSingle value
Table 4. Design parameters and results.
Table 4. Design parameters and results.
Configuration 1Fault Offset
(m)
Spacing between Wells
(m)
Open Hole Length
(m)
Vstim
(km3)
Pstim
(MPa)
Pcirc
(MPa)
Qrec
(L s−1)
Wloss
(%)
H
(MPa L−1 s−1)
Tdrawdown
(°C/year)
I*R*I*R*
Fracture network 1
A0.107006000.8252810−552.918.10.280.16
A0.017006000.825289−547.417.60.300.12
A0.0015006000.829356−449.015.50.200.10
B0.107005000.423287−522.420.40.540.60
B0.017005000.425307−528.818.20.420.73
Fracture network 2
A0.107006000.87142055.613.70.040.03
A0.017006000.8492−156.42.10.050.03
A0.0016006000.85102−150.15.30.060.02
B0.107006000.4241−122.811.80.090.50
B0.017006000.4351−126.13.40.080.35
Fracture network 3
A0.107006000.823248−450.78.00.240.07
A0.016006000.823246−446.89.90.210.05
A0.0015006000.826326−451.85.00.190.04
B0.107006000.422243−322.32.90.270.45
B0.017006000.422244−322.920.50.310.56
1. The reader is referred to Figure 4 for further details on the configuration. Vstim—stimulation volume, Pstim—stimulation pressure, Pcirc—circulation pressure, I*—injection well, R*—recovery well, Qrec—recovered flow rate, Wloss—water loss, H—hydraulic impedance, Tdrawdown—temperature drawdown. The reader is referred to Section 4.2 for further details on the fracture network.
Table 5. Design parameters and results.
Table 5. Design parameters and results.
Configuration 1Fault Offset
(m)
Spacing between Wells
(m)
Open Hole Length
(m)
Vstim
(km3)
Pstim
(MPa)
Pcirc
(MPa)
Qrec
(L s−1)
Wloss
(%)
H
(MPa L−1 s−1)
Tdrawdown
(°C/year)
I*R*I*R*
Fracture network 2
A0.107006000.8474414−293.30.40.170.47
A0.017006000.8474416−496.610.60.210.58
A0.0017006000.8484516−488.218.20.230.90
B0.107006000.843453−129.89.70.130.87
B0.017006000.843454−232.22.70.190.73
1. The reader is referred to Figure 4 for further details on the configuration. Vstim—stimulation volume, Pstim—stimulation pressure, Pcirc—circulation pressure, I*—injection well, R*—recovery well, Qrec—recovered flow rate, Wloss—water loss, H—hydraulic impedance, Tdrawdown—temperature drawdown. The reader is referred to Section 4.2 for further details on the fracture network.
Table 6. Power plant and other surface facilities cost for each simulated design.
Table 6. Power plant and other surface facilities cost for each simulated design.
DesignPower Plant and Other Surface Facilities Cost
(M$)
Factor of 2 (Likely Scenario)Factor of 5 (Pessimistic Scenario)
MinimumMeanMaximumMinimumMeanMaximum
Configuration AQ = 52.990100110225250275
Q = 47.4798897198220242
Q = 49.08696106216240264
Q = 55.697108119243270297
Q = 56.497108119243270297
Q = 50.190100110225250275
Q = 50.78696106216240264
Q = 46.8798897198220242
Q = 51.88696106216240264
Q = 93.3687684171190209
Q = 96.6728088180200220
Q = 88.2616875153170187
Configuration BQ = 22.440444899110121
Q = 28.8505662126140154
Q = 22.840444899110121
Q = 26.1455055113125138
Q = 22.338424695105116
Q = 22.940444899110121
Q = 29.8222426546066
Q = 32.2252831637077
The cost is in USD. The reader is referred to Figure 5 for further details on the configuration. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Table 7. Capital cost for each simulated design.
Table 7. Capital cost for each simulated design.
DesignCapital Cost
(M$)
Factor of 2Factor of 5
MinimumMeanMaximumMinimumMeanMaximum
Configuration AQ = 52.9104114124244269294
Q = 47.491102112215239261
Q = 49.0100110121233259283
Q = 55.6111122134260289315
Q = 56.4111122134261289315
Q = 50.1103114125244269295
Q = 50.799110121233259283
Q = 46.892102112216239260
Q = 51.899110121235259284
Q = 93.3819099189209228
Q = 96.68594104199219239
Q = 88.2748290170189208
Configuration BQ = 22.4535863116129142
Q = 28.8637077142159174
Q = 22.8535863115129142
Q = 26.1586470129144158
Q = 22.3515661112124135
Q = 22.9535863115129140
Q = 29.8353842707986
Q = 32.2384246798998
The cost is in USD. The reader is referred to Figure 5 for further details on the configuration. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Table 8. Levelized cost of energy for each simulated design.
Table 8. Levelized cost of energy for each simulated design.
DesignLevelized Cost of Energy
($ MWh−1)
HeatCHP
OptimisticLikelyPessimisticOptimisticLikelyPessimistic
Configuration ABest-case Q = 52.9136247584134244576
Q = 47.4121221517120218510
Q = 49.0132239561130236553
Q = 55.6145265625143262616
Q = 56.4145265626143262617
Q = 50.1136247584134244577
Q = 50.7131239561130235553
Q = 46.8123221517122218510
Q = 51.8132239563130235555
Base-caseQ = 93.3111195453---------
Q = 96.6114205475---------
Q = 88.2102178410---------
Configuration BQ = 22.475126280---------
Q = 28.888152343---------
Q = 22.875126279---------
Q = 26.182139312---------
Q = 22.373121268---------
Q = 22.975126278---------
Q = 29.85483170---------
Q = 32.25891192----------
CHP—combined heat and power. The reader is referred to Figure 5 for further details on the configuration. The reader is referred to Table 4 and Table 5 for further details on the flow rates.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Miranda, M.; Raymond, J.; Willis-Richards, J.; Dezayes, C. Are Engineered Geothermal Energy Systems a Viable Solution for Arctic Off-Grid Communities? A Techno-Economic Study. Water 2021, 13, 3526. https://doi.org/10.3390/w13243526

AMA Style

Miranda M, Raymond J, Willis-Richards J, Dezayes C. Are Engineered Geothermal Energy Systems a Viable Solution for Arctic Off-Grid Communities? A Techno-Economic Study. Water. 2021; 13(24):3526. https://doi.org/10.3390/w13243526

Chicago/Turabian Style

Miranda, Mafalda, Jasmin Raymond, Jonathan Willis-Richards, and Chrystel Dezayes. 2021. "Are Engineered Geothermal Energy Systems a Viable Solution for Arctic Off-Grid Communities? A Techno-Economic Study" Water 13, no. 24: 3526. https://doi.org/10.3390/w13243526

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop