Thermophysical properties of surficial rocks: a tool to characterize geothermal resources of remote northern regions

The energetic framework of Canadian remote communities relies on fossil fuels. This has adverse environmental and energy security issues. In order to offset diesel consumption, the search for local, sustainable and carbon-free energy sources is of utmost importance. Unfortunately, in such remote regions, subsurface data to evaluate the geothermal potential is often nonexistent. This raises a key question: how to characterize geothermal resources associated to petrothermal systems based on surface data? Answering this question is the purpose of this work highlighting how outcrops can be used as deep subsurface analogues. The variability induced by laboratory methods to characterize thermophysical properties is further evaluated in the estimation of the present-day temperature at depth. The community of Kuujjuaq, Canada, is used as an example where guidelines are defined to evaluate the steady-state geotherm. Rock samples were collected and analyzed with a guarded heat flow meter and an optical scanner to determine thermal conductivity. Radiogenic elements concentration was evaluated with gamma-ray and mass spectrometry. 2D temperature models were built taking into account the regional geology and the results obtained from the different laboratory methods. A base-case temperature of 57–88 °C at 5 km is predicted below Kuujjuaq. This range is based on different methods used to evaluate both thermal conductivity and internal heat generation. The work conducted in Kuujjuaq shows that the combination of gamma-ray spectrometry and optical scanning gives lower base-case temperature predictions when compared to mass spectrometry combined with the guarded heat flow meter. Despite the nonexistence of deep temperature measurements in northern regions, the assessment of thermophysical properties from outcrops is shown to be a useful tool for a preliminary assessment of geothermal resources in remote areas facing critical energy issues.

and offset the use of fossil fuels, the search for local sources of environmentally friendly energy is of fundamental interest. Amongst the renewable energy options, geothermal resources have the advantage of providing continuous heating and base-load power generation regardless of the weather conditions.
In Canada, the utilization of geothermal energy is growing annually. From 1990 to 2013, the ground-source heat pump (GSHP) market experienced a significant increase from 450 to 8250 installed units. In total, more than 110 000 GSHP units were installed throughout the southern part of the country until 2013 (Raymond et al. 2015). Assuming a linear growth of this market, more than 180 000 units might have been installed by 2019. Deep geothermal resources have been the target of recent research (e.g., Bédard et al. 2018;Ferguson and Grasby 2014;Grasby et al. 2012;Hofmann et al. 2014;Majorowicz and Grasby 2010;Minea 2012, 2015a;Majorowicz and Moore 2014;Nasr et al. 2018), but no geothermal power plant is yet producing electricity.
Geothermal investigations with a focus on the Canadian northern communities facing critical energy challenges have additionally been carried out (e.g., Comeau et al. 2017;Giordano and Raymond 2019;Grasby et al. 2013;Gunawan et al. 2020;Majorowicz and Grasby 2014;Majorowicz and Minea 2015b;Minnick et al. 2018). These studies indicate promising geothermal energy development for the off-grid communities due to the cold climate and the high energy cost. Shallow geothermal resources are seen as viable shortterm alternative solutions whereas deep geothermal development can be a long-term objective. However, the uncertainty about the depth and temperature of geothermal resources in northern Canada is considered as the main obstacle for their exploitation. Due to the lack of heat flow data in such remote regions, it is difficult to accurately assess the extent of the geothermal resources. Majorowicz and Minea (2015b) presented an evaluation of the geothermal resources for northern Québec based on sparse heat flow data and assumptions regarding the subsurface thermophysical properties. These authors mention that a temperature of 100 °C can be reached at a depth greater than 5 km, making electricity generation difficult. However, in such remote areas, where fossil fuels are the only source of energy, it is crucial to carry out detailed and local studies to avoid extrapolating sparse heat flow data over a large territory. As additionally highlighted by Eppelbaum et al. (2014), the use of average published literature values of rock thermophysical properties can indeed lead to miscalculations of the geothermal potential of a target area. Perhaps an accurate knowledge of rock thermophysical properties from surface outcrops can help better assess temperature at depth and provide the critical information to directly use deep geothermal resources for space heating and to offset fossil fuel consumption.
Deep boreholes with equilibrium temperature measurements and thermal properties assessment located near northern communities of Canada are rarely available. For example, the territory of Nunavik covering 507 000 km 2 and enclosing parts of the Superior and Churchill geological provinces only has three locations at which heat flow has been evaluated : Raglan mine, Asbestos Hill mine and Coulon (Fig. 1).
These locations are far from the communities. Raglan and Asbestos Hill mines lie at a distance of almost 500 km from Kuujjuaq and Coulon at a distance of 420 km (Fig. 1). Additionally, two other heat flow measurements located in Nunavut (Nielsen Island; Fig. 1) and Newfoundland and Labrador (Voisey Bay; Fig. 1) lie at a distance of approximately 500 km and 430 km, respectively, from Kuujjuaq. Similarly, Nunavut is facing a data gap challenge as outlined by Minnick et al. (2018), that conducted a geothermal potential assessment for this region. This raises the question: how to evaluate the depth and temperature of geothermal resources when only surface data is available? This evaluation, even if only preliminary, might trigger the interest of stakeholders for deep drilling nearby the communities. An option for overcoming the subsurface data gap problem is to use outcrops as surface analogues. Although exposed to weathering and erosion processes, outcrops can be easily accessed at low cost and provide reliable data to build geothermal conceptual models (e.g., Bauer et al. 2017;Blázquez et al. 2017a, b;Cumming 2009;Homuth and Sass 2014;Weydt et al. 2018). In the absence of any sign of geothermal activity (e.g., thermal springs) and in the presence of low-permeability crystalline rocks, conduction is the main heat transfer mechanism expected for petrothermal systems of the Canadian Shield. Therefore, the estimation of the depth and temperature of geothermal resources is mainly controlled by the heat flux, the surface temperature and the thermal conductivity and internal heat generation of the geological materials.
The work presented in this study is focused on Kuujjuaq, where geothermal evaluation has been conducted to define guidelines for other communities facing the same data gap and exploration challenges. Thermal conductivity of the field samples collected in this community was evaluated by steady-state and transient methods. The concentration of heat-producing elements was determined in the laboratory by mass and radiometric spectrometry. Additionally, hydraulic properties have been evaluated in the laboratory with transient methods to confirm the petrothermal regime assumption. The subsurface temperature distribution was then evaluated numerically by 2D steady-state heat conduction simulations. The use of the aforementioned laboratory methods was essential to assess the variability induced by the laboratory analyses on the temperature extrapolations, and, thus, define lower and upper bounds for the temperature field at depth.  (Jessop 1968), Voisey Bay (Mareschal et al. 2000), Camp Coulon (Lévy et al. 2010), Raglan and Asbestos Hill mining sites (e.g., Comeau et al. 2017)

Geology
Kuujjuaq is the administrative capital of Nunavik (Fig. 1). This village is also the largest Inuit community in Nunavik, with about 2750 inhabitants. Given its high energy demand compared to the remaining smaller communities, it is the most suitable target to carry out a geothermal energy feasibility study in northern Québec.
The main lithological units outcropping nearby the community of Kuujjuaq are paragneiss and diorite, with smaller occurrences of gabbro, tonalite and granite (Fig. 2). These lithologies are of two main origins: metamorphic/metasedimentary, comprising the paragneiss rocks; and igneous, enclosing diorite, gabbro, tonalite and granite lithologies. These rocks belong to the Southeastern Churchill Province of the Canadian Shield (e.g., Wardle et al. 2002) and are Archean to Paleoproterozoic in age (SIGÉOM 2019).
The paragneiss unit is described as a biotite-rich migmatitic paragneiss, with occurrences of millimetric garnet minerals (Lafrance et al. 2014). The diorite unit is described by Simardet al. (2013) as a very foliated to mylonitic granoblastic diorite and quartz diorite rich in hornblende and actinolite. The gabbro unit is, according to Simard et al. (2013), an amphibolite-rich granoblastic gabbro and diorite. The tonalite unit is defined as a white color tonalite and granite of mobilize type (Simard et al. 2013). The granite intrusions are described as two-mica pink-color granite to pegmatitic granite (Simard et al. 2013).
The Lac Pingiajjulik fault (Fig. 2) is described in Simard et al. (2013) and SIGÉOM (2019) as a regional thrust fault with dextral movement. The fault is characterized by a broad zone of highly recrystallized mylonite, separating different lithologies (Poirier 1989). The dextral movement was determined from pressure shadows around the porphyroblasts (Simard et al. 2013 and references therein).

Methods and techniques
A total of 24 samples were collected in the Kuujjuaq field area (Fig. 2) and core plugs with 20-mm-radius and thicknesses of 20 to 30 mm were drilled from the specimens, taking into account the observed anisotropy. Core plugs were used for laboratory experiments to infer thermal and hydraulic properties at INRS. Different laboratory methods were used to evaluate the influence of the chosen laboratory approach on the numerical modeling of the temperature distribution at depth. In all cases, the mean value of the relative error (MRE; %) between laboratory results with different methodologies was calculated as: where X stands for the parameter under evaluation and the subscript for the applied methodology.
A total of 14 samples were selected as representative of the different lithologies and analyzed for the radioisotopes 238 U, 232 Th and 40 K to estimate the radiogenic heat production. These analyses were carried out at the University of Coimbra. The samples selection took into account the concentrations of U and Th obtained through ICP-MS, their geographical position (Fig. 2) and the weathering degree.

Thermal conductivity
Thermal conductivity of dry samples was evaluated at room temperature using the guarded heat flow meter (GHFM) technique (e.g., Raymond et al. 2017;Ruuska et al. 2017; TA Instruments 2019), with a FOX50 Heat Flow Meter from TA Instruments having an accuracy of 3%. The instrument has two plates, two heat flow meters and two protective casings to prevent heat losses. The analysis is made when temperature across the sample reaches steady state. A temperature difference is imposed on both plates and successive data acquisition cycles grouped in blocks are run until the temperature of the upper and lower plates and transducer signals satisfy all the necessary equilibrium criteria to declare the sample in thermal equilibrium. Then, thermal conductivity is evaluated. Each plate must meet each equilibrium criterion independently. These criteria are (TA Instruments 2019): between a block and a previous one must change its sign or be equal to zero. Only when this final criterion is met, the equilibrium is declared, and the results are calculated.
A film of silicone paste of about 0.1 mm was smeared on both samples surfaces to improve the contact between rock sample and heating plates.
Additionally, 18 samples from the same lithological units as the core plugs were analyzed by the optical scanning technique (Popov et al. 2016 and references therein) with an infrared thermal conductivity scanner (TCS) from LGM Lippmann having an accuracy of 3%. These measurements took into account the observed anisotropy. Thermal conductivity is evaluated transiently based on solutions of the heat conduction equation for a quasi-stationary temperature field in a movable coordinate system OXYZ (Popov et al. 2016). The main elements of the TCS are a focused, mobile and continuously operated optical heat source mounted on an array of three infrared temperature sensors (Popov et al. 2016). The cold sensor passes first and records the temperature of both standards and sample before the thermal perturbation. Then, the optical heat source disturbs the temperature and, finally, the two hot sensors record the temperature after the perturbation. The sample thermal conductivity is evaluated taking into account this temperature variation and the thermal conductivity of both standards.

Radiogenic elements and heat production
The concentration of naturally occurring radioisotopes 238 U, 232 Th and 40 K was determined by gamma-ray spectrometry (GRS) using a NaI(Tl) detector (7.62 × 7.62 cm) connected to a multichannel pulse-height analyzer (1024 channels) equipped with a spectrum stabilizer for automatic compensation of gain shift. The detector is surrounded by a 5-cm-thick lead shield to smoothen background gamma-radiation (ORTEC 2015). The isotope concentrations are measured using the three-window method. This involves the detection of the gamma-radiation emission on the decay of 214 Bi ( 238 U decay chain), 208 Tl ( 232 Th series) and 40 K. The analyzed portion of the gamma-ray spectrum ranges in energy from 0 to 3000 keV. 40 K as an energy peak of 1460 keV. 214 Bi has an energy peak of 1764 keV. 208 Tl has the most energetic peak with 2614 keV (e.g., Lamas et al. 2017). The system is calibrated with standard solutions certified by the International Atomic Energy Agency (IAEA) for K, U and Th activity measurements. The potassium calibration standard is extra-pure potassium sulphate (99.8%) and uranium and thorium content lower than 0.001 mg kg −1 and 0.01 mg kg −1 , respectively. The uranium standard is U-ore diluted with silica, containing a negligible amount of K (< 0.00234 mg kg −1 ) and Th (< 1 mg kg −1 ). The thorium standard is Th-ore diluted with silica, with trace content of uranium and potassium. The background gamma-radiation subtraction is performed for each measurement. The counting time for each sample was increased for 24 h due to the results from the trace elements concentration (Table 1). The elemental concentration on U, Th and K is then calculated based on the daughter isotopes activity.
Additionally, GRS results were compared with those obtained from ICP-OES/MS (Table 1), as mass spectrometry methods are referred to be more accurate by six orders of magnitude than radiometric methods (Hou and Roos 2008). In ICP, the chemical elements contained in the sample solution are decomposed into their atomic constituents in an inductively coupled argon plasma. Then, the positively charged ions are extracted from the inductively coupled plasma into a high vacuum via an interface. These ions are then separated by mass filters and finally measured by an ion detector (e.g., Hou and Roos 2008). Quality control procedures were considered during the analyses to guarantee the accuracy of laboratory measurements. First, the ICP was calibrated using a blank and the working calibration standard. Then, the initial calibration verification standard was run, and the percent of recovery was ± 10%. Following the initial calibration verification, the initial calibration blank was analyzed. The concentration was verified to be less than the reporting limit for each element. The reporting limit standard was run, followed by the spectral interference check solution, the continuing calibration verification (CCV) and the continuing calibration blank (CCB). Afterwards, the method blank, the laboratory control samples, and the samples themselves were analyzed. The CCV/CCB was run every 10 samples. At the end of the analytical sequence, a final CCV/CCB was analyzed. The following certified reference materials (CRM) were used: • WPR-1a-CRM for a peridotite with rare earth and platinum group elements; • QLO-1-CRM for quartz latite; • SGR-1-CRM for green river shale; • SY-4-CRM for diorite gneiss.
Radiogenic heat production (A; W m −3 ) was calculated afterwards knowing that: where ρ (kg m −3 ) is the density, P (wt %), H 0 (W kg −1 ) and c (mg kg −1 ; %) are the specific abundance, heat generated and concentration of each radioisotope, respectively. The constants P and H 0 for each element were estimated with Rybach (1976) approach, transforming Eq. (2) in the following empirical function (Rybach 1988): where potassium concentration was converted from its oxide form to the elemental form by:

Hydraulic properties
Porosity was evaluated using the combined gas permeameter-porosimeter AP-608 from Core Test following Boyle's law (e.g., Coretest Systems, Inc. 2019; Raymond et al. 2017). This law states that the pressure exerted by a given mass of an ideal gas is inversely proportional to the volume it occupies  and references therein).
Permeability was also evaluated using the AP-608. The analysis follows the transient pressure decay method and the permeability is inferred by Darcy's law. Klinkenberg correction is then applied to convert the gas to liquid permeability  and references therein). Density of solid grains was evaluated with the AP-608 grain volume chamber. (

Temperature field at depth
Preliminary estimation of the present-day temperature at depth was solved numerically using the COMSOL Multiphysics software assuming 2D steady-state heat conduction: where λ (W m −1 K −1 ) is the thermal conductivity, T (°C) is the temperature and x (m) and z (m) are the spatial variables.
Geological models of the lithosphere stratigraphy and thickness in Kuujjuaq were built based on available regional literature data (e.g., Bourlon et al. 2002;Christensen and Fountain 1975;Hall et al. 1995;Mareschal et al. 1990;Moorhead et al. 1989;Poirier 1989;Seguin and Goulet 1990;St-Onge et al. 2002;Telmat 1998;Telmat et al. 1999;Vervaet 2016;Wardle and van Kranendonk 1996). Moho is placed at 37.5 km depth. The upper crust is mainly composed by paragneiss and has a thickness of 26.5 km. The lower crust is assumed to be made of granulitic facies rocks with a thickness of about 11 km. Its radiogenic heat production is estimated to be 0.45 × 10 −6 W m −3 (Ashwal et al. 1987). The Moho contribution to the heat flow is estimated to be, on average, 15 × 10 −3 W m −2 (e.g., Mareschal 2007, 2014;Jaupart et al. 2016;Mareschal and Jaupart 2013).
The temperature distribution was modeled for a rectangular geometry with a width of 8 km and a depth of 10 km (Fig. 3). The center of the 2D model corresponds to a change in lithology between paragneiss and diorite (Fig. 2). The diorite layer thickness is considered to vary within 1 to 5 km. This geometry was selected to evaluate how lithological changes with different thicknesses of the diorite can influence the temperature at depth. A constant surface temperature of − 1 °C (climate normals 1981-2010; Comeau et al. 2017; Environment and Climate Change Canada 2019) was set as the upper boundary condition. The lower boundary condition is the heat flux at 10 km depth. This parameter is calculated as: where Q (W m −2 ) is the heat flux, z (m) is depth and h (m) is thickness. The subscript 10 stands for 10-km depth and M for Moho. LC and UC are, respectively, lower and upper crustal layers. Both lateral boundaries are assumed adiabatic. Thermal conductivity and radiogenic heat production were defined according to lithologies, assuming uniform values. Several scenarios were studied with the goal of assessing the influence of the different laboratory methods in the prediction of the temperature field at depth. These are: (1) scenario GHFM-GRS, (2) scenario TCS-GRS, (3) scenario GHFM-ICP-MS, and (4) scenario TCS-ICP-MS. The thermophysical properties of each geological material were, thus, varied according to the applied method. It is important to highlight that only the samples evaluated by the four methods were considered for the numerical simulations to avoid a misinterpretation of the results. The worst and best temperature predictions were calculated within each scenario considering the statistical distribution of thermal properties.
Long-lasting changes in surface temperature, which can occur during glacial episodes, have an effect on the subsurface temperature distribution due to downwards thermal diffusion (e.g., Beardsmore and Cull 2001;Bédard et al. 2018 and references therein; Beltrami et al. 2005;Birch 1948;Chouinard and Mareschal 2009;Chouinard et al. 2007;Jessop 1990;Mareschal et al. 1999;Rath et al. 2012;Suman and White 2017). The temperature predictions were therefore corrected to evaluate by how much the temperature at 5 km depth can have been misestimated by assuming a constant Dirichlet condition in the numerical simulations. The temperature correction was based on Carslaw and Jaeger's (1959) solution: where ΔT (°C) is the departure from original equilibrium temperature at depth z (m) and time t (s) after an instantaneous change in surface temperature T 0 (°C), α (m 2 s −1 ) is the thermal diffusivity of the geological materials (assumed as 1 × 10 −6 m 2 s −1 ) and erfc(x) is the complimentary error function.
The duration and surface temperature perturbations of each Pleistocene climate events in Northern Canada is still debatable and several models have been proposed (e.g., Beltrami et al. 2003;Birch 1948 and references therein; Chouinard and Mareschal 2009;Jaume-Santero et al. 2016;Majorowicz et al. 2005;Mareschal et al. 1999;Pickler et al. 2016). Therefore, three different ground surface temperature histories (GSTH) were used based on Birch (1948) and Jessop (1990) model of the Pleistocene glaciations (Fig. 4). For the first scenario, the temperature during a glacial cycle was considered as 10 °C colder than today. The second assumes a temperature of 5 °C colder and, for the third, a temperature of 1 °C colder than current times was considered (Fig. 4).

Rock samples description and geochemistry
The paragneiss samples (P1-P8; Fig. 2) collected present high content in amphibole and biotite minerals, in a fine-to medium-grained matrix. The feldspar minerals are weathered showing a brown to light-brown tone. The main mineral phases identified in thin sections are, on average, 16% quartz, 29% feldspar and 55% mafic minerals. The diorite samples (D1-D5; Fig. 2) are a fine-to medium-grained matrix and do not have evidence of weathering. They have an averaged content of 17% quartz, 23% feldspar and 59% mafic minerals. The gabbro samples (G1-G4; Fig. 2) have a fine-to medium-grained texture and their feldspar minerals have signs of weathering given by their brownish tone. The analyzed samples from this unit have 25% quartz, 24% feldspar and 51% mafic minerals content. The tonalite samples (T1-T4; Fig. 2) collected are coarse-grained, and their color varies from white to pinkish. The samples analyzed are characterized by a higher content of quartz (41%) and feldspar (51%) than of mafic minerals (8%). The granite samples (Gr1-Gr3; Fig. 2) are composed, on average, of 35% quartz, 55% feldspar and 10% mafic minerals. One of the granite samples collected presents foliation. Based on the similar texture and mineralogical content, the rock samples collected can be classified in three main groups: paragneiss, diorite-gabbro and tonalite-granite. The content of major and trace elements of rock samples were analyzed by inductively coupled plasma (ICP) techniques, more precisely, optical emission spectrometry (OES) and mass spectrometry (MS) at INRS (

Thermal properties
The paragneiss samples are characterized by an average thermal conductivity of 2.32 W m −1 K −1 when evaluated with the steady-state method (GHFM) and 2.52 W m −1 K −1 with the transient method (TCS; Table 2; Fig. 5). For the igneous mafic group (diorite-gabbro), both steady-state and transient methods give a similar value of Fig. 4 Timeline of Pleistocene and Holocene climate events and temperature step. Solid line-temperature during glacial episode assumed as 10 °C colder than today; pointed line-temperature during glacial episode assume 5 °C colder than today; dashed line-temperature during glacial episode 1 °C colder than today 2.83 W m −1 K −1 and 2.84 W m −1 K −1 , respectively. Thermal conductivity of igneous felsic samples evaluated by the steady-state method is, on average, 3.08 W m −1 K −1 and 3.36 W m −1 K −1 when evaluated by the transient method. The lower MRE is obtained for the igneous mafic group (diorite-gabbro), with a value of − 2%, considering the steadystate value at the denominator. Regarding the paragneiss samples and the felsic igneous rocks, the MRE is − 15% and − 13%, respectively. On average, the relative error between steady-state and transient methods for all the lithologies is − 9.8%, ranging from a maximum of 30% to a minimum of − 58%.

Radiogenic elements and heat production
Results obtained for radiogenic elements from radiometric (GRS) and mass (ICP-MS) spectrometry methods allowed to estimate the internal heat generation with Eq. (2) for the main lithologies in Kuujjuaq (Table 3; Fig. 6). The results of GRS reveal that, on average, paragneiss samples are characterized by a uranium concentration of 1.50 mg kg −1 , a thorium concentration of 3.75 mg kg −1 and a potassium concentration of 0.50%. Lower values are found for diorite-gabbro group with average concentration in uranium, thorium and potassium of 0.38 mg kg −1 , 1.94 mg kg −1 and 0.47%, respectively. The igneous felsic group is characterized by higher uranium and potassium concentration of 1.3 mg kg −1 and 0.88%, respectively, than the igneous mafic, but a similar concentration in thorium of 1.82 mg kg −1 . The evaluation of the heat-producing elements by ICP-MS leads to a decrease on uranium concentration of about 60% and 20% for the paragneiss and tonalite-granite groups, respectively. The igneous mafic group reveals an increase of 50% for uranium. Thorium concentration decreased by 11% for the paragneiss samples when evaluated by ICP-MS, whereas for both igneous groups, this element increased up to 70%. Potassium concentration, in turn, increased by more than 60% when evaluated by ICP-MS for all the lithological groups under study.

Table 3 Heat-producing elements, density and radiogenic heat production statistics of the main lithologies in Kuujjuaq
A radiogenic heat production, ρ density  The different concentrations obtained by the two methods have an influence on the radiogenic heat production (Table 3; Fig. 6). On average, the paragneiss samples generate about 9% less internal heat when the ICP-MS method is used to evaluate the radiogenic elements. In turn, both igneous groups show an increase of more than 50% of internal heat generation when their radioactive isotopes are analyzed by ICP-MS.

Hydraulic properties
Paragneiss samples are characterized by an averaged porosity value of about 6%, whereas the igneous specimens have an averaged porosity of 4% for the diorite-gabbro group and 3% for the tonalite-granite. The matrix permeability for all the samples analyzed is below 10 −19 m 2 , which is beyond the detection limit of the instrument used.

Temperature field at depth
Thermal conductivity and internal heat generation were varied according to the analysis method (Table 4) to study the variability induced by laboratory methods on the steadystate temperature extrapolation. Moreover, due to uncertainty on the diorite-gabbro thickness, this layer was considered to vary from a minimum of 1 km to a maximum of 5 km (Table 4).
Results reveal that a high temperature is found when combining the minimum value obtained for thermal conductivity with the maximum of internal heat generation. This corresponds to the best-case scenario of Table 6. The worst-case scenario is obtained by using the maximum thermal conductivity with the minimum heat production. The median values of the thermophysical properties were used to calculate the base-case scenarios. This choice took into account the population standard deviation calculated for each thermophysical parameter (Tables 2, 3). The median was considered to be statistically more robust than the population average value.
Using GHFM to evaluate thermal conductivity and GRS for the internal heat generation, the temperature at 5 km can vary from a minimum of 31-39 °C to a maximum of 129-168 °C. The base-case scenario of this simulation points towards temperatures of 66-85 °C (Table 6). These temperature values change when TCS is used instead of GHFM. For the scenario TCS-GRS, at 5 km depth, temperature is predicted to vary from  (Table 6). Another interesting aspect is the influence of the diorite-gabbro layer. A thin diorite-gabbro layer of 1 km has a negligible influence on the temperature at 5 km. However, for a 5-km-thick layer, the isotherms show a slight decrease in this layer when compared to the adjacent paragneiss (Fig. 7). This is due to the contrast in the thermal conductivity between the paragneiss and diorite-gabbro. The only exception occurs for the base-case scenarios with TCS-GRS and TCS-ICP-MS, since both median thermal conductivity values are similar (Table 4). However, the simulations carried out indicate that the thickness of the diorite-gabbro layer induces a variability of about 4% only on the temperature field at depth. This variability is 12% to 14% when the two different techniques to evaluate the internal heat generation and the thermal conductivity are considered. A higher variability, of more than 50%, is obtained as a result of    the lithological intrinsic heterogeneity associated to the statistical distribution of rock thermal properties.
Considering the base-case scenarios, the laboratory methods can be ranked in terms of temperature field at depth such that (Table 6;   However, considering the variability induced by the intrinsic heterogeneous character of each lithological unit, the rank is: Pleistocene climate events appear to have disturbed the temperature in the shallower part of the crust by up to 4.5 °C when considering Eq. (7) to correct the simulated temperature profiles (Fig. 8). At 5 km depth, the influence of the thermal perturbation induced by the GSTH becomes minimal with correction factor up to 1 °C. However, this correction is dependent on the assumed temperature steps between climate events (Fig. 8). For example, considering the glacial episodes with a temperature of 10 °C lower than nowadays, the necessary correction is 1 °C. However, a temperature difference of 1 °C for the glacial episodes leads to a correction factor of 0.03 °C (Fig. 8).

Thermal properties
The thermal conductivity of the different lithologies (Table 2) is within the values presented in the database of Eppelbaum et al. (2014, and references therein) and Schön (2011, and references therein). Moreover, the thermal conductivity results are in TCS-ICP-MS < GHFM-ICP-MS < TCS-GRS < GHFM-GRS. Fig. 8 Subsurface temperature perturbation caused by Pleistocene climate events. Solid line-temperature during glacial episode assumed as 10 °C colder than today; pointed line-temperature during glacial episode assume 5 °C colder than today; dashed line-temperature during glacial episode 1 °C colder than today accordance with the mineralogical composition of the rock samples. Specimens with a higher percentage of quartz minerals, i.e., tonalite-granite group, reveal higher thermal conductivity (3.1 to 3.3 W m −1 K −1 ; Table 2) than the samples with lower content in quartz (paragneiss and igneous mafic group; Table 2).
Two methods were applied to evaluate the thermal conductivity: steady state (GHFM) and transient (TCS). The average MRE between the two methods is − 9.8% for all lithologies (Table 2; Fig. 5).  compared the same methods using a different dataset and obtained a mean absolute value of the relative error of 9.8%, ranging from a minimum value of 1.7% to a maximum of 23.1%. These authors proposed that the differences can be caused by sample preparation and inherent heterogeneity associated with the rock itself, rather than the intrinsic accuracy of the device. It should also be noticed that the transient method is to evaluate the thermal conductivity on a larger sample length but on a smaller volume than the steady-state method. For the latter, a volume of rock sample of 2.4 to 1.6 × 10 −3 cm 3 was used. Taking into account the results of this work, both steady-state and transient methods can be compared together to avoid a biased evaluation of the thermal conductivity.

Radiogenic elements and heat production
Studies found in the literature mentioned that both gamma-ray and mass spectrometry should give similar results when evaluating the radiogenic elements concentration (Chiozzi et al. 2003;Zhu et al. 2017). However, we obtained MRE between − 61 and 49% for U, − 11 to 72% for Th and more than 60% for K. The standards for each method were analyzed before and during the analyses to assure the calibration of both devices. These differences may be related to the analytical error of each method and to the low concentration of the radioisotopes common to Precambrian rocks. Few samples fell within the minimum detectable activity of gamma-ray spectrometry. The analytical error for gamma-ray spectrometry for the analyzed dataset is, on average, 6% for K, 35% for Bi and 21% for Tl. Regarding ICP-MS, the analytical error varies within 2-3%, increasing up to 30% if closed to the detection limit. However, the most appropriate method cannot be determined from the sole results of this work. A larger dataset including equilibrium temperature profile is needed to do so, but this falls outside of the scope of the present work. Moreover, even if ICP-MS technique has been reported as more accurate by six orders of magnitude than radiometric methods (Hou and Roos 2008), the results obtained by the NaI(Tl) system (GRS) are better correlated with the range of values found for similar Canadian Shield rocks (e.g., Perry et al. 2006;Phaneuf and Mareschal 2014;Rolandone et al. 2002). Additionally, the radiogenic heat production evaluated by the two laboratory methods is in accordance with the rock samples' geochemistry. An increase of both U and Th as a function of the SiO 2 content (Table 1) is observed. Based on these results, it can be highlighted that both methods must be compared together in order to prevent a misestimation of the internal heat generation, mainly in such old Precambrian rocks common to the Canadian northern regions.

Hydraulic properties
Porosity and permeability of the analyzed samples show low values, less than 6% and 10 −19 m 2 , respectively. These values justify the use of dry samples for the thermal conductivity analysis. A matrix permeability lower than 10 −19 m 2 and a thermal conductivity in the range of 2.3 to 3.3 W m −1 K −1 (Table 2) constrain the heat transfer mechanisms naturally occurring in the subsurface to conduction-dominated. Taking into account the thermofacies concept proposed by Sass and Götz (2012), the characterization of thermophysical properties confirms a petrothermal regime. According to Moeck (2014) classification, the study area is a conduction-dominated geothermal play of basement type. Therefore, technologies such as Engineered Geothermal Systems (EGS) and deep borehole heat exchangers (DBHE) are more adequate to exploit the geothermal resources. Permeability needs to be increased to induce advection and operate an EGS, while DBHE can be operated taking advantage of heat conduction only.

Temperature field at depth
The lack of subsurface geological knowledge (e.g., geophysical data, well logs) and detailed geochemical analyses of the mineral phases makes it difficult to accurately constrain the thickness of the diorite-gabbro layer. However, the temperature simulations reveal that the influence of this layer is minimal (about 4%) when compared to the variability originating from the laboratory methods (12-14%) and the heterogeneity of each lithological unit (more than 50%).
The heat flux inferred at surface (Table 5) is within the range of typical values found in the Superior, Nain and Churchill geological provinces north of the 55 th parallel (Fig. 1). Nielsen Island and Camp Coulon are in the Superior geological province, while Kuujjuaq, Raglan and Asbestos Hill are in the Churchill province. Voisey Bay belongs to the Nain geological province. These few heat flux assessments suggest a value that is about 10 mW m −2 higher in the Churchill province than in Superior and Nain. This can be related with the age, structure and composition of the geological provinces (e.g., Jaupart and Mareschal 2007). Within Churchill province, Raglan and Asbestos Hill are located in the Ungava orogen while Kuujjuaq is in the limit of Labrador Trough, where heat flux has never been assessed through borehole measurements. Based on the obtained results, both orogens appear characterized by similar heat flux values.
The evaluation of both thermal conductivity and radiogenic heat production is influenced by the laboratory method (Tables 2, 3), which, consequently, induces variability on the predictions of the temperature field at depth (Table 6). Without deep boreholes to evaluate the prevailing temperature at depth, it is critical to consider this variability. For the Kuujjuaq dataset, TCS and GRS lead to lower temperature values at 5 km than GHFM and ICP-MS. In the former scenario, the base-case temperature ranges from 57 to 71 °C, whereas for the latter it varies between 68 and 88 °C (Table 6; Fig. 7). This corresponds to a difference of about 18%. Nevertheless, the highest variability on the temperature prediction is undoubtedly caused by the intrinsic heterogeneous character of each lithological unit. Scenario GHFM -GRS presents the maximal value (77%) whereas TCS-ICP has the minimal variability (55%).
Surface temperature variations indeed disturbed the subsurface temperature at shallow depths, but their effect diminishes by about 80% at 5 km depth, becoming almost negligible (Fig. 8). At such depths, the paleoclimate effect is minimal when compared with the uncertainty related to the analytical methods, or to the intrinsic heterogeneous character of each lithological unit. In the absence of borehole temperature data, the assumption of a constant surface temperature seems justifiable for a preliminary characterization of heat flux and subsurface temperature.
Considering the worst-case scenarios, the exploration of deep resources is nonviable (Fig. 9a). Temperatures of 50 °C are predicted to be found at depths greater than 5 km. On the other hand, the best-case scenarios provide temperature estimates that are suitable for electricity generation representing an option for further exploration (Fig. 9c). Assuming that the base-case scenarios are the most likely to occur, then direct heat production is the most logical option to develop geothermal energy in the community of Kuujjuaq (Fig. 9b).  Base-case scenario z (km) 3.5-6.5 4-7 4-8 4-8 3-6 3.5-7 4-7.5 4-8 An average drilling depth of 4 to 7 km is estimated to reach the temperature interval 50-100 °C (Table 7). This can have an impact on the potential geothermal project's rate since drilling cost increases nonlinearly with depth (e.g., Augustine et al. 2006). For example, in Nunavut territory, northern Canada, the drilling of a full-size production well can cost approximately $12 million USD for 4 km depth and up to $30 million USD if the well has a length of 8 km (e.g., Minnick et al. 2018).
This new evaluation of the temperature field at depth based on local rock sampling and outcrop analogues points towards potential development of deep geothermal resources for direct use. Moreover, this work demonstrates that the study of outcrop analogues is essential to build geothermal conceptual models, especially in remote areas away from locations with heat flow assessment. However, there is still an epistemic uncertainty in the temperature predictions that can only be overcome with deep temperature measurements. This is a key to advance to the stage of geothermal resource exploration.

Conclusions
The energetic framework of remote and off-grid communities in northern Canada, heavily relying on fossil fuels, needs to change. Geothermal energy is considered a potential solution, but a profound data gap exists. While it can be difficult to accurately define geothermal resources in such remote regions, the critical energy situation raises a key question that was addressed in Kuujjuaq. How can the depth and temperature of geothermal resources be evaluated when only surface data is available?
The guidelines drawn from the work conducted in Kuujjuaq are an example to other remote communities of Nunavik, Nunavut and Nunatsiavut facing the same energy development challenges. The estimation of the steady-state temperature field at depth requires an evaluation of the thermal conductivity and internal heat generation of the geological materials. In the lack of subsurface data, rock samples are collected from outcrop analogues. Often, the evaluation of the thermal conductivity is carried out either by steady-state or transient methods. Similarly, radiogenic heat production is usually estimated based on gamma-ray or mass spectrometry. The variability induced to the prediction of the temperature field at depth and originating from the laboratory methods was assessed in this study. The use of the four aforementioned methods indicates that thermal conductivity and radiogenic heat production are affected by laboratory analysis. The results revealed that thermal conductivity evaluated by TCS is 8% higher than by GHFM. Similarly, radiogenic heat production estimated based on ICP-MS results is, on average, 36% higher than based on GRS. This, consequently, influences the temperature predictions. In this case, the TCS combined with GRS gives lower base-case temperature predictions at depth, while upper base-case temperature predictions are found with the GHFM and ICP-MS. However, the variability induced by the laboratory methods is smaller (less than 15%) than the one resulting from the intrinsic heterogeneity of each lithological unit (more than 50%). The simulation TCS-ICP-MS reveals lower variability (55%) between the worst-and best-case scenarios when compared to GHFM-GRS (77%).
Regional geophysical data and local geological mapping are useful to build geological conceptual models, essential to carry out the temperature simulations at depth. The absence of local subsurface geological model is another source of uncertainty such that the thickness of the diorite-gabbro layer cannot be accurately constrained. Nevertheless, more detailed geochemistry analyses on the mineral phases can overcome this data gap. Geobarometry allows to infer the pressure at which a mineral or mineral assemblage formed (e.g., Mukherjee 2011;Wendlandt 1999). If reasonable assumptions are made to convert pressure to depth, then, this methodology together with the erosion rate can help narrowing the range of possibilities for the present-day thickness of the dioritegabbro layer (e.g., Burbank 2002). Such analyses were out of the scope of the present work but will be envisioned for further studies. However, as the simulations reveal, its influence on the temperature field is minimal (about 4%). In Kuujjuaq, the preliminary steady-state simulations reveal a base-case temperature varying between 57 and 88 °C at 5 km depth, pointing to a minimum drilling depth of 3-4 km to reach the necessary temperature to use the geothermal resources for space heating.
In conclusion, this work demonstrates that, even with the arising uncertainty due to the lack of subsurface temperature data and geophysical studies, it is possible to characterize the temperature of deep geothermal resources in remote regions. The study of outcrop analogues is essential to provide a reliable evaluation of the temperature conditions that can prevail in petrothermal systems hosted in ancient, here Precambrian, rocks. Further numerical modeling will address if deep enhanced geothermal systems are technically suitable to replace the fossil fuels consumption and provide higher energy security for the communities north of the 55th parallel.