Geothermal resource assessment of remote sedimentary basins with sparse data: lessons learned from Anticosti Island, Canada

Anticosti Island is located in the Anticosti sedimentary basin, an Ordovician/Silurian carbonate platform. This platform is mainly composed of limestone and shale with some dolomite and sandstone and reaches up to 5 km depth in the southwest. It overlies a Precambrian basement of the Grenville Province made of magmatic and metamorphic rocks. Like most remote and off-grid regions in Canada, it relies heavily on fossil fuels for energy supplies. An assessment of deep geothermal resources was achieved in this area with the objective of diversifying energy resources to help develop renewable energy for villages deserved by micro-grid systems. Despite sparse and low-quality bottom-hole data (15 wells of 1111 m to 2762 m depth), a 3D temperature model was developed for this sedimentary basin and its underlying Precambrian basement up to 40 km (mantle depth). Quantifying confidence intervals for thermal parameters, namely bottom-hole temperature, thermal conductivity, heat generation rate and mantle heat flux, was paramount to obtain a reliable range of temperature predictions. A high variability of modeled temperature, up to 41% at the base of the sedimentary basin and 70% at mantle depth, remains when trying to constrain input parameters. The lack of equilibrium temperature measurements at depth affects the temperature predictions, both in the sedimentary basin and the Precambrian basement. It is an important issue to solve in further studies. Furthermore, knowledge of the thermal properties of the Precambrian basement of the Grenville Province and its geometry is poor. In addition, there is a wide confidence interval on thermal conductivity of specific lithologies in the Anticosti sedimentary basin. It has a significant impact on temperature predictions at depth and should be improved for studies focusing on electricity production. Despite a wide confidence interval on temperature predictions, geothermal electricity generation from reservoirs at 120 °C or more appears difficult in the current technical and economic context. Electricity generation at a low temperature with an inlet of 70 °C could be achieved at a reservoir depth of 2–4 km, but with a net efficiency of 10–11% (considering a flow rate of 40 l s−1 and a cooling temperature of 5 °C). Direct use of geothermal heat from the deepest part of the sedimentary basin seems to be the most realistic option, provided that sufficiently permeable horizons can be found.


General context
Interest in renewable energy technologies is growing worldwide as awareness of climate change increases. Effort to develop geothermal energy, for either direct heat use or electricity generation is being made, as it provides baseload and low-carbon energy (e.g., Hammons 2011;Hou et al. 2018;Lund et al. 2008). However, development of geothermal projects requires extensive site-specific studies of the underground thermal regime, which can be predicted with a certain degree of confidence.
In the Province of Québec (Eastern Canada), sites showing a fair potential for geothermal energy production were identified in sedimentary basins of the St Lawrence Lowlands and the Gaspé Peninsula (Bédard et al. 2018;Chabot-Bergeron et al. 2016;Majorowicz and Minea 2013;Nasr et al. 2018). Isolated sedimentary basins, like the Anticosti Island, figure only in studies conducted at provincial or national scale, most likely because available data are limited (Jessop et al. 1984;Minea 2013, 2015a, b). In such off-grid regions, entirely relying on diesel for both electricity and heat generation, the need for renewable energy development is crucial. A re-evaluation of geophysical and geological data available in the context of deep geothermal energy may provide a better understanding of the type of technology to be considered for the exploitation of geothermal resources.
This study aims to provide a first assessment of the deep temperature distribution of the Anticosti sedimentary basin and the underlying Precambrian basement of the Grenville Province. The only village in this area, Port-Menier, had 218 year-round inhabitants in 2016, according to Statistics Canada (2017). The objective was to help develop renewable energy resources as an alternative to fossil fuels. Over the last 60 years, 31 oil and gas exploration wells have been drilled. Once data were made public, Bédard et al. (2014) built a 3D geological model of the Anticosti sedimentary basin including eight distinct sedimentary rock units. However, the available dataset is sparse and only contains 25 non-equilibrium temperature measurements at depth. Similar studies aiming to provide an initial assessment of geothermal resources with conductive heat transfer calculations while having access to sparse temperature data only were achieved in neighboring regions (Majorowicz and Minea 2013), but the sensitivity to input data still needs to be better quantified.
Such models commonly need a large volume of data for temperature predictions to be representative. Examples can be found with numerous studies focused on sedimentary basins benefiting from large and/or high-quality temperature datasets (Calcagno et al. 2014;Ebigbo et al. 2016;Fuchs and Balling 2016;Hofmann et al. 2014;Lenkey et al. 2017;Sippel et al. 2013;Wellmann and Reid 2014). Studies of these sedimentary basins showed that low-quality input data can make modeling results unreliable, leading to false conclusions (Fuchs and Balling 2016;Gray et al. 2012).
The hypothesis behind our work was that a conductive heat transfer model to define the thermal state of a sedimentary basin can still be built from sparse oil and gas exploration data, but it can only be reliable if the impact of the parameter confidence interval is taken into account. For this purpose, a 3D numerical conductive heat transfer model of the Anticosti Island was developed from the improved 3D geological model of Bédard et al. (2014). Thereby, a range of temperature/depth isotherms is given rather than a single value in order to avoid misleading conclusions. Furthermore, local 1D sensitivity of each thermal parameter is studied independently to highlight those that should be prioritized in further data acquisition. A cross-sensitivity study combining minimum and maximum values for these parameters (considered independent) was also conducted in 1D. While the values are specific to the Anticosti sedimentary basin, this provides a useful perspective for studies facing similar challenges associated to remote regions.

Geological setting
Anticosti Island, located in Eastern Québec, offers a window through the Anticosti Platform in the Gulf of St. Lawrence (Fig. 1). This sedimentary basin is mainly composed of carbonate rocks from the Ordovician to Silurian ages, resting unconformably on a basement of Precambrian rocks of the Grenville Province (Fig. 2). The sedimentary succession is not affected by major tectonic deformation, but it is cut by synsedimentary normal faults shaping the lower part of the succession in the subsurface, where the most important is the Jupiter fault (Bordet et al. 2010;Fig. 2). The sedimentary basin is 4000 m thick in the south, whereas it is only one kilometer thick to the north (Bédard et al. 2014;Castonguay et al. 2005;Desrochers et al. 2012). Sedimentary succession represents a transition from a passive margin environment to a foreland basin (Chi et al. 2010;Desrochers et al. 1988;Desrochers et al. 2012;Lavoie et al. 2005;Long 2007;Pinet et al. 2012). The 3D geological model of Bédard et al. (2014) divides the sequence into eight sedimentary thermal units, as well as one undifferentiated thermal unit for the entire Precambrian basement of the Grenville Province ( Fig. 2) (Bordet et al. 2010;Desrochers et al. 1988;Desrochers 2006Desrochers , 2009Desrochers , 2010Desrochers , 2012Lavoie et al. 2005;Pinet et al. 2012).

Methods
A 3D numerical conductive heat transfer model was built to assess the deep temperature distribution of the Anticosti Island, where the workflow is presented in Fig. 3. It was developed with the following assumptions: (1) available geothermal energy in the Earth's crust comes from the mantle and from the decay of radioactive elements in the crust (Perry et al. 2010); (2) the pre-existing 3D geological model of Bédard et al. (2014), which was modified to include recent geological data, is representative of the large-scale structure of the Anticosti sedimentary basin; (3) mantle depth is constant and 40.5 km underneath the Anticosti Island (average of three points assessments from LITHOPROBE, CNSN Polaris and USGSD Database programs; Schetselaar et al. 2017); (4) only conductive heat transfer is considered since radiation and natural convection can be neglected at the basin scale; (5) the impact of seasonal temperature variations can be considered negligible deep down at basin scale.
An analytical and iterative method, where temperature of each layer is calculated individually using the upper layer as a boundary condition, was initially used to extrapolate 1D temperature profile at depth. The 1D temperature models were built using following formulas which account for vertical conductive heat transfer and heat generation in the lithosphere (sedimentary basin and Precambrian basement). These equations require the knowledge of key parameters such as thermal conductivity, heat generation rate and temperature data (on the surface and at depth) to constrain the models.
where i and i+1 are superimposed cells. ef corresponds to the effective mean calculated from the surface to BHT depth (z). T (°C) is the temperature; e (m) is the thickness; λ (W m −1 K −1 ) is the thermal conductivity; Q 0 (W m −2 ) is the surface heat flux; Q i (W m −2 ) is the heat flux at the bottom of the cell I; T z (°C) is temperature data at depth; and A f (W m − 3 ) is the heat generation rate of the associated unit.
The geothermal exploration of Anticosti Island is still at an early stage, with a limited amount of available data. In these conditions, assigning a single "mean" value to all these parameters can yield important errors, as parameter variability is high. Therefore, a range for each parameter was defined in order to calculate their sensitivity with respect to temperature predictions. 1D temperature models using minimum and maximum values defined for each parameter were then computed. First, the parameters were varied one by one using references values for the others (1D local sensitivity analysis). Then the minimum and maximum values for the different Workflow for the modeling of 3D temperature distribution parameters were combined (1D cross-sensitivity analysis). It allowed to determine which of these parameters should be prioritized for further characterization and a temperature range in 1D. Finally, two combinations of minimum and maximum values for parameters were defined in a way to minimize and maximize temperature predictions at depth. They were used to give a range for the 3D temperature distribution.

Thermal conductivity of thermal units
Bulk thermal conductivity values were calculated theoretically for each thermal unit of the Anticosti sedimentary basin. It was based on a thickness-weighted harmonic mean, using literature values for lithologies of the sedimentary thermal units (Alishaev et al. 2012;Barkaoui et al. 2014;Birch and Clark 1940;Canakci et al. 2007;Cermak and Rybach 1982;Fuchs and Förster 2010;Gornov 2015;Romero et al. 2016;Sayed 2011;Schütz et al. 2012;Vélez et al. 2018;Vosteen and Schellschmidt 2003;Zheng et al. 2016). An upper and lower boundary for theoretical thermal conductivity were calculated by adding or subtracting 20% to the lithology thermal conductivity mean, as shown in Table 1. This percentage was chosen because mean value plus 20% fits the lab measurements that were conducted on samples from Anticosti (24 limestone samples and 2 sandstone samples). These laboratory results were not considered representative enough to be directly used as mean values as the crumbliest layers, associated with lower thermal conductivity, could not be collected. Four lithological zones were defined based on their distinct lithostratigraphic sequence since there are important lateral lithological variations from east to west. A synthetic log was compiled for each zone (Fig. 4, Table 2), using geological reports (Bertrand 1987;Globensky 1993) and core descriptions (INRS-Pétrole 1973a, b, c, 1974. Definition of relative proportion of main lithologies for each zone and thermal units gives a quantitative estimation of the unit's equivalent thermal conductivity.
Comparison of well logs helped to define boundaries between all four zones, which were kept parallel to the unit's dip. Local 1D sensitivity to zonation, and synthetic lithological logs are discussed in the following sections.
A single thermal conductivity value for each thermal unit and lithological zones was assigned for the low and high scenarios. An example of the thermal conductivity obtained according to the lithological composition of the thermal units is presented for the West Center zone (Table 3).    Bédard et al. (2014), Bertrand (1987), INRS-Pétrole (1974) Gascuel et al. Geotherm Energy (2020 Grenville Province (outcropping north of Anticosti Island) is mainly composed of different metamorphic and plutonic rocks (Moukhsil et al. 2017), but its lithological distribution under the Anticosti Island is unknown. For this reason, a single thermal conductivity value was assigned to the Basement thermal unit. Minimum and maximum values of 2.70 and 4.71 W m −1 K −1 were considered, based on the thermal conductivity average of metamorphic rocks rich or poor in quartz (Clauser 2006). The Basement thermal unit is 35-40 km thick in the 3D model. Therefore, it was considered that the effect of potential not taken into account intermediary layers of extreme thermal conductivity would average due to the thickness of Basement thermal unit. A value of 3.00 W m −1 K −1 , corresponding to mean value for granite and granitic gneiss according to Cermak and Rybach (1982) was taken as the reference scenario for the thermal conductivity of the Precambrian basement of the Grenville Province. Finally, thermal conductivity is known to vary depending on temperature. This was taken into account by using the empirical approach provided by Sass et al. (1992).

Heat generation rate of the sedimentary thermal units
Rocks produce internal heat by the decay of their radioactive elements, mainly uranium-238 and 235, thorium-232 and potassium-40 (Hamza 1973), heterogeneously present in all rock types (Jaupart and Mareschal 2011). This study considers internal heat production by assigning a different bulk heat generation rate (A) to each thermal unit. Heat generation rate of thermal unit was calculated via two different methods, then compared in "Results". Firstly, heat generation rate profiles were derived from uranium, thorium and potassium concentration obtained from spectral gamma ray logs of vintage oil and gas wells. For this calculation, thermal unit densities in the Ch-Ju-GR, Lower Vauréal and Upper and Lower Mingan thermal units were assumed similar to that of the Tr-BR-Ch unit of the St. Lawrence Lowlands sedimentary basin, while density of rocks from Macasty thermal unit was assumed similar to density of the Utica unit from Bédard et al. (2018). These units can be considered as analogues as they are similar in lithologies and age (time for compaction and diagenetic processes). A was calculated all along the spectral gamma ray logs (Eq. 4) using the equation of Bucker and Rybach (1996), then an average value for each thermal unit was calculated:

is the density of the thermal unit; [U] (ppm) is its concentration in uranium, [Th] (ppm) is its concentration in thorium and [K] (%) is its concentration in potassium.
Secondly, the heat generation rate was calculated from theoretical lithology values considering lithological logs (Fig. 4), as concentration in radiogenic elements is correlated with lithologies (Roque and Brenha 1996). A thickness-weighted arithmetic mean was calculated using a similar approach that was used for assigning thermal conductivity values to thermal units. The theoretical values, shown in Table 4, were based on the work of Hasterok et al. (2017). In the Anticosti sedimentary basin, shales in thermal units other than the Macasty are described as greenish and calcareous and were, therefore, assimilated to iron-rich calcareous shales. Heat generation rates for these intervals were calculated as a mix between iron-rich shale (75%) and limestone (25%). Only the Macasty thermal unit was considered similar to Hasterok et al. (2017) shale category.

Temperature data
Temperature at shallow depth is impacted by seasonal variations and snow cover insulating the ground in winter. This phenomenon is attenuated at depth and temperature remains constant in a zone named the undisturbed ground temperature by  Ouzzane et al. (2015), at about 10 m below the surface, depending on the subsurface thermal properties. As this study targets the deepest part of the sedimentary basin and its underlying Precambrian basement, seasonal temperature variations were not simulated. Thus, the undisturbed temperature at 10 m depth was considered as the ground surface temperature T 0 for this study. Values assigned are from a surface temperature assessment for the Province of Québec by Comeau et al. (2017), who used climate normal and the empirical relationship of Ouzzane et al. (2015) to calculate the undisturbed ground temperature. Bottom-hole temperature (BHT) data from 25 oil and gas exploration wells on Anticosti Island (Fig. 2) were available via the Ministère de l'Énergie et des Ressources naturelles (2019) web site. Some of the BHT data were not taken into account for this study, as they were considered unreliable (measured before 1975 or with a measured temperature lower than 28 °C, which could be a record of surface temperature as the measurements were conducted in summer). 15 BHT were selected for this study (Fig. 4). These BHT were measured a few hours to a few days after drilling. Since drilling operations cause thermal disturbances on surrounding rocks that can last for ten times as long as the drilling time (Beardsmore and Cull 2001;Jessop 1990), a correction is thus needed to estimate the equilibrium temperature.
An empirical correction established by Waples and Ramly (2001) was used in this study. It calculates thermal disturbance due to drilling according to depth and time since the end of the drilling mud circulation, which were the only parameters available: where T c (°C) is the BHT corrected for drilling disturbance, T 0 (°C) is the assumed undisturbed ground temperature, T mes (°C) is the measured BHT and f s (−) is a correction factor determined with: where TSC (h) corresponds to the time since the end of mud circulation and z (m) to the depth of the measured BHT.
Paleoclimatic variations due to glacial and interglacial periods also affect the temperature at depth (Beck 1977;Chouinard and Mareschal 2009;Jaupart and Mareschal 2011;Jessop 1990). This can significantly impact the evaluation of the geothermal gradient and heat flux. Since most of Anticosti BHT measurements were made at a depth shallower than 2000 m, paleoclimatic corrections were applied for all BHTs. The paleoclimatic correction theory was originally developed by Birch (1948), then supplemented by others (e.g., Beck 1977;Crain 1968;Jaupart and Mareschal 2011;Turcotte and Schubert 2014;Westaway and Younger 2013). In most studies, the paleoclimatic correction is directly applied to the geothermal gradient derived from a temperature profile. Here, the only data points available are from BHT data. Therefore, it was chosen to calculate a theoretical equilibrium BHT, as it would be if the climate had always been the same (Bédard et al. 2018;Jessop 1990): where ΔT (°C) is the paleoclimatic correction to be applied on temperature value at depth z (m); T t (°C) is the variation between the mean surface temperature during the considered period and the present mean surface temperature; erf is the error function; s (m 2 s −1 ) is the thermal diffusivity of rock [1.2 × 10 −6 m 2 s −1 ; (Bédard et al. 2018)]; t t1 (s) is the end time of the period considered and t t2 (s) is the beginning time of the period considered. Glacial history for the Anticosti Island was assumed similar than the nearby Sept-Îles area where ground surface temperature variations have been deduced from a temperature profile in a borehole that is over 2000 m deep.

Surface heat flux
Surface heat flux Q 0 was determined for each well with BHT data, with iterative computations of the thermal conductivity of the thermal units according to the temperature at depth using the method of Sass et al. (1992) (Fig. 5). Because rock's porosity in Anticosti sedimentary basin is generally low to non-existent, influence of pressure on porosity and therefore thermal conductivity was not considered. Surface heat flux was used in the 1D vertical temperature simulations as an initial value for heat flux when calculating temperature from the surface to the bottom of the model (based on Eqs. 1, 2, 3). It was also used to calculate theoretical values for the heat generation rate of the Precambrian basement of the Grenville Province used in the 1D models and the 3D numerical heat transfer model.

Heat generation rate of the basement thermal unit
Theoretical Basement thermal unit heat generation rates (A Pr ) were calculated from the resulting surface heat flux, sedimentary thermal units heat generation rate and mantle heat flux in 1D for each well (Step 2 in Fig. 3; Eq. 8). It corresponds to global averages for all the various lithologies of the Basement thermal unit located at depth, whose geometry is unknown.  (Mareschal and Jaupart 2013). It was taken as 15 mW m −2 for the simulation of the reference scenario. The heat generation rates calculated in 1D for each well with BHT data with the optimistic and pessimistic sets of parameters were averaged. The average for the optimistic scenario plus standard deviation was assigned as the heat generation rate for the Basement thermal unit in the optimistic 3D model. The average for the pessimistic scenario minus standard deviation was used for the pessimistic 3D model.

1D sensitivity of thermal parameters
Local sensitivity analysis and cross-sensitivity analysis were performed manually in 1D to better understand temperature predictions and possible sources of errors.
Vertical 1D temperature profiles were computed analytically for each well with BHT to interpolate temperature from the surface to the BHT measurement point and to extrapolate temperature into the Basement thermal unit. As a result, temperature and heat flux were calculated analytically along the 1D lines with a spatial resolution varying from 1 to 100 m from the surface toward mantle depth at well locations (based on Eqs. 2 and 3). The sensitivity of temperature predictions with respect to thermal parameters variability (A: lithology of the sedimentary thermal units, B: thermal conductivity assigned to the lithologies of sedimentary thermal units. C: thermal conductivity of the Precambrian basement of the Grenville Province. D: effective heat generation rate assigned to the entire sedimentary basin E: BHT correction. F: mantle heat flux) was evaluated in 1D for 15 wells by assigning a reference parameter set and varying these parameters one by one within their confidence interval. The reference parameter set comprises the synthetic lithology corresponding to the wells location (Fig. 4), the average thermal conductivity and heat generation rates estimated according to the lithologies, the corrected BHT and a mantle heat flux of 15 mW m −2 . Thermal parameters were considered to be independent as there was not enough data to study their correlations. This sensitivity analysis conducted in 1D allowed to identify the parameters with most impact to consider for further simulations defining temperature uncertainty at depth.
A cross-sensitivity analysis combining the minimum and maximum values for the six parameters identified above was then made for well DZ017 to illustrate the variability of temperatures at depth in 1D. Well DZ017 was considered the most representative out of the available dataset because the difference in corrected or uncorrected temperature gradient was close to the mean of all the wells. Parameters were considered independent. A total of 64 temperature scenarios were thus calculated using all possible combinations minimum and maximum values for parameters. The confidence interval of the parameters with most impact (i.e., corrected BHT, thermal conductivity of sedimentary lithologies and thermal conductivity of the Basement thermal unit) were combined to generate optimistic and pessimistic scenarios for predicting a range of temperature at depth. These scenarios were computed in 1D for all wells with BHT data, generating theoretical values for the heat generation rate of the Basement thermal unit.
Minimum and maximum values for parameters with most impact were then combined to produce optimistic and pessimistic 3D simulation scenarios. These two scenarios were defined to include all possible temperature distribution at depth when modeled in 3D.

3D temperature distribution
Based on the 3D geological model of the Anticosti sedimentary basin, a semi-structured finite element mesh containing about 1400,000 elements was built using FEFLOW triangle mesh generator. Each thermal unit present in the 3D geological model was meshed as one layer in the FEFLOW model. The thickness and depth of the layers vary according to the geological model, with the constraint of a minimum thickness of 0.1 m. The Basement thermal unit was extended up to a depth of 40.5 km (mantle depth; Shetselaar and Snyder 2017) and was divided into 13 layers of increasing thickness (up to 5500 m between the deepest slices). Because it causes a significant difference in thermal units' depths and thicknesses over a short distance (ex.: up to 500 m difference in depth of the Basement thermal unit), Jupiter fault was meshed as a fault zone where meshing is refined. Element diameter varies from 203 to 4895 m on the island and from 3083 to 6497 m in the surrounding zone.
No relevant information to constrain possible heat convection was available in the study area. For that reason, purely conductive heat transfer was assumed, and groundwater flow was not simulated. The simulation was run under steady state, neglecting effects of past climate change. Two temperature boundary conditions were imposed: a first type condition at the top of the model (near-surface undisturbed ground temperature equivalent to the 1D models), and a second type at the bottom of the model (mantle heat flux, originating from the value used to calibrate the heat generation rate in the 1D models). The lateral boundaries were defined as adiabatic or insulating.
Model properties were chosen according to the optimistic and pessimistic scenarios. The Basement thermal unit heat generation rate was adjusted in 1D for each well to reproduce BHT data for both the optimistic and pessimistic scenarios. The average 1D value for the pessimistic scenario minus the standard deviation was assigned as Basement thermal unit heat generation rate for the 3D pessimistic simulation. Similarly, the average value plus the standard deviation calculated in 1D was assigned for the 3D optimistic simulation. Thermal conductivity in both scenarios was then corrected iteratively for temperature effect running the model and using the approach of Sass et al. (1992), as previously done in 1D (Table 5).

Thermal conductivity of thermal units
Thermal conductivity values calculated for the optimistic and pessimistic scenarios are presented in Table 6. The dolomitic thermal unit of Romaine is the most conductive, while the shaly Macasty thermal unit presents the lowest thermal conductivity. Furthermore, for the Bs-EB-UV and Lower Vauréal thermal units, the East zone is more conductive than other zones due to lateral variation from limestone and shale to sandstone (Fig. 4).

Heat generation rate of sedimentary thermal units
Heat generation rates calculated from theoretical values were compared to those calculated from spectral gamma ray logs for each lithological zone and unit where it was possible (East Center and West Center zones; Figs. 6, 7 and 8). Results of both methods are very similar in such a way that theoretical values were then used for this study.

Temperature data
Temperature data were corrected for drilling and paleoclimatic effects (Fig. 9). BHT data were compared with available temperature data from drill stem test (DST; Fig. 10), considered to give a more realistic equilibrium temperature approximation since it is from production fluids (Förster et al. 1997;Harrison et al. 1983;Peterson 2013;Waples and Ramly 2001). Simple linear regression from DST data with a fixed surface temperature gives an expected temperature at BHT depth comprised between uncorrected and corrected BHT values. This seems to indicate that the correction used for drilling effect overestimates temperatures at depth.
As a result, corrected BHT (for both drilling disturbance and paleoclimatic effect) and uncorrected BHT values were taken as an upper and lower boundary for the calculation of maximum and minimum surface heat flux values. Minimum values stem from uncorrected BHT data and give a global gradient of 17.4 °C km −1 . Maximum values stem from corrected BHT data and show a global gradient of 25.4 °C km −1 (Table 7 and Fig. 9).

Table 5 Glacial history considered for paleoclimatic correction of the Anticosti Island
Correspond to the scenario a) of a study conducted near Sept-Îles by Mareschal et al. (1999)

Surface heat flux
Average 1D surface heat flux is 68.9 mW m −2 for the optimistic scenario, with a standard deviation of 9.4 mW m −2 (15 1D values), while the average is 33.9 mW m −2 for the pessimistic scenario with a smaller standard deviation of 2.9 mW m −2 (Table 7). According to Minea et al. (2011), mean heat flow in Québec Province is about 56.9 mW m −2 , with a standard deviation of 17.4 mW m −2 . The higher standard deviation in surface heat flux obtained for the optimistic scenario could be linked to the corrections applied to BHT data. These corrections differ from well to well depending on depth and TSC and tend to increase the standard deviation in temperature gradient (Table 7).

Heat generation rate of the basement thermal unit
Theoretical heat generation rates calculated in 1D for the optimistic scenario present an average of 1.36 µW m −3 and standard deviation of 0.28 µW m −3 (15 1D values), while the average is 0.48 µW m −3 for the pessimistic scenario with standard deviation of 0.07 µW m −3 . Because of the poor knowledge of the Basement thermal unit properties and spatial variability, it is difficult to assess the realism of these values. Though the heat generation rate calculated for the Basement thermal unit falls within values used by other authors for the Precambrian basement of the Grenville Province such as Liu et al. (2018), values of 1.64 µW m −3 for the optimistic scenario and 0.41 µW m −3 for the pessimistic scenarios were assigned to the Basement's heat generation rate in the 3D model.

1D sensitivity of thermal parameters
Results from the 1D temperature computation are available for all wells used in this study (see Additional files 1, 2 and 3). Confidence interval differs from well to well, depending on TSC (impacts BHT correction), depth of the BHT data and depth of the base of the sedimentary basin. The most sensitive parameter in 1D models depend on the depth of the temperature predictions considered (Fig. 11). Sensitivity to the lithological zone was tested to know what would be the impact of an error in the zone's delimitation (Fig. 11a). It shows that impact on modeled temperature is only significant if a well is wrongly placed in the East zone or vice versa. This is due to the higher sandstone content in the East, leading to higher thermal conductivity values and high heat flux, impacting in turn the extrapolation of temperatures at depth. Thermal conductivity of the different lithologies in the sedimentary basin further impacts the calculated heat flux and therefore temperature in the Precambrian basement (Fig. 11b). The temperature modeled in the Precambrian Fig. 10 Comparison of available DST temperature data with uncorrected and corrected BHT data. The regression line was fixed at 0 m depth according to the undisturbed ground temperature basement is also sensitive to its own thermal conductivity (Fig. 11c). The basin's heat generation rate does not have a significant impact on modeled temperatures and needs not be considered in further studies (Fig. 11d). The model is calibrated to reproduce BHT, using the surface temperature as a second control point (the Basement thermal unit heat generation rate is consequently the variable). Thus, in the upper part of the model, temperature is mostly sensitive to BHT correction. Sensitivity to BHT correction increases as temperature is extrapolated at depth (Fig. 11e). With our modeling approach, variations of surface heat flux are attributed to the Basement thermal unit heat generation rate that is calculated. Thus, temperature sensitivity to mantle heat flux is null (Figs. 11f, 12).
Diagrams in Fig. 13 highlight which parameters are the most sensitive to temperature predictions depending on the depth of the predictions. For direct use (at 2 km depth), priority should be given to obtaining more reliable value for equilibrium BHT. When trying to identify suitable depth for power generation with a minimum temperature of 120 °C, model results are mostly sensitive to BHT correction and thermal conductivity of rocks of the sedimentary basin and the underlying Precambrian basement. Knowledge of those three parameters should be improved in accordance with opportunities.
The sensitivity study conducted with multiple wells shows that the uncertainty on BHT and the resulting impact on simulated temperature decrease when TSC increases, as the measured temperature is closer to equilibrium. There is a linear correlation factor of 65% between differences in modeled temperature at 5 km depth obtained with minimum and maximum values for BHT and TSC. The impact of the thermal conductivity of the sedimentary and Basement thermal units varies with the thickness of the sedimentary basin. There is a linear correlation factor of 68% between differences in modeled temperature Table 7 BHT data corrections and input parameters ∇: refers to the temperature gradient; 0 refers to the uncorrected value; 1 to the correction for the drilling effect (true temperature at well depth) only; and 2 to the corrections for both drilling and paleoclimatic effects (theoretical temperature if the surface temperature had stayed constant through time). Q o min, and Q o max are the surface heat flux calculated for the optimistic and pessimistic scenarios. σ is for standard deviation at 5 km depth for minimum and maximum thermal conductivity of sedimentary thermal units and basin thickness. The correlation factor is 74% when varying thermal conductivity of the Basement thermal unit instead. The cross-sensitivity analysis achieved with all possible combinations of the minimum and maximum values defined for parameters shows that temperature variability generally increases with depth (Fig. 13).
Other sources of errors in the estimation of temperature at depth are related to the 3D geological model geometry (expected to have a small impact) assumption of uniform thermal properties for the entire Basement thermal unit, and assumption of conductive heat transfer only (see "Discussion" section).

Fig. 11
Temperature sensitivity with respect to thermal parameters considered computing analytical 1D profiles, using DZ017 well. Impact of a lithology from the synthetic logs of other zones are applied to the calculation of the thermal conductivity and heat generation rate in the sedimentary thermal units; b the thermal conductivity assigned to the lithologies of sedimentary thermal units; c the thermal conductivity assigned to the Basement thermal unit; d heat generation rate assigned to the entire sedimentary basin; e BHT corrections, and f mantle heat flux were evaluated separately Confidence intervals for BHT data and thermal conductivity of both sedimentary basin and the underlying Precambrian basement were combined to generate pessimistic and optimistic scenarios for the 3D numerical model (Table 6).

Fig. 12
Average relative sensitivity for all wells used in this study considering each thermal parameter to compute analytical 1D temperature at a 2 km depth and at b depth of the 120 °C isotherm. Impact is calculated from pessimistic and optimistic scenarios for each thermal parameter, while quantifying difference in temperature or isotherm depth (Fig. 11)   Fig. 13 1D cross-sensitivity analysis between all thermal parameters for temperature evaluated at 2 and 5 km depth in well DZ017 using all possible combinations of pessimistic and optimistic scenarios. The orange points around 1400 m depth are the corrected and uncorrected BHT. Blue points represent computed temperature in the different scenarios

3D temperature distribution
Temperature modeled numerically with the optimistic and pessimistic scenarios were, respectively, compared to corrected and uncorrected BHT data. The optimistic simulation scenario fits the upper boundary for BHT temperature with a root mean square error (RMSE) of 9.38 °C, with simulated temperature being higher than BHT for most wells. However, the temperature simulated for wells DZ019 and D020 remains lower than their BHT. Pessimistic simulation scenario fits the lower boundary for BHT with an RMSE of 3.62 °C. For most wells, modeled temperatures are lower than BHT, however, wells DZ009 and D013 show the opposite. This difference between modeled temperature and BHT in both scenarios is mainly due to the addition or subtraction of standard deviation to mean theoretical Basement thermal unit heat generation rate (calculated in 1D) before implementation in the 3D model. Differences for wells DZ019 and D020 in the optimistic scenario, and DZ009 and D013 in the pessimistic one, can be explained by the spatial variability in heat generation and thermal conductivity of the Basement thermal unit, not considered in numerical simulations. Rocks of the Precambrian basement on the coast north of Anticosti Island, used as an analogue to the Basement thermal unit, are indeed known to have a great variability in heat generation rate, with prospects for uranium mining (Lavallée 2010). Moreover, a magnetic body in the upper Precambrian basement has been identified near well DZ019 (Pinet et al. 2012), which can be linked to the higher than expected BHT.
A good agreement was found between the result of the pessimistic 3D simulation and the pessimistic temperature profile computed for the 1D cross-sensitivity analysis for well DZ017 (41 °C and 40 °C at 2 km depth, respectively). However, the optimistic 3D simulation resulted in hotter temperature when compared to the 1D cross-sensitivity analysis (63 °C and 55 °C at 2 km depth, respectively). The same trend is observed at greater depth. This can be explained by the homogeneous Basement thermal unit heat generation rate assigned for the 3D model. Heat generation value used in the pessimistic 3D simulation is close to the one computed for well DZ017 (0.410 µW m −3 in 3D compared to 0.406 µW m −3 for well DZ017). However, the value used for the 3D optimistic scenario (1.64 µW m −3 ) is greater than the one calculated in 1D for DZ017 (1.29 µW m −3 ). Different options can be considered for further geothermal energy development on Anticosti Island. Producing geothermal electricity from high-to-medium temperature (over 120 °C) would require very deep wells and enhanced geothermal systems (EGS) or closed loop, as 120 °C is only reached between 3.8 km (optimistic scenario) and 11 km depth (pessimistic scenario) inside the Basement thermal unit (Fig. 14). Confidence intervals need to be reduced by further analysis, but such systems would range from very costly to not feasible in current technical conditions. Electricity generation with a low temperature system can be considered. According to Climeon et al. (2017), electricity can be generated from temperature as low as 70 °C. Assuming a cooling temperature of 5 °C, it would yield an electricity output per module of less than 80 kW when having a flow rate of 40 l/s for a net efficiency less than 11%. The optimistic simulation scenario indicates that 70 °C can be reached in the Romaine thermal unit in the southwest part of the island. As the basal Romaine thermal unit is known to have secondary porosity and permeability due to diagenetic processes (secondary porosity up to 30% in some meter-scale strata (Lavoie et al. 2005), it may be permeable enough to implement a geothermal doublet without EGS. However, at the current state of knowledge, implementation of such a system presents a double risk. Wells in the Romaine thermal unit may not be permeable enough and/or the 70 °C isotherm may be reached in the upper Basement thermal unit only (4 km depth at the shallowest in the pessimistic scenario), requiring EGS technology in this case (Figs. 14,15).
District heating in Port-Menier (main village on the island) seems possible in both optimistic and pessimistic scenarios with a temperature of 40-55 °C reached at the base of the sedimentary basin at this location (Fig. 16). As the temperature is moderate, heat exchangers and heat pumps may be necessary to reach a useful temperature (Glassley 2010). Geothermal doublets or closed loop systems could be developed depending on permeability. This geothermal resource could also be used for greenhouses heating, generating new economic activity and further reducing the island's dependence on imports.

Discussion
The assessment of temperature prediction sensitivity with respect to thermal parameters variability was shown to be important while estimating the geothermal resource potential of a sedimentary basin. This step is, however, rarely conducted in such modeling studies. Most often, sensitivity to modeled temperature is only estimated when a large and/or good-quality database is available (1243 temperature data points in Lenkey et al. (2017); an equilibrium temperature profile and 148 BHT from 44 wells in Fuchs and Balling (2016); 7 equilibrium temperature data for Ebigbo et al. (2016)). Moreover, sensitivity analysis made with regional models is commonly conducted in 3D. Several recent studies used a geostatistical approach to produce equiprobable distributions of thermal properties, which were then used for sensitivity analysis (Camp et al. 2018;González-Garcia and Jessell 2016). Models can additionally be calibrated by automatic optimization to reproduce measured temperature data (Fuchs and Balling 2016). Such approaches can help define realistic models allowing for a small confidence interval, but require a prior knowledge of the spatial variability of the rock's properties. An important amount of computation power and time is additionally needed, especially if several parameters are varied. Such complex sensitivity analysis can only be made at a late stage of geothermal exploration and could not be realized at Anticosti Island. The rock properties or boundary conditions were varied one by one in other studies to evaluate their sensitivity (Della Vedova et al. 2008;Ebigbo et al. 2016;Frick et al. 2015;Noack et al. 2012). While this approach can be more similar to what was done in the present study, it differs by the absence of calibration process in the work flow. The temperature data at depth were used for comparison with modeled temperature only, while in the present study, temperature data are used to calibrate the Basement thermal unit heat generation rate for each scenario. Only one to three parameters were varied in these studies, except for the work of Della Vedova et al. (2008), which can be due to the important time needed to complete a full-sensitivity analysis in 3D. As a consequence, a new approach had to be defined for Anticosti Island to rigorously consider parameter sensitivity while facing the challenges caused by data sparsity common to remote regions. The originality of the present work relies on the sensitivity analysis conducted at the early stage of geothermal exploration on multiple thermal parameters. Sensitivity to thermal conductivity (in both sedimentary basin and Precambrian basement), heat generation rate in the sedimentary basin, BHT data used to constrain the Precambrian basement heat generation rate and the mantle heat flux was first studied in 1D, allowing to rapidly define temperature uncertainty. Based on this first analysis, the parameters with most impact were identified and the change of sensitivity with depth as well as the impact of TSC, sedimentary basin thickness and BHT depth was quantified. Minimum and maximum values of the parameters having the most impact (BHT correction, thermal conductivity of sedimentary rocks and thermal conductivity of the Basement thermal unit) were then combined to build optimistic and pessimistic simulation scenarios. Simulations using these scenarios were run in 3D providing results enclosing all possible temperature variations at depth.
Results of this sensitivity study can be used as a guide for other studies suffering from a similar lack of data. It highlights the parameters that must be prioritized for in depth characterization depending on the depth of target resources. Furthermore, it shows caution by giving a range of possible temperature modeled at an early stage of exploration, rather than providing a unique scenario.
Variability of the 3D geological model geometry, constrained by 24 wells only, was, however, not considered in the evaluation of temperature sensitivity. Despite a seemingly simple architecture of the sedimentary basin, significant adjustment was necessary when the latest well data were added to the model of Bédard et al. (2014). It stands to reason that this 3D geological model is not an exact representation of subsurface architecture and could be improved by adding further well data.
Thermal conductivity values for lithologies of sedimentary thermal units were taken from literature, adding ± 20% variability (Table 6). Sensitivity analysis reveals an important impact on temperature predicted at depth (Fig. 11b). This corroborates the work of Rauch et al. (2018), which indicates that previous studies on Pennsylvania Appalachian Basin overestimated the heat flow by as much as 50% because of inaccurate extrapolation of thermal conductivity. The confidence interval could be reduced by laboratory analysis on each lithology identified in well logs. This would take a fairly large number of samples as lithological variations are observed in each unit. Another option would be to use well log analysis, relating well log signals to thermal conductivity, as was done in several recent geothermal studies (e.g., Fuchs 2018;Fuchs and Balling 2016;Lenkey et al. 2017;Nasr et al. 2018;Sippel et al. 2013). Accurate calibration with associated drilling cores would be needed to avoid bias.
Assigning different values of heat generation rate to the different sedimentary thermal units and zones appears unnecessary at this stage of geothermal exploration. It was indeed shown in the 1D local sensitivity analysis that even when assigning single extreme values for the whole sedimentary basin, impact on temperature predictions was much less than that resulting from the other thermal parameters (Figs. 11, 12).
The assumption of uniform thermal properties for the Basement thermal unit is greatly over-simplified, as it is known that the Grenville Province shows high variability in thermal properties Moukhsil et al. 2017). The amount of data points obtained at each well was judged insufficient to determine a reliable distribution of Basement thermal unit heat generation rate with geostatistical methods. We added or subtracted standard deviation to the average value calculated in 1D to take this spatial variability into account by widening the range of predictions. One option to improve representation of the Precambrian basement would be to combine information obtained at wells from 1D temperature profile to gravimetric and magnetotelluric data analysis to define the spatial variations of its properties. Moreover, vertical averages for heat generation rate in the sedimentary basin were calculated in 1D for each BHT point. Yet, the effect of heat generation rate variations on the resulting surface heat flux value decreases when the depth of the heat source increases. Variability in surface heat flux is therefore linked to changes of thermal properties in the upper part of the Basement thermal unit. This calculation could be revised for the spatial variability observed in surface heat flux to be linked with variations of heat generation rate in the upper Basement thermal unit only. Such work definitely needs to be supported by geophysical data analysis such as gravity, magnetotelluric fields and seismic velocity to further improve Precambrian basement knowledge (Iovenitti et al. 2016).
The lack of equilibrium temperature data, or information to use for a more robust correction of drilling disturbance effect constitutes one of the main reasons the confidence interval on temperature predicted in this study is so wide. Recording equilibrium temperature profiles must become a priority when assessing the geothermal potential of a region for which good-quality temperature data at depth is not already available. If it is not possible, then at least obtaining more information on the conditions of the drilling (duration and speed of the drilling mud circulation, mud thermal properties, etc.) and temperature measurement (depth, TSC, confidence range of the instrument) and/or having access to several measurements per well would allow for more reliable temperature corrections.
The 1D cross-sensitivity analysis was run using only minimum and maximum values for parameters. While it gives a range of possible temperature at depth, it combines extreme values for parameters, generating scenarios that are unlikely to occur in real life. Furthermore, the parameters were considered independent. Similarly, scenarios simulated in 3D were defined with a combination of extreme values, considered independent. There are some correlations between parameters (e.g., BHT correction and thermal conductivity), but they could not be quantified in this study. Further study can include statistical repartitions and correlations of parameters to generate more realistic scenarios.
Results of this study show a high sensitivity to thermal parameters at depth. Both optimistic and pessimistic scenarios adopted for the 3D simulations are unlikely as they combine extreme values. Automatic calibration of model parameters used in several geothermal studies (Fuchs and Balling 2016;Hardwick et al. 2014;Wellmann and Reid 2014) should be used in a next step to reduce the range of parameters considered in the optimistic and pessimistic scenarios, especially for the Precambrian basement, unattainable by direct characterization methods. However, this can only be possible with equilibrium temperature data to improve BHT correction.

Conclusions
This study presents a method to analyze information available from a small and sparse dataset to predict temperature distribution at depth for remote sedimentary basins and assess the geothermal resource potential. It also allowed to evaluate the temperature sensitivity with respect to the variability of input thermal parameters. It can be used as a template for other studies suffering from the scarcity of data. It highlights the parameters to focus on, depending on the objectives, and shows why it is important to consider confidence intervals of thermal parameters and their impact on modeled temperatures.