North American regional climate reconstruction from ground surface temperature histories

Within the framework of the PAGES NAm2k project, 510 North American borehole temperature–depth profiles were analyzed to infer recent climate changes. To facilitate comparisons and to study the same time period, the profiles were truncated at 300 m. Ground surface temperature histories for the last 500 years were obtained for a model describing temperature changes at the surface for several climate-differentiated regions in North America. The evaluation of the model is done by inversion of temperature perturbations using singular value decomposition and its solutions are assessed using a Monte Carlo approach. The results within 95 % confidence interval suggest a warming between 1.0 and 2.5 K during the last two centuries. A regional analysis, composed of mean temperature changes over the last 500 years and geographical maps of ground surface temperatures, show that all regions experienced warming, but this warming is not spatially uniform and is more marked in northern regions.


Introduction
The energy imbalance between incoming and outgoing radiation in the upper atmosphere due to increased concentrations of greenhouse gases is well documented (e.g., Hansen et al., 2011;von Schuckmann et al., 2016).The redistribution of the excess energy between climate subsystems, the atmosphere, the oceans, and the solid Earth drives changes in globaland regional-scale climate.As the consequences of climate change are expected to be negative for natural ecosystems and society, it is necessary that the projected changes in climate be established with sufficient details and certainty to provide the framework for policy directives intended to mitigate, adapt, and build resilience at the community scale.Although there are multiple measures of climate change, surface air temperature (SAT) is the most common indicator because of the availability of data over the postindustrial period and also because it represents, in one way or another, the thermal conditions near the ground where people live.
The great majority of information on the future character and dynamics of the climate system comes from experiments with general circulation models (GCMs).GCMs are useful tools to assess future climate scenarios under different representative concentration pathways (RCPs).However, because of the limited resolution of GCMs, many climatically relevant processes operating at less than the GCM gridsize scale are parametrized differently among model teams, such that GCM simulations for the same RCP yield a climate state with a wide range of variability.Thus, GCM simulations must be compared with data to assess the validity of their climate change projections (PAGES 2k-PMIP3 group, 2015;Smith et al., 2015).
However, these proxy indicators are responses to a complex dynamical system and do not represent a direct measure of climate variability.While they allow for the determination and comparison of past climate trends, each of these methods of paleoclimatic reconstruction has different resolution, advantages, disadvantages, and uncertainties.
Furthermore, due to spatial and natural limitations, the significance of the global and regional climate reconstructions decreases as it extends back in time.Calibration disparities and different reconstruction methods among these proxies give rise to a diverse range of weaknesses and strengths, making each paleo-indicator better suitable for a specific timespan.From a large set of natural phenomena, those sensitive to temperature variations can be used as climate indicators to reproduce past temperature histories.
Collaborative efforts have been conducted under the "2k Network" of the Past Global Changes (PAGES) project to produce a global array of regional climate reconstructions for the past 2000 years using proxy datasets derived from different natural sources (PAGES 2k Consortium, 2013).It is within this multidisciplinary framework that geothermal data measured in boreholes can contribute with low-frequency trends retrieved from anomalies of the underground thermal regime.
Temperature-depth profiles measured in boreholes have commonly been used to study the magnitude and spatial variability of the flow of heat from the interior of the Earth (Bullard, 1939;Benfield, 1939;Jaupart and Mareschal, 2015, and references therein).It has been known since the times of Fourier and Kelvin that underground temperatures are affected by past surface conditions.Assuming a coupling between ground surface temperate (GST) and SAT, borehole temperature reconstructions can be used as climate indicators for hundreds to thousands of years before present.Lane (1923) and Hotchkiss and Ingersoll (1934) were the first to use temperature-depth profiles for paleoclimatic studies in an attempt to determine the timing of the last glacial retreat.It was only in the 1970s that studies to infer past climate from borehole temperature profiles (BTPs) became more systematic, developing into the field of borehole climatology (Cermak, 1971;Sass et al., 1971;Beck, 1977).
Following the work of Lachenbruch and Marshall (1986), and because of concern about climate change, paleoclimatic reconstructions from borehole temperature data have become widespread and have yielded local, regional, and global analyses (see Lewis, 1992;Bodri and Cermak, 2007;González-Rouco et al., 2009).However, the majority of the data are from the Northern Hemisphere.
In North America, several local and regional analyses have been performed (e.g., Beltrami and Mareschal, 1992;Guillou-Frottier et al., 1998;Chouinard et al., 2007).However, very few studies so far have addressed the entire North American continent.
In this paper, and within the framework of the PAGES NAm2k project, we aim to estimate regional trends in the GST change of the past 500 years in North America from a dataset containing almost twice the number of data and larger depth range ( > 300 m) compared to previous analyses.The dataset analyzed here contains 510 borehole temperaturedepth profiles distributed over the North American continent.

Methodology
The thermal regime of Earth's subsurface is governed by the outflow of heat from the Earth's interior and by the temporal variations in the GST.For a homogeneous subsurface with no internal heat sources and with no GST variations, the temperature in the subsurface increases linearly with depth.This profile can be considered as in a quasi-steady state relative to the timescale of recent climatic variations, since it depends solely on heat flux from Earth's interior, which varies over much longer timescales.Persistent temporal changes in GST propagate into the subsurface and are recorded as transient perturbations to this geothermal quasi-steady state.Because of heat diffusion, the amplitude of the subsurface anomalies is proportional to the duration and magnitude of the GST perturbations and decreases with time since their occurrence.Since these temperature fluctuations diffuse downward, only the low-frequency climate signals are preserved.To reconstruct the temporal evolution of the GSTs, the variation in the subsurface temperature as a function of depth is measured in boreholes following the procedure described in Sect.2.4.The transient perturbation is then retrieved from the borehole temperature profile (BTP) and inverted as described in Sect.2.3 in order to reconstruct the temporal GST changes.
Furthermore, borehole climatology assumes that the GST changes track long-term variations in SAT.That is, it is assumed that GST and SAT are coupled.This coupling has been confirmed by model simulations (e.g., González-Rouco et al., 2006;García-García et al., 2016), as well as data from continuous monitoring of air and ground temperature variations (Putnam and Chapman, 1996), and by comparing BTPs with meteorological records at nearby stations (Harris and Chapman, 1998).However, the relationship between SAT and GST can also be altered by transient effects in the surface conditions such as land use and associated hydrological, snow, and vegetation cover changes (Lewis and Wang, 1998;Gosselin and Mareschal, 2003a;Bartlett et al., 2004).Thus, changes in GST are not necessarily related to climate.Some of these perturbations of the surface environment can be observed at the time of measurement and should be considered prior to interpretation.When all non-climatic effects have been ruled out, the interpretation of the perturbations of the temperature profiles allows us to reconstruct the past temperature changes at the surface.

Temperature-depth equation
In order to interpret the temperature-depth profiles, we must be able to describe quantitatively the thermal regime of subsurface and also how it is affected by changes in surface temperature.This requires the solution of the heat diffusion equation for a continuous medium given by (Carslaw and Jaeger, 1959 where ρ is the density, c p is the specific heat of the medium at constant pressure, λ is the thermal conductivity, ∇ is the vector differential operator and Qs is the heat production rate per unit volume.Because heat production rates in crustal rocks are small (on the order of 1 µW m −3 ) and the effect of heat production is negligible for holes that are only a few hundred meters deep (< 1 mW m −2 ), we have neglected heat production in this study.
Assuming that heat production can be neglected ( Qs ≈ 0), that there is no advection of heat (v • ∇T = 0) and that Earth is interpreted as a homogeneous half-space, the temperature at a depth z is given by the superposition of the steady-state profile and the transient perturbation due to time variations in surface temperature: where T 0 is the long-term surface temperature, q 0 is the quasi-steady-state heat flux and R(z) is the thermal depth defined as (Bullard, 1939) where λ(z ) is the thermal conductivity at depth z .For constant conductivity, Eq. ( 2) is written as where 0 = q 0 /λ is the quasi-steady-state temperature gradient.
If thermal conductivity can be assumed constant for the measured depth interval (λ(z) = λ), the transient component of temperature is calculated from the one-dimensional heat conduction equation (Carslaw and Jaeger, 1959).
where κ = λ ρc p is the thermal diffusivity, also assumed constant for all cases (κ ≈ 10 −6 m 2 s −1 or κ ≈ 31.6 m 2 yr −1 ).The main reason to use an average value is because thermal diffusivity measurements were not made on rock samples for most of the boreholes.Equation ( 5) must be solved with initial and boundary conditions: the temperature perturbation at the surface, T (z = 0, t) = T 0 (t), no perturbation for z → ∞, T (z = ∞, t) = 0, and T (z, t = 0) = 0.The use of the onedimensional Eq. ( 5) is valid if the surface temperature variations have much larger spatial scale than their penetration depth (Clauser and Mareschal, 1995).Equation ( 5) also shows that the diffusivity determines the scaling relationship between time τ and depth L, scaling as τ ∝ L 2 /κ.Periodic surface temperature variations propagate as a damped wave with skin depth δ = √ κT /π (Jaupart and Mareschal, 2011).For standard values of κ for rocks, the amplitude of the wave associated with the annual temperature cycle is 10 % of its surface value at 10 m depth.For 100-and 1000-year cycles, the amplitude of the wave is 10 % its surface value at 100 and 300 m, respectively.

Parametrization of the temperature anomaly
Assuming that Earth's underground thermal regime is at equilibrium and there are negligible diffusivity (κ) changes in the subsurface, the transient perturbation temperature T t (z) = T (z, t = 0) defined over a semi-infinite half-space with surface temperature T (z = 0, t) = T 0 (t) at time t before present is given by (Carslaw and Jaeger, 1959) For an instantaneous temperature change T at time t before present, integrating the Eq. ( 6) yields (Carslaw and Jaeger, 1959) where erfc is the complementary error function: In order to approximate GST changes, we assume that GST can be replaced by its average value over time intervals of several years, so that the daily, annual, and solar activity cycles are removed.
Defining the ground temperature changes as T k during K time steps (i.e., T k for t k−1 < t < t k where k = 1, . .., K), we find that the transient perturbation is the sum of the contributions for each time step: F. Jaume-Santero et al.: North American regional climate reconstructions Equation ( 9) gives the temperature anomaly T t (z) due to a sequence of GST changes T k for K time intervals.The problem consists in determining the GST history from the temperature versus depth anomaly, T t (z), at a given site.This is routinely done using inversion techniques.

Inversion
Combination of Eqs. ( 2) and ( 9) yields a linear equation with the parameters T 0 , 0 , and T k for each depth with temperature data.Thus, the inversion consists of solving the resulting system of linear equations.Obtaining the solution, however, is never straightforward because the system is "illconditioned", i.e., its solution is unstable (a small change in the data causes a very large change in the solution) and, for all practical purposes, the solution is non-unique.Different methods have been developed to solve inverse problems: the Backus-Gilbert method (Parker, 1977(Parker, , 1994)), singular value decomposition (SVD) (Lanczos, 1961;Jackson, 1972), Bayesian inversion (Tarantola and Valette, 1982), Tikhonov regularization (Tikhonov and Arsenin, 1977), and Monte Carlo simulations (Mosegaard and Tarantola, 1995).One of the first applications of inversion to borehole temperature data was based on the Backus-Gilbert method (Vasseur et al., 1983); Shen and Beck (1991) proposed an algorithm based on the Bayesian approach, while Mareschal and Beltrami (1992) used SVD.Because of the very small number of parameters, these methods of inversion are not computationally intensive.The Monte Carlo method, which has been used by Mareschal et al. (1999) and Kukkonen and Jõeleht (2003), explores the entire parameter space and requires larger computational resources than the other methods.In this study, we have used SVD to find the GST history because of its simplicity and then used a Monte Carlo procedure to determine the range of model parameters that satisfy the data within some error bounds.

Subsurface temperature anomaly
In this study we determined the long-term surface temperature and quasi-steady-state geothermal gradient by linear regression to the lowermost 100 m of the measured temperature profile.This linear regression represents the geothermal quasi-steady state (Eq.2) from which the subsurface temperature anomalies are estimated.The anomaly T t (z) is obtained by subtracting this quasi-equilibrium thermal profile from the measured temperature profile.The least-squares regression also yields an estimate of the maximum error on slope and intercept estimates (95 % confidence interval).These error bounds represent the upper and lower limits for the quasisteady-state temperature profile, hereafter referred to as the extremal geothermal steady states.Figure 1 shows an example of a measured temperature profile and its estimate subsurface temperature anomaly, near Lynn Lake, Manitoba.

Singular value decomposition
After removal of the quasi-steady-state component of the temperature profile, we are left with a system of linear equations between J temperature anomalies T t (z j ) = T j for each depth and the K parameters of the surface temperature history T k : where the A j k are given by Eq. ( 9) The number of equations J could be greater, equal, or less than the number of parameters K.In general, this number is larger than the number of parameters, but this does not ensure that the system of Eq. ( 10) has a unique solution.
Written formally, the matrix of Eq. ( 10) where is the data vector, A is the rectangular (J × K) matrix containing the coefficients of the equations, and x is the vector of unknown coefficients.SVD decomposes the matrix as (Lanczos, 1961) where U is an (J × J ) orthonormal matrix in data space, V is an (K × K) orthonormal matrix in parameter space and is a J × K rectangular matrix with only non-zero values, called "singular values" λ l (l = 1, . ..L) on the diagonal, with L ≤ min(J, K).The singular values are the square root of the eigenvalues of the J × J symmetric matrix (A A).If L < J , the system is overdetermined, and if L < K, it is underdetermined.Regardless of whether the system is overdetermined, underdetermined, or both, it admits a generalized solution given by where −1 is a K × J rectangular matrix with L elements 1 λ l on the diagonal completed with zeros.This provides a solution which is usually not very meaningful (Mareschal and Beltrami, 1992) because it is unstable and dominated by noise.The instability of the solution comes the presence of very small singular values λ l .In the case of borehole temperature profiles, the fifth largest singular value is 0.01 times the largest one, and the tenth is < 10 −8 times the largest one, that is, less than numerical noise.In order to stabilize the solution, we eliminate the part associated to the very small singular values.This is done by replacing with 0 the inverse of all the singular values less than a "cut-off value", typically on the order of 10 −2 .This means that the solution is obtained as a linear combination of four orthogonal vectors in parameter space.Each vector represents a surface temperature history, and the vectors selected are those that have the largest impact on the data.By eliminating the small singular values, we choose to neglect the part of the solution that has little or no effect on the data and therefore cannot be determined.In general, the selection of a cutoff value is done by trial and error, by increasing the number of singular values and inspecting the solution for signs of instabilities and loss of resolution, i.e., large non-physically meaningful fluctuations or no useful information.For this study, we used a cut-off of 0.03, which resulted in four singular values being retained for all profiles except for CU-C-357 measured in Cuba, where only 3 singular values were retained.
The choice of a proper parametrization is useful to reduce the number of parameters to be estimated.This can be achieved by increasing the duration of the GST history model time intervals.For very long reconstructions a logarithmic distribution has been used (e.g., Mareschal et al., 1999).For the present study, we have used a model consisting of a series of 10 time intervals of varying duration after testing with different parametrizations and verified that similar results were obtained (see Appendix A).Their temporal length is smaller for the near (past 100 years) than for the remote past.The distribution used here is , 25, 50, 75, 100, 150, 200, 250, 300, 400, 500}. (15) When regional averages are made, the GST histories are shifted in time to account for the date when they were logged (i.e., years before present is the year of measurement).
As an example, Fig. 2 shows the result of inversion of the subsurface temperature anomaly for the Fox Mine site, and the results from the inversions of the two extremal geothermal steady states.

Forward model
GST histories can be forward-modeled using Eq. ( 9) to assess the fit of the SVD inversion with respect the initial anomaly profile.A Monte Carlo procedure was applied (Mareschal et al., 1999;Kukkonen and Jõeleht, 2003;Chouinard et al., 2007) by randomly perturbing the model parameters to find the range of GST histories that fit the data within a maximum root mean square (RMS) error less or equal than the difference between the forward-modeled SVD reconstruction and the anomaly.Using the Monte Carlo approach to invert the temperature profiles is particularly inefficient because it requires a very large number of simulations to explore the entire parameter space.It requires at least 10 7 -10 8 longer computational time than using the SVD inversion.However, this can be alleviated by using a priori information or the result of an existing GST history from inversion to reduce the region explored in parameter space.After the Monte Carlo inversion, the mean and standard deviation of all the accepted models are estimated to show the trend of all the solutions with a same or better fit than the inversion for four singular values.For the present study, we halted the calculations after 500 models are accepted or after 5 million forward model comparisons.
This is illustrated in Fig. 3, which shows the results of the Monte Carlo inversion for the Fox Mine temperature profile.

Data
We have compiled from different sources (Table 1) a set of temperature-depth profiles for North America.Thousands of borehole temperature profiles have been measured in North America, but the majority of them are not suitable for climate reconstructions.For instance, bottom hole temperatures, commonly measured during oil exploration drilling, are not measured at equilibrium and are affected by errors several times larger than the signals we want to detect.Water wells are usually too shallow to be useful and likely to be affected by water flow.Many holes were drilled for geothermal energy in the western US, but these are often perturbed by water circulation.For heat flow or climate studies, the most useful boreholes are those that have been drilled by mining companies for exploration or development purposes.Oil exploration wells cannot be used for several reasons: holes that are not put in production must be cemented and they are not accessible for steady-state measurements.In addition, oil exploration boreholes have a large diameter and are susceptible to perturbations due to convection in the hole.Furthermore, sedimentary rocks are permeable and often affected by convection as well.Hence, their temperature profiles are not suitable for climate studies.Drilling perturbs the thermal regime of the subsurface around the drill site and some time is needed for thermal re-equilibration.As a rule of thumb, the time to return to equilibrium is ∼ 5-6 times the duration of drilling.The temperature in the hole is measured with a calibrated thermistor.The probe is lowered in the hole and measurements are made at fixed intervals along the length of the hole, which results in varying depth intervals as most boreholes are inclined.The sampling interval is usually 10 m, and sometimes 50 ft for US and old Canadian temperature logs.Continuous measurements can be obtained, but these are not common because they require heavy equipment.Measurements made above the water table are rarely equilibrated; consequently, the upper 20 or 30 m of the temperature logs must be discarded.This is also done in order to eliminate the annual temperature variation signal.In heat flow studies, core samples must be obtained to determine the underlying rock's thermal conductivity and heat production.Changes in thermal conductivity are thus included in the interpretation of these data.

Data selection
Different criteria have been applied in selecting the temperature profiles.Temperature profiles must be at least 300 m deep to contain the signal to allow for the reconstruction of the climate of the past 500 years.Profiles must include at least 10 measurements, as well as measurements in the uppermost 100 m.Profiles that meet these conditions are then visually inspected to detect discontinuities, signs of water flow, or other perturbations that make them unsuitable for interpretation.The vertical temperature gradient profile amplifies the noise and usually provides a better diagnostic for the level of noise in the measurements.Although we have not established a quantitative criterion for selecting profiles based on the noise level, we have examined the vertical gradients to eliminate obviously unsuitable profiles.
After selection process, we retained 510 profiles.These data will be available in a public database in Figshare (Jaume-Santero et al., 2016).Borehole locations are not uniformly distributed across the continent (Fig. 4).Several regions are very poorly sampled because they are very difficult to access (Alaska and most of Canada, north of 56 • ).Furthermore, in the northernmost regions, drill holes cannot be routinely logged because of permafrost.Temperature logging in frozen ground requires special equipment to be emplaced at the end of drilling and is very costly.The southern part of the Canadian Shield is the region most extensively sampled because of the mining activity and because the temperature profiles are less likely to be perturbed in the crystalline rocks of the Shield.In contrast, numerous drill holes are available in the southwestern US, but most of them cannot be used because they are perturbed by water flow.The sedimentary cover in many regions of the US explains why no suitable holes have been found for many states, including Texas and Oklahoma and the southeastern US.This very uneven distribution of suitable boreholes is demonstrated in Table 2, which shows the number of temperature profiles for each of the regions defined for Pages2k (McKay, 2014).

Results and discussion
All 510 borehole temperature-depth profiles were inverted individually to reconstruct the GST histories for the past 500 years.The model consisted of a series of 10 temperature change intervals of varying temporal duration following the distribution (Eq.15).For the inversion, we used the SVD inversion with a cutoff of 0.03, retaining four singular values.We also used the Monte Carlo methodology to estimate the range of parameter values consistent with the data.The means of the GSTs obtained by Monte Carlo are similar to the solution by SVD inversion.With the condition that the   RMS difference between model and data be no larger than the misfit for the SVD, the 2σ range of accepted models is no larger than 0.44 K.

North American GST change
We have calculated the variation in GST for North America by averaging all the Monte Carlo inversions.The averaging was done on a yearly basis because the logging dates vary between boreholes from 1958 to 2014 (Fig. 5).
Figure 5 shows the individual Monte Carlo inversions together with their average.We believe our results are consistent because similar mean North American GST histories were obtained from different parametrizations (see Appendix A).However, individual inversions in Fig. 5 exhibit a wide variability due to the large range of latitudes (∼ 80 to ∼ 18 • N) in the dataset of GST reconstructions.
Nevertheless, a clear warming transition is observed from the preindustrial era (1500-1800) to the postindustrial era .The temperature difference between the preindustrial mean (1500-1700) and the mean between the years ) is 1.1 K.Because of the marked warming of the past 50 years, the total change of the average GST is 1.8 K between preindustrial time and the year 2000.
These results agree with findings of other GST reconstructions (Huang et al., 2000;Harris and Chapman, 2001;Beltrami and Bourlon, 2004;Pollack and Smerdon, 2004).Furthermore, they agree with instrumental data, CRUTEM4 (Jones et al., 2012;Morice et al., 2012), and pollen and treering reconstructions (PAGES 2k Consortium, 2013;Trouet et al., 2013).All the reconstructions are presented as departures from the 1904-1980 temperature mean (Fig. 6).However, the reconstructed GST warming signal for the past 200 years is greater than results from pollen reconstructions, being consistent with findings of PAGES 2k-PMIP3 group (2015).Furthermore, multi-centennial temperature reconstructions for North America and the Northern Hemisphere, based on multiproxy records, showed trends similar to temperature-depth reconstructions: an unclear coldwarm trend followed by a clear increase in temperature for the past two centuries (Moberg et al., 2005;Mann et al., 2008;PAGES 2k-PMIP3 group, 2015).This warming has also been recorded by instrumental data for the last century (Hansen et al., 2010).However, the difference between the long-term preindustrial temperature mean and the recent past trend is larger in the GST histories than in the pollenbased and tree-ring reconstructions.These disparities among different proxy-based reconstructions can be attributed to a combination of factors as discussed in Pollack and Smerdon Tree rings (Trouet et. al., 2013) CRUTEM4 NAmGSTH (300 m) Pollen (Trouet et. al., 2013) Pollen (Viau et. al., 2006) Pollen (Viau et. al., 2012) Figure 6.Mean North American GST history (blue) and maximum temperature range of accepted models (∼ 0.44 K) obtained from the Monte Carlo method (blue shade).Also shown are proxy-based SAT reconstructions for North America from 1500 to 2000 CE.All anomalies are displayed as departures from the 1904-1980 mean.
(2004).For instance, while a significant part of boreholes are located in higher latitudes (eastern and central Canada), treering data are mainly obtained in lower latitudes (western US).Therefore, the spatial distribution of proxies could explain colder temperatures.Other possible reasons for those differences are the seasonal bias of the proxies and the limitation of borehole climatology in resolving short-term variability.The Little Ice Age (LIA) is not resolved because the boreholes were truncated at 300 m, which is too shallow to allow for a clear LIA signal in most of the borehole profiles as can be shown with synthetic models (Mareschal and Beltrami, 1992) and was confirmed in several studies (Guillou-Frottier et al., 1998;Chouinard et al., 2007;Pickler et al., 2016b).Some profiles, such as the Fox Mine shown in Fig. 2, may indeed show the LIA cooling, but the majority of them do not.In addition, because the LIA signal may vary both in time and in amplitude between regions, a marked signal cannot be expected from averaging weak and inconsistent signals.

Regional averages
The PAGES NAm2k working group divided the North American continent into seven subregions for paleoclimate studies (McKay, 2014).The distribution of boreholes between these regions is extremely uneven as shown in Table 2, with only four regions appearing adequately sampled (central and eastern Canada, midwestern US, Arctic, and Pacific northwest).Furthermore, the sampling in the Arctic and the Pacific northwest is very biased because all the boreholes are close to the coast (Fig. 4).For the three other regions, the sampling is insufficient to obtain robust climate trends.
A warming by ∼ 1.8 K for the past 200 years is observed in the Arctic (Fig. 7a), but the histories show wide variability.This variability suggests the need for smaller-scale regional analysis such as the pollen-based reconstructions of Gajewski (2015) and Viau and Gajewski (2009).Their findings illustrate that recent Arctic increases in temperature have exceeded natural climate variability, which is consistent with borehole GST reconstructions.
The region of the Pacific northwest (western Canada and northwestern US) shows an increase in temperature of ∼ 0.8 K with a 95 % variability range of ∼ 3.4 K for the last two centuries (Fig. 7b).This warming is consistent with previous findings (Majorowicz and Safanda, 2001).
An average warming of ∼ 1.1 K with a 95 % variability range of ∼ 2.2 K for the past two centuries is observed for central and eastern Canada (Fig. 7c), agreeing with previous studies (Beltrami et al., 1992;Guillou-Frottier et al., 1998).
The western US GST mean shows a small increase in temperature of ∼ 0.2 K ± 1.8 K (Fig. 7d).This could be the result of strong irrigation processes and water flow at the sampling locations, but the number of borehole temperature profiles available in the region are insufficient to verify this.The limited number of useful borehole temperature profiles for the western US (only nine) was logged in the 1960s, the most recent of which was measured in 1970.Thus, it is not possible to reconstruct the past 40 years, when the increase in temperature recorded in weather stations was more marked.
The average reconstruction for the Midwestern US suggests a warming of ∼ 1.3 K ± 2.0 K for the last 50-year average (Fig. 7f).This recent warming has also been observed in previous GST reconstructions as well as SAT records (Skinner and Majorowicz, 1999) and could reflect the significant land use change in the region.
A warming of ∼ 1.0 K ± 1.0 K has been reconstructed for the last 200 years in the eastern United States (Fig. 7e).However, due to the rejection of borehole profiles affected by topography and water flow, the number of reconstructions made is too small to describe climate trends of the region with confidence.
There is a warming trend of ∼ 3.0 K ± 3.6 K until the mid-1960s in the Caribbean (Fig. 7g).Due to the low number of profiles sampled in Mexico (0) and the Caribbean (4), it is not possible to obtain a robust reconstruction for this region.

Geographical representation
A North American regional analysis of GST changes is presented as six geographical maps for different 50-year time intervals during the last 300 years (Fig. 8).
Trends prior to 1681 are not shown because they did not yield significant information.However, a small (∼ 0.5 K) cooling is observed in certain regions.Previous small-scale regional analyses have reconstructed a LIA signal during this period (e.g., Beltrami and Mareschal, 1992;Chouinard et al., 2007).Furthermore, the regional variability of the cooling is consistent with previous studies, illustrating that not all regions of North America present a LIA signal (Gosselin and Mareschal, 2003b;Mann et al., 2009).However, due to the truncation at 300 m of the temperature-depth profiles analyzed here, a clear LIA signal cannot be resolved.
Figure 8 indicates a warming trend of ∼ 1-2 K in most parts of North America during the last 200 years.This is consistent with previous studies (Huang et al., 2000;Harris and Chapman, 2001;Beltrami et al., 2003).A cooling trend is observed in central California.Stevens et al. (2008) show how this differs from the output of the ECHO-G model and postulates that it is the result of intensive irrigation in California's central valley, which could drive a regional cooling signal (Kueppers et al., 2007).A similar cooling signal is observed in British Columbia which might be associated with irrigation in the Fraser Valley.On the Canadian east coast, Newfoundland presents little to no change with respect to the long-term mean.This agrees with meteorological data for the region (Gullett and Skinner, 1992).The absence of temperature profiles along the Gulf coast and Mexico does not allow for any determination of climate trends.The southwestern US is also a region where the number of boreholes is not enough for reliable reconstructions.For these regions, a multi-proxy approach would be necessary to improve the reconstruction of regional past climate in regions with an insufficient number of borehole profiles.

Conclusions
The average North American GST change reconstructed from 510 boreholes deeper than 300 m suggests a warming of ∼ 1.8 K for the last 200 years.However, these temperatures exhibit a wide range of spatial variability among all regions.For instance, reconstructed regional GST changes for seven climate distinct regions, defined within the PAGES NAm2k project, suggest a warming range of ∼ 0.5 to ∼ 2.0 K with a standard deviation (?) no smaller than 0.5 K. Furthermore, regional variations in GST yield a warming range of 1 to 2 K between 1780 and 1980.These warming trends are consistent with multi-proxy reconstructions.
Although the number of borehole temperature profiles for North America has been notably increased in our study, it is still insufficient to guarantee a non-spatially biased regional analysis because their distribution is not sufficiently uniform.Nevertheless, despite spatial and natural limitations, subsurface thermal profiles obtained from boreholes provide robust long-term GST histories which could be used to improve climate multi-proxy-based reconstructions.Those enhanced reconstructions would bring out worthwhile information for a straightforward assessment of past climate GCM outputs.

Data availability
The sources of all the data used in this study are listed in Table 1.
North American borehole temperature profiles valid for climate reconstructions were uploaded to Figshare (https:// figshare.com/s/0a1d213c3814024c4333)and published with doi:10.6084/m9.figshare.2062140.These profiles present a set of depths (z) with their associated temperature T (z).Attached to the profiles there are important supplementary metadata such as coordinates, logging date, or the person who measured the profile.The dataset is presented in two formats, as comma-separated value (CSV) files and as TABULAR text files.Example codes to load the different elements are also included in the CODE folder.

Figure 1 .
Figure 1.Temperature profile measured at Fox Mine (CA-9519), Lynn Lake, northern Manitoba, Canada.Main panel: measurements are shown in circles T (z), the red line represents the geothermal steady state, obtained by linear regression of the lowermost 100 m, and extrapolated to the surface (z = 0).Blue and green lines represent the 95 % confidence interval from the linear regression.Inset: transient perturbation or anomaly relative to the geothermal steady state (red line) and the 95 % confidence interval (blue and green lines).For this site, the geothermal steady state is given by o z+T 0 = (10.51K km ±0.34 K km )×z+ (1.44 • C ± 0.19 • C) (z in km).

Figure 2 .
Figure 2. Ground surface temperature history for CA-9519 (Fox Mine, 1995).The red line represents the GST history reconstructed from inversion.The blue and green lines are the GSTs for the anomalies estimated from the 95 % uncertainty limits of the quasisteady-state profile.

Figure 3 .
Figure3.CA-9519 (Fox Mine, 1995)  mean GST history (red) and 2σ uncertainty intervals (blue) from the Monte Carlo inversion.The grey lines represent all the perturbed models within an interval determined by the RMS misfit from the SVD inversion.

Figure 4 .
Figure 4. Location of the 510 selected boreholes.The colors represent the maximum depth of each borehole.

Figure 5 .
Figure 5. Mean North American GST change (black).Shown in blue are the 510 GST reconstructions inferred from the Monte Carlo inversion.

Figure 7 .
Figure 7. Mean GST histories (black), the blue shaded areas represent the 95 % confidence interval associated with the climate variability of each area.Regional mean temperatures are shown until the year of measurement of the most recent thermal profile in each region.(a) Arctic (78 sites), (b) Pacific northwest (78 sites), (c) central and eastern Canada (220 sites), (d) western US (21 sites), (e) eastern US (9 sites), (f) midwestern US (100 sites), (g) Caribbean (4 sites).

Figure 8 .
Figure 8. Spatial variability in the GST variation (in kelvin) from 1681 to 1980.Each panel shows a regionally interpolated mean GST over 50 years.The surface has been masked for zones without at least one datum within a radius of 400 km.Ground surface temperature changes are presented as departures from long-term mean surface temperatures prior to 1500 CE.

Table 1 .
Sources of the temperature-depth profiles.