A numerical approach to infer terrestrial heat flux from shallow temperature profiles in remote northern regions

Highlights

• Heat flux evaluated numerically with a 1D inverse heat conduction problem.

• Time-varying upper boundary condition reproducing paleoclimate events.

• Temperature and pressure-dependent thermal conductivity considered in the model.

• Heat flux in the range 31.8 – 69.4 mW m⁻² depending on paleoclimate and thermophysical properties conditions.

• Low-cost method to help advancing geothermal exploration in remote northern regions.

Abstract

Assessing the geothermal potential in remote northern regions is challenging. First, due to remoteness and, second, due to data gaps. Often, Earth heat flux is inferred from deep equilibrium temperature profiles. However, in remote areas where deep boreholes are nonexistent, an adapted methodology is needed. Therefore, this study presents a new versatile and reliable numerical approach that simulates climate events and infers heat flux from a shallow temperature profile (80 m depth). Thus, addressing the challenge linked with heat flux assessment in northern regions and offering an original contribution to advance geothermal research in remote areas facing energy development challenges.

Introduction

Heat flux, together with the thermophysical properties of the geological materials and the geothermal gradient, plays a major role in the evaluation of the subsurface temperature distribution and, hence, in assessing the geothermal potential of a target area (e.g., Bédard et al., 2018; Gascuel et al., 2020). Two main methods have been developed to infer heat flux from temperature-depth (T-z) profiles: the interval/product method and the Bullard method (e.g., Beardsmore and Cull, 2001; Jessop, 1990; Powell et al., 1988). However, the application of these methods implies T-z profiles that are deep enough to observe the linear increase in temperature with depth and that are corrected for both artificial and natural disturbances (e.g., drilling, true-vertical depth, free thermal convection, groundwater flow, topography, climate change, etc.; Beardsmore and Cull, 2001; Jessop, 1990; Powell et al., 1988). Neglecting the corrections for these features, mainly topography and climate, may lead to important misestimations of heat flux (Beardsmore and Cull, 2001; Jessop, 1990; Powell et al., 1988). For example, ignoring the effects of paleoclimate has been shown to underestimate the heat flux by 10 % or more (e.g., Beck, 1977; Beck and Shen, 1989; Bédard et al., 2018; Birch, 1948; Crain, 1968; Jessop, 1971; Majorowicz and Wybraniec, 2011; Mather et al., 2018; Suman and White, 2017; Westway and Younger, 2013). This is due to the downward diffusion of the surface temperature variations with larger amplitudes and duration than daily or annual cycles (e.g., Beardsmore and Cull, 2001; Bodri and Cermak, 2007; Jessop, 1990). Solutions to quantitatively correct paleoclimate effects and to reconstruct the ground surface temperature history (GSTH) have received attention since the first observations of the impact of climate events on T-z profiles (e.g., Birch, 1948; Clauser, 1984; Crain, 1969; Harris and Chapman, 1998; Huang et al., 1996; Mareschal and Beltrami, 1992; Nielsen and Beck, 1989; Shen and Beck, 1992). However, the application of such methods to invert both the heat flux and the GSTH requires a temperature profile that is deep, at least 300 m, to reach a depth were paleoclimate perturbations are less important (e.g., Beltrami et al., 2011). In remote territories, where only shallow boreholes are available (e.g., boreholes drilled for groundwater monitoring), a different methodology is needed to provide a first reliable estimate of heat flux and justify further work to drill deep geothermal exploration boreholes. Therefore, a numerical approach was developed in this work to correct a shallow temperature profile for past and recent surface temperature variations and, thus, infer the present-day heat flux. The approach is based on the work of Márquez et al. (2019) and has been improved to numerically consider both quaternary glaciations and Holocene climate events as an upper boundary condition to the model. The method was used to infer heat flux in the community of Kuujjuaq (Nunavik, Canada) where a temperature profile was measured in a groundwater monitoring borehole having 80 m depth that has been drilled prior to this study. The present numerical approach does not address topography effect since changes in surface elevation in the study area are negligible, neither groundwater advection effect since no groundwater flow perturbation was detected in the temperature profile. Moreover, it does not deal with sedimentation and/or erosion through time that affect the surface and induce terrain effect. Furthermore, the present study does not intend to reconstruct the GSTH since the temperature profile measured is unfortunately too shallow to carry out such reconstructions. The goal of this work is to provide a tool for a first-order assessment of the present-day terrestrial heat flux by numerically correct a shallow temperature profile for literature-based surface temperature variations considering heat conduction only. The idea is to provide new geothermal exploration methods using means that are affordable to remote communities and adapted to northern regions facing data gap challenges, trying to advance geothermal exploration to the benefit of the indigenous population.

The heat flux is solved numerically with the finite element method (FEM) as a 1D inverse heat conduction problem (IHCP), considering known transient surface temperature and thermophysical properties. A shallow T-z profile acquired during a field campaign is used to find the control variable, i.e. the heat flux at the base of the model, that minimizes the objective function. In other words, the goal is to find the heat flux for simulated temperature to best match measured temperature. An arising problem of carrying out climate corrections is linked with the imperfect knowledge about the duration and temperature of each Pleistocene event. Emiliani (1955) and Flint (1947) proposed a fourfold chronology for the Quaternary glaciations that has been used by Jessop (1971, 1990) to correct heat flow measurements in Canada. The same chronology has been used, for instance, by Bédard et al. (2018) and Márquez et al. (2019) to correct T-z profiles in southern Québec. Although more recent literature proposes a different timeframe with the beginning and end of the glaciations at an older time (e.g., Ehlers and Gibbard, 2011; Jennings et al., 2013), the same fourfold stratigraphic framework was used in this work for the numerical simulations. The basal temperature during a glaciation episode can vary regionally. For instance, the surface temperature during the last glacial maximum (LGM) has been evaluated 10 °C colder than present in eastern Canada, while in central Canada studies point 5 °C, or less, colder than current times (e.g., Chouinard and Mareschal, 2009; Majorowicz et al., 2012; Majorowicz and Safanda, 2015; Mareschal et al., 1999; Pickler et al., 2016; Rolandone et al., 2003; Sass et al., 1971; Sugden, 1977). Hence, a sensitivity analysis was undertaken to assess the impact of the different temperature signals during glacial periods on the heat flux estimate. Moreover, an extensive literature review was carried out to define the timeframe and temperature of the major Holocene events that can influence the subsurface temperature. Additionally, the effect of temperature, pressure and water-saturation on thermal conductivity was considered and implemented in the numerical model. These thermal conductivity conditions have been shown to influence the heat flux (e.g., Förster et al., 2020; Harlé et al., 2019; Lerche, 1991; Nasr et al., 2018; Norden et al., 2020). Furthermore, the influence of the statistical distribution of the thermophysical properties on heat flux cannot be neglected (Miranda et al., 2020) and this uncertainty was considered in the numerical model as well.

The methodology described in this work was used to estimate heat flux in the context of geothermal energy sources assessment. In northern Canada, 239 remote communities spread over a territory of more than 3 500 000 km² rely exclusively on diesel fuel for electricity and oil for space heating (Arriaga et al., 2017; Grasby et al., 2013). Fuel is transported by boat, plane or trucks and sometimes stored in old facilities increasing risks of spill (Arriaga et al., 2017). This unsustainable energetic framework needs to change and local and clean sources of energy, such as geothermal, might be a solution (Giordano and Raymond, 2019; Grasby et al., 2012, 2013; Gunawan et al., 2020; Kanzari, 2019; Kinney et al., 2019; Mahbaz et al., 2020; Majorowicz and Grasby, 2010a, b, 2014, 2020; Majorowicz and Minea, 2015; Minnick et al., 2018). However, a large data gap exists in such remote areas. For example, in Nunavik, a 507 000 km² territory, only 3 locations have boreholes deep enough to evaluate heat flux in a conventional manner: Raglan mine, Asbestos Hill mine and Coulon exploration camp, and these lie at a distance of approximately 420 km (camp Coulon) to 500 km (Raglan and Asbestos Hill mining sites) from Kuujjuaq (Fig. 1). Two other deep boreholes, one in Nunavut (Nielsen Island; Fig. 1) and another in Newfoundland and Labrador (Voisey Bay; Fig. 1), are also located around 500 km and 430 km, respectively, from Kuujjuaq.

This data gap highlights the need to adapt methodologies for the data sources that are available in a small radius (< 4 km) from the northern remote communities. The numerical approach described was applied to Kuujjuaq (Nunavik, Canada) as a case study, but is versatile enough to be utilized in other remote northern communities facing the same data gap and geothermal exploration challenges.