Climate trends in northern Ontario and Quebec from borehole temperature profiles

The ground surface temperature histories of the past 500 years were reconstructed at 10 sites containing 18 boreholes in northeastern Canada. The boreholes, between 400 and 800 m deep, are located north of 51N, and west and east of James Bay in northern Ontario and Quebec. We find that both sides of James Bay have experienced similar ground surface temperature histories with a warming of ∼1-2 K for the last 150 years, similar to borehole reconstructions for the southern portion of the Superior Province and in agreement with available proxy data. A cooling period corresponding to the little ice age was 5 found at only one site. Despite permafrost maps locating the sites in a region of discontinuous permafrost, the ground surface temperature histories suggest that the entire region is and was free of permafrost for the past 500 years. This could be the result of air surface temperature interpolation used in permafrost models being unsuitable to represent the spatial variability of ground temperatures along with an offset between ground and air surface temperatures due to the snow cover.


Introduction
Earth's subsurface thermal regime is governed by the outflow of heat from the interior and by temporal variations in ground surface temperature (GST).The heat flux from the interior of the Earth varies on timescales of the order of a few millions of years in active tectonic regions and several 100 Myr in stable continents.It can be considered as steady state relative to the timescale of climatic surface temperature varia-tions.To determine the heat flow from the Earth's interior, temperature-depth profiles are measured in boreholes.In homogeneous rocks with no heat production, the steady-state temperature profile linearly increases with depth.Persistent temporal changes in the ground surface energy balance cause variations of the ground temperature that diffuse downwards and are recorded as temperature anomalies superimposed on the linear steady-state geotherm (e.g.Hotchkiss and Ingersoll, 1934;Birch, 1948;Beck, 1977).The extent to which the ground surface temperature changes are recorded is proportional to their duration and amplitude and inversely proportional to the time when they occurred.For periodic oscillations of the surface temperature, the temperature is propagated downward as a damped wave.The amplitude of the wave decreases exponentially with depth over a length scale δ (skin depth) proportional to the square root of the period (δ = √ κT /π), where κ is the thermal diffusivity of the rock, ≈ 10 −6 m 2 s −1 or ≈ 31.5 m 2 yr −1 .This damping removes the high-frequency variability that is present in meteorological records and allows for the preservation of the long-term climatic trends in the ground temperature signal (e.g.Beltrami and Mareschal, 1995).
From the interpretation of temperature-depth profiles, it is possible to infer centennial trends in Earth's surface temperature variations.The first attempt to infer climate history from temperature-depth profiles was the study by Hotchkiss and Ingersoll (1934), who estimated the timing of the ice retreat at the end of the last glaciation.It was, however, not until the 1970s that systematic studies were undertaken to infer past climate from such profiles (e.g.Cermak, 1971;Sass et al., 1971;Vasseur et al., 1983).In the 1980s, with increasing concern over global warming, use of borehole temperature-depth profiles to estimate recent (< 300 years) climate change became widespread following the study of Lachenbruch and Marshall (1986).This has lead to many local, regional, and global studies (e.g.Huang et al., 2000;Harris and Chapman, 2001;Gosselin and Mareschal, 2003;Beltrami and Bourlon, 2004;Pollack and Smerdon, 2004;Chouinard et al., 2007;Pickler et al., 2016).
Because of the availability of suitable temperature-depth profiles in the Canadian Shield, many studies have been undertaken in central and eastern Canada.The majority of these studies have used temperature-depth profiles from the southern portion of the Superior Province of the Canadian Shield (∼ 45-50 • N), where many mining exploration holes are readily available and the crystalline rocks are less likely to be affected by groundwater flow than sedimentary rocks.These studies have shown a warming signal of ∼ 1-2 K over the last ∼ 150-200 years, following a period of cooling about 200-500 yr BP associated with the Little Ice Age (LIA) (e.g.Beltrami and Mareschal, 1992;Wang et al., 1992;Guillou-Frottier et al., 1998;Gosselin and Mareschal, 2003;Chouinard and Mareschal, 2007).
For logistical reasons, mining exploration has been restricted to the southernmost part of the Shield and the few holes that have been drilled in northern regions cannot be measured because they are blocked by permafrost.Nevertheless, a few studies were conducted at higher latitudes.Majorowicz et al. (2004) reconstructed the GST history for 61 temperature-depth profiles between 60 and 82 • N in northern Canada.They found strong evidence that GST warming started in the late 18th century and continued until present.Simultaneous inversion of their data yielded a warming of ∼ 2 K for the last 500 years.Studies on Ellesmere Island (above 60 • N) have shown varying trends, confirming that temperatures do not increase uniformly over Arctic regions.Taylor et al. (2006) reconstructed the 500-year GST history from three boreholes and found a 3 K warming since the LIA minimum, ∼ 200 yr BP, which is consistent with Beltrami and Taylor (1995) results and the oxygen isotopes studies on ice cores from the region (Fisher and Koerner, 1994).Chouinard et al. (2007) used three temperature-depth profiles in a region with continuous permafrost at the northernmost tip of Québec to infer the GST history.They found a very strong and recent warming of ∼ 2.5 K, with the largest part of this warming occurring in the preceding 15 years, i.e. much later than on Ellesmere Island.Because of lack of adequate borehole temperature depth profiles in eastern Canada between 51 and 60 • N, the large region between the Canadian Arctic and the southern part of the Canadian Shield has not been studied and the climate trends of the last 500 years for this region remain unclear except for boreholes at Voisey Bay, at 56 • N on the east coast of Labrador, which show almost no climate signal (Mareschal et al., 2000).
The first motivation of this study is to reduce the gap in data between the Arctic and southeastern Canada.We shall examine 18 temperature-depth profiles measured at 10 sites from eastern Canada to reconstruct the GST histories for the last 500 years.The sites are located in the poorly sampled region north of 51 • N, west and east of James Bay in northern Ontario and Québec.They are to the north of the previous eastern Canada studies and south of the Arctic ones, in a part of the Superior Province where heat flux is extremely low (< 30 mW m −2 ) (Jaupart et al., 2014).
The second motivation of the study is to assess whether borehole temperature profiles can be used to retrace the evolution of permafrost in northern Ontario and Québec.Permafrost maps locate the boreholes in a region of discontinuous permafrost (Brown et al., 2002).All sites, excluding Noront, lie in discontinuous isolated patches of permafrost, i.e.where less than 10 % of the ground is frozen.Noront lies near the southern boundary of extensive discontinuous permafrost, i.e.where permafrost affects between 50 and 90 % of the ground.In regions with an absence of ground temperature measurements, such as northern Ontario and Québec, permafrost maps are estimated from surface air temperature and their contour lines (Heginbottom, 2002).The −2.5 • C mean annual surface air temperature (SAT) contour line for the period 1950-1980 crosses the southern part of our study region in Ontario, and most of the Québec sites are located between −2.5 and −5 • C SAT contour lines (Phillips, 2002).However, permafrost was not encountered during sampling of the Québec or Ontario boreholes.It is also worth pointing out that the ground is covered by thick snow cover during several months (from mid-December to late April) in the regions above 50 • N. Studies demonstrated that the ground surface temperatures are strongly affected by the duration of the snow cover and are offset from SAT (Bartlett et al., 2005;Zhang, 2005;González-Rouco et al., 2006, 2009;García-García et al., 2016).In these regions with extensive snow cover, the borehole temperature profiles are affected by changes in both SAT and snow cover.Meteorological and proxy data indicate that there is more snowfall and longer snow cover on the ground in the Québec region than in Ontario (Bégin, 2000;Brown and Mote, 2009;Brown, 2010;Environment Canada, 2010;Nicault et al., 2014).This points to possibly warmer present ground surface temperatures and smaller permafrost extent in northern Québec than in Ontario and the prospect for different ground surface temperature histories between the regions.

Theory
Assuming Earth is a half-space where physical properties only vary with depth, the temperature at depth z, T (z), can be written as (Jaupart and Mareschal, 2011) where T 0 is the reference surface temperature, Q o the reference surface heat flux, the integral accounts for the vertical distribution of heat producing elements H (z), and T t (z) is the temperature perturbation at depth z due to time-varying changes to the surface boundary condition.The thermal depth R(z) is defined as where λ is the thermal conductivity.
The temperature perturbation can be calculated by the following equation (Carslaw and Jaeger, 1959): where κ is thermal diffusivity and T o (t) is the surface temperature at time t before present.For a step change in surface temperature, T , at time t before present, the temperature perturbation T t (z) is given by Carslaw and Jaeger (1959): where erfc is the complementary error function.If the GST perturbations are approximated by their mean values T k during K intervals (t k−1 , t k ), the temperature perturbation is written as follows: T k is the average difference between the ground surface temperature during the time interval (t k−1 , t k ) and the reference surface temperature T 0 .

Inversion
To reconstruct the GST history for each temperature-depth profile, we must invert Eq. ( 5).The inversion involves solving for the parameters T o , Q o , and T k of the temperaturedepth profile.Equation (5) yields a system of linear equations in the unknown parameters for each depth where temperature has been measured.If N temperature measurements were made in the borehole, a system of N linear equations with K + 2 unknowns, T o , Q o , and the K values of T k is obtained.However, this system of equations is illconditioned and its solution is unstable to small perturbations in the temperature data; i.e. a small error in the data results in a very large error in the solution (Lanczos, 1961).Different inversion methods are available to stabilize (regularize) the solution of ill-posed problems (Backus-Gilbert method, Tikhonov regularization algorithm, Bayesian methods, singular value decomposition, Monte Carlo methods).All these inversion techniques have been applied to reconstruct the GST history (e.g.Vasseur et al., 1983;Nielsen and Beck, 1989;Shen and Beck, 1991;Mareschal and Beltrami, 1992;Clauser and Mareschal, 1995;Mareschal et al., 1999).In this paper, we have used the singular value decomposition because it is a very simple method to reduce the impact of noise and errors on the solution (Lanczos, 1961).This technique is well documented for geophysical studies (Jackson, 1972;Menke, 1989) and its application for inversion of the ground temperature history is straightforward (Mareschal and Beltrami, 1992).
For sites including several boreholes with similar surface conditions, the data are inverted simultaneously because it is assumed that they have experienced the same surface temperature variations and therefore consistent subsurface temperature anomalies.It was expected that consistent trends in the temperature profiles would reinforce each other while errors and random noise would cancel each other.However, the resulting improvement in the signal to noise ratio remains marginal unless a sufficiently large number of profiles with the same GST history are available, which is almost never the case.Simultaneous inversion is described in detail and discussed by Beltrami and Mareschal (1992), Clauser andMareschal (1995), andBeltrami et al. (1997), among others.

Description of data
Figure 1 shows the locations of the thirteen sites including 25 boreholes across northern Ontario and Québec.The heat flow of these sites has previously been studied and a detailed description of the measurement techniques and sites can be found in the heat flow publications (Jessop, 1968;Jessop and Lewis, 1978;Lévy et al., 2010;Jaupart et al., 2014).Their temperature-depth profiles can be found in Jaume-Santero et al. (2016).All the sites are located north of 51 • N, west and east of James Bay, and the boreholes range in depth between 400 and 800 m.All the holes were cased in the upper ∼ 10-15 m and, excluding Otoskwin and Nielsen Island, were drilled for mining exploration purposes.The temperature was measured at 10 m intervals using a calibrated thermistor.The sampling rate is higher at Otoskwin with temperature measurements every 1 m, while Nielsen Island was measured every 30 m.The overall accuracy is estimated on the order of 0.02 K with a precision of greater than 0.005 K. Thermal conductivity was measured on core samples by the method of divided bars (Misener and Beck, 1960).Radiogenic heat production measurements were also made on core samples but are not needed for corrections because the holes are not deep and the heat production rate is low.Only 18 holes proved suitable for inversion of the ground surface temperature.Their location and depth can be found in Table 1.Three northern Ontario sites are located in a region of discontinuous isolated patches of permafrost: Musselwhite, Thierry Mine, and Otoskwin.One site, Noront, lies near the southern edge of a region with extensive discontinuous permafrost.Thierry Mine (0605, 0606, 0608) and Noront (1012, 1013, 1014, 1015) include several boreholes.Six sites are located in northern Québec in a region of discontinuous isolated patches of permafrost (Nielsen Island, La Grande, Eastmain, Eleonore, Corvet, Camp Coulon), with Eastmain (0803, 0804) and Camp Coulon (0712, 0713, 0714) having multiple boreholes.Systematic variations in thermal conductivity observed at Nielsen Island were corrected by using the thermal depth (Bullard, 1939).Some measurements were not used for this study for different reasons (Table 2).Boreholes less than 300 m deep were rejected for being too shallow.When the mean distance to a lake was less than the depth of the hole, or less than 300 m, they were rejected.Furthermore, boreholes were deemed too steep and were rejected when they had a slope of 5 % or more over a distance comparable to the depth.The profile at Miminiska Lake (Ontario) is too shallow to be inverted, the boreholes at Clearwater (Québec) are plunging under a lake that affects the temperature pro-files, and the borehole at Poste Lemoyne is on the side of a very steep hill and the profile is seriously perturbed by the topography.We also discarded one of the temperature profiles at the Eleonore site because the borehole was plunging under a recently filled water reservoir and one of the profiles at the La Grande site because it was a few metres away from the edge of a 30 m cliff.The borehole temperature-depth profiles at sites with multiple boreholes were truncated at the depth of the shallowest borehole to ensure that the same period of time was being studied (Thierry Mine at 530 m, Noront at 400 m, Eastmain at 400 m, and Camp Coulon at 400 m) (Beltrami et al., 2011).The temperature anomaly for each site was calculated by subtracting from the data the estimated steady-state temperature obtained by least-square fitting of a linear function to the bottom 100 m of the profile (Figs.2-3).Tests were made to show that, below 300 m, the heat flux does not vary with the selected depth interval (Lévy et al., 2010;Jaupart et al., 2014) and that the reference temperature profile is stable.

Results
The temperature-depth profiles from the 10 sites were inverted to reconstruct the GST histories for the last 500 years divided in intervals of 20 years (Figs.4-6).Simultaneous inversion was used at sites with multiple boreholes: Thierry Mine, Noront, Eastmain, and Camp Coulon.The cutoff value or number of eigenvalues determines which part of the solution is eliminated to reduce the impact of noise.A lower cutoff value results in higher resolution in the reconstruction of the GST but at the expense of stability (Mareschal and Beltrami, 1992).Three eigenvalues (0.2 cutoff) were retained for all the sites except Otoskwin and Corvet, where four eigenvalues (0.08 cutoff) were retained.The results of the inversions are summarized in Table 3.Although the ground surface temperature histories differ in their details, they consistently show a trend of warming relative to the reference temperature (i.e.temperature 500 years before logging).Only one site shows indications that the GST was affected by the LIA, a cold period that occurred between 200 and 500 yr BP.
In northern Ontario, the trends of the inferred GST differ between sites.Otoskwin is the only site to show a LIA signal, with a cooling of ∼ 0.5 K with respect to the reference temperature (500 yr BP).Evidence for this cooling can be found in the temperature anomaly at ∼ 200 m but is also observed at Musselwhite and TM0608, where no LIA is reconstructed (Fig. 2).Moreover, there is a noticeable change in the Otoskwin temperature gradient at ∼ 200 m, which cannot be correlated to variations in thermal conductivity measured on 80 samples from the borehole.The recent warming at Musselwhite and Otoskwin occurred around the same time but it is observed earlier (∼ 250-300 yr BP) at Thierry Mine and Noront (Fig. 4).The total amplitude of warming differs greatly between the sites: 0.50 K with respect to refer-Table 1. Location and technical information concerning the boreholes used in this study, where rue depth is the depth corrected for the dip of the borehole, λ is the thermal conductivity, Q is the heat flux, and Q corr is the heat flux corrected for post-glacial warming.The heat flux for the site is represented by the bold values.The heat flux for each individual hole is in roman.

Site
Log ID Latitude Longitude Dip True depth Unlike for northern Ontario, a LIA signal was not found for any of the northern Québec sites (Figs.5-6).A LIA signal was expected because pollen data have suggested that the cooling during the LIA (up to −0.3 • C for North America) was strongest in northern Québec (Gajewski, 1988;Viau and Gajewski, 2009;Viau et al., 2012).However, a cooling signal at ∼ 200 m, which could be associated with the LIA, is observed in the temperature anomaly of CC0713 (Fig. 3).Discontinuities are observed in the temperature anomaly of CC0712 between 100 and 300 m.These are also observed in the temperature gradient and could be due to small water flows.The onset of the recent climate warming is the same for all the sites (∼ 100-150 yr BP), except Eleonore, where it began ∼ 200-300 yr BP.The amplitude of the warming varies between 0.5 and 2 K (Figs. 5-6), with the largest warming occurring at Corvet (2.18 K).

Discussion and conclusions
Borehole temperature profiles in northern Ontario and Québec consistently show a ground surface temperature increase of 1.51 ± 0.76 K above the reference temperature.Most of this increase took place for the period of 1850 to 2000.Thierry Mine shows larger than average warming signals: 2.85 K.The area around the Thierry Mine boreholes (0605, 0606, 0608) was cleared in the 1940s after the first opening of the mine and a satellite image locates all three boreholes ∼ 300 m from a lake.Lakes disturb a profile if they are at a distance less than the depth of the boreholes (Lewis and Wang, 1992).The proximity to the lake along with the change in vegetation cover could explain the enhanced warming signal (Lewis and Wang, 1998;Lewis, 1998).This illustrates the significant influence of non-climatic effects on ground surface temperature reconstructions from borehole temperature-depth profiles.
A cooling period corresponding to the LIA was found for only one site, Otoskwin (ON), which exhibits marked perturbations of the temperature profile.While spatial and temporal variation in the LIA have been noted (Matthews and Briffa, 2005), the absence of a consistent LIA signal in northern Ontario and Québec deserves some discussion.The LIA cooling The black line represents the best linear fit, while the pink and blue lines are the upper and lower bounds, respectively, of the 2σ confidence intervals.For N1014, TM0606, and Otoskwin, the upper and lower bounds of the confidence interval are not visible due to the temperature scale.The temperature anomaly at Musselwhite was cut at 650 m. period has been inferred from different proxies and selected borehole temperature-depth profiles in eastern Canada (e.g.Archambault and Bergeron, 1992;Beltrami and Mareschal, 1992;Wang and Lewis, 1992;Chouinard et al., 2007;Bunbury et al., 2012).For example, pollen data indicate a pronounced LIA cooling in Québec (Viau and Gajewski, 2009;Gajewski, 1988).The lack of LIA signal in the majority of the borehole inversions could be related to a combination of several factors.One is the limited resolution of the inversion of borehole temperature profiles.In the presence of noise, a period of weak cooling between 500 and 200 yr BP followed by strong warming is difficult to resolve.Resolution at Otoskwin is better because the singular value cutoff The anomaly is obtained by subtracting the estimated steady-state geotherm obtained by the least-square fit of a straight line to the bottom 100 m of the borehole temperature-depth profile.The black line represents the best linear fit, while the pink and blue lines represent the upper and lower bounds, respectively, of the 2σ confidence intervals.For N1014, TM0606, and Otoskwin, the upper and lower bounds of the confidence interval are not visible due to the temperature scale.The temperature anomaly at Nielsen Island was cut at 600 m. was lowest at this site.Furthermore, Mareschal and Beltrami (1992) showed that resolution decreases when noise and errors must be filtered and a higher singular value cutoff required to reduce the impact of noise but retain the gross features of the solution.Profiles and anomalies of Musselwhite and TM0608 are noisy, this could explain the absence of LIA signal in the reconstructions.However, the CC0713 temperature anomaly is less noisy and shows a mild cooling of ≤ 0.2 K.A test run with 1 K cooling between 1600 and 1800 and varying singular value cutoff showed that, in noisefree synthetic data, a 1 K cooling cannot be resolved with less than five singular values.This explains why a mild cooling, such as that observed at CC0713 (≤ 0.2 K) and in pollen data (≤ 0.3 K) (Gajewski, 1988;Viau and Gajewski, 2009; Years CE Figure 4. GST histories for the northern Ontario sites determined by inversion of the anomalies.For multiple holes at a given site (Thierry mine and Noront), simultaneous inversion was used.The pink and blue lines represent the inversions of the upper and lower bounds of the anomaly.For Otoskwin, the three lines are superposed.Viau et al., 2012), could not be resolved.Also, Chouinard and Mareschal (2007) suggested that the LIA could have started ∼ 100 years earlier in northern Québec than in southern Canada.Resolving the LIA would require a borehole deeper than ∼ 400 m, which is not the case of all the holes.Let us also point out that the sampling resolution at Nielsen Island is low, with measurements only every 30 m, which is not sufficient to resolve a LIA signal.While the absence of LIA signal in Québec was unexpected, its absence in northern Ontario confirms the findings of Gosselin and Mareschal (2003), who found only two sites with a LIA signal among 33 temperature-depth profiles from northwestern Ontario.They hypothesized that the lack of LIA signal could be due to the influence of Lake Superior because the two sites with LIA signals were above 50 • N and the furthest from Lake Superior.This is not supported by the present study because the four Ontario sites are several hundreds of kilometres away from Lake Superior.It is also possible that the LIA signal is masked by other physical effects, such as an advance and retreat of permafrost or a change in the precipitation regime and the duration of the ground snow cover during the LIA.
No geographic trends in the GST histories were observed, despite different SAT conditions.Meteorological data from the NOAA weekly dataset and eight general circulation models (GCMs) for the period of 1970-1999 display a longer snow cover duration in northern Québec than in northern Ontario (Brown and Mote, 2009).The higher precipitation is confirmed by proxy reconstructions of lake levels and tree forms (Bégin, 2000;Lavoie and Payette, 1992).Because of the greater snowfall and longer snow cover, the present ground surface relative to air surface temperatures in Québec are warmer than in Ontario.However, these dissimilar conditions have not resulted in noticeable discordance between the GST histories between northern Ontario and Québec, suggesting that the same differences in precipitation persisted throughout the period reconstructed.
The magnitude of the recent warming is about the same as the ∼ 1-2 K warming for the period of 1850 to 2000 inferred from several studies in the southern portion of the Superior Province (Beltrami and Mareschal, 1992;Shen and Beck, 1992;Chouinard and Mareschal, 2007) and less than the very pronounced warming in the eastern Canadian Arctic  (Beltrami and Taylor, 1995;Taylor et al., 2006;Chouinard et al., 2007).The sites are located in a region described as discontinuous permafrost, where ground temperatures are slightly below freezing, at least according to the Canadian and world permafrost maps.No sign of permafrost was found at any of the measured sites nor at the sites that were excluded (see Table 2); however, permafrost has been reported in the James Bay lowlands near Eleonore, Eastmain, and La Grande (Thibault and Payette, 2009).Not only are the present average ground surface temperatures well above the freezing point of water, but, except for Nielsen Island, the ground surface temperature histories retrieved from inversion also re-veal that the temperature has remained well above the freezing point for the last 500 years, indicating that the potential for permafrost was minimal to absent.We have also noted that during logging of more than 100 holes in the regions of Manitoba, Saskatchewan, and northern Ontario classified as discontinuous permafrost, permafrost has been encountered at only one hole, north of the town of Lynn Lake in northern Manitoba (Guillou-Frottier et al., 1998).Clearly, the spatial distribution of permafrost outlined in the available permafrost maps is questionable, possibly because they are not based on sufficiently deep ground temperature measurements but estimated from interpolated sparse records of SAT (Heginbottom, 2002;Gruber, 2012).The discrepancy between permafrost maps and direct field observations reveal that SAT interpolations are unsuitable to estimate the spatial variations of ground temperatures.This is likely because the maximum thickness of snow exceeds 1 m at the end of the winter remaining on the ground from mid-December to mid-April, resulting in a large offset between the GST and SAT (Zhang, 2005;Grosse et al., 2016).Our study suggests that borehole temperature profiles could be used in the future to assess the reality of the permafrost retreat assumed to have occurred after the LIA (Halsey et al., 1995;Schuur et al., 2008).Furthermore, borehole temperature profiles might be a better means for determining the southern extent of areas of past and present permafrost than current permafrost maps and a useful tool for validation of climate models.

Appendix A: Detailed site description
Ten sites, including 18 boreholes, were utilized to reconstruct the ground surface temperature history of northern Ontario and Québec for the past 500 years.A detailed description of the rock type and geological unit of each site can be found in Table A1 and heat flow studies of the region (Jessop, 1968;Jessop and Lewis, 1978;Lévy et al., 2010;Jaupart et al., 2014).
Four sites (Musselwhite, Thierry Mine, Otoskwin, Noront) comprising nine boreholes are located in northern Ontario.The Musselwhite (0601) site is located in a clearing of ∼ 60 m in diameter with a lake ∼ 330 m to the west and ∼ 180 m to the east.Three boreholes (0605, 0606, 0608) are located at Thierry Mine, ∼ 500 m away from a large clearing.This clearing is associated with the development of the nearby mine in 1934-1950.Furthermore, all three sites are ∼ 300 m from a lake.0605 and 0608 are found in a clearing due to drilling of ∼ 80 m in diameter.The Otoskwin borehole was measured in 1985 and is located ∼ 180 m from the Otoskwin River.The final four Ontario boreholes (1012,1013,1014,1015) are found in Noront.The sites are fairly flat, swampy, and muddy.They are in the McFaulds Lake project, ∼ 300 km north of the town of Geraldton and within a region referred to as the Ring of Fire of the James Bay Lowlands.N1012 and N1013 are ∼ 40 m apart, ∼ 500 from N1015 and ∼ 400 m from N1014.The complex lithological column has led to noisy temperature-depth profiles (Jaupart et al., 2014).Northern Québec is home to the remaining nine boreholes, found in six distinct sites.Nielsen Island is the northern most site of this study.It was logged in 1977 and is located on an island in the Hudson Bay.The La Grande borehole ( 0405) is in a large clearing for drilling, in fairly flat and swampy region.It is located ∼ 400 m from a main Hydro Québec power line and power line clearing.Eleonore ( 0502) is ∼ 200 m from the Opinaco reservoir and is dipping towards the reservoir.Two boreholes (0803, 0804) are in Eastmain, ∼ 220 km from the mining camp of Matagami.0803 is ∼ 50 m south of a road but is dipping away from the road.Corvet ( 0716) is located in a fairly flat area.The final three boreholes (0712, 0713, 0714) are in Camp Coulon.0713 is ∼ 20 km north of the hydroelectric station of Laforge-Deux.0712 and 0714 are ∼ 450 m apart and located on top of a small relief, with 0712 found in a clearing of ∼ 30 m diameter.

Figure 1 .
Figure 1.Map of Ontario and western Québec showing the location of sites (red dots).For sites with several boreholes (Camp Coulon, Eastmain, Thierry Mine, and Noront), the number of profiles available is enclosed in parenthesis.Black diamonds show the locations of sites that were discarded.

Figure 2 .
Figure 2. Temperature anomalies for the northern Ontario boreholes.Holes TM0605, TM0606, and TM0608 are from the Thierry Mine site; holes N1012, N1013, N1014, and N1015 belong to the Noront site.The anomaly is obtained by subtracting the estimated steady-state geotherm obtained by the least-square fit of a straight line to the bottom 100 m of the borehole temperature-depth profile.The black line represents the best linear fit, while the pink and blue lines are the upper and lower bounds, respectively, of the 2σ confidence intervals.For N1014, TM0606, and Otoskwin, the upper and lower bounds of the confidence interval are not visible due to the temperature scale.The temperature anomaly at Musselwhite was cut at 650 m.

Figure 3 .
Figure 3. Temperature anomalies for the northern Québec boreholes.CC0712, CC0713, and CC0714 are the boreholes from Camp Coulon; Ea0803 and Ea0804 are the boreholes from Eastmain.The anomaly is obtained by subtracting the estimated steady-state geotherm obtained by the least-square fit of a straight line to the bottom 100 m of the borehole temperature-depth profile.The black line represents the best linear fit, while the pink and blue lines represent the upper and lower bounds, respectively, of the 2σ confidence intervals.For N1014, TM0606, and Otoskwin, the upper and lower bounds of the confidence interval are not visible due to the temperature scale.The temperature anomaly at Nielsen Island was cut at 600 m.

Figure 5 .
Figure 5. GST histories for the northern Québec sites.Simultaneous inversion was used for Eastmain, which includes two holes.The pink and blue lines represent the inversions of the upper and lower bounds of the anomaly.

Figure 6 .
Figure 6.GST histories for the northern Québec sites.Simultaneous inversion was used for Camp Coulon, which includes more than one hole.The pink and blue lines represent the inversions of the upper and lower bounds of the anomaly.

Table 2 .
Location and technical information concerning boreholes not suitable for this study, where T o is the reference surface temperature and Q o is the reference heat flux

Table 3 .
Summary of GST history results, where T o is the reference surface temperature, Q o is the reference heat flux, and T is the difference between the maximal temperature and the reference temperature 500 years before logging.

Table A1 .
Geological unit and rock type concerning the boreholes used in this study.