Numerical Studies on the Impact of the Last Glacial Cycle on Recent Borehole Temperature Profiles: Implications for Terrestrial Energy Balance

Reconstructions of past climatic changes from borehole temperature profiles are important independent estimates of temperature histories over the last millennium. There remain, however, multiple uncertainties in the interpretation of these data as climatic indicators and as estimates of the changes in the heat content of the continental subsurface due to long-term climatic change. One of these uncertainties is associated with the often ignored impact of the last glacial cycle (LGC) on the subsurface energy content , and on the estimate of the background quasi steady-state signal associated with the diffusion of accretionary energy from the Earth's interior. Here, we provide the first estimate of the impact of the development of the Laurentide ice sheet on the estimates of energy and temperature reconstructions from measurements of terrestrial borehole temperatures in North America. We use basal temperature values from the data-calibrated Memorial University of Newfound-land glacial systems model (MUN-GSM) to quantify the extent of the perturbation to estimated steady-state temperature profiles, and to derive spatial maps of the expected impacts on measured profiles over North America. Furthermore, we present quantitative estimates of the potential effects of temperature changes during the last glacial cycle on the borehole reconstructions over the last millennium for North Amer-ica. The range of these possible impacts is estimated using synthetic basal temperatures for a period covering 120 ka to the present day that include the basal temperature history uncertainties from an ensemble of results from the calibrated numerical model. For all the locations, we find that within the depth ranges that are typical for available boreholes (≈ 600 m), the induced perturbations to the steady-state temperature profile are on the order of 10 mW m −2 , decreasing with greater depths. Results indicate that site-specific heat content estimates over North America can differ by as much as 50 %, if the energy contribution of the last glacial cycle in those areas of North America that experienced glacia-tion is not taken into account when estimating recent subsur-face energy changes from borehole temperature data.


Introduction
The past 2000 yr are an important time period for datamodel comparisons as a means of assessing past variability and change, as well as evaluating the performance of the Published by Copernicus Publications on behalf of the European Geosciences Union.
same general circulation models (GCMs) that are used for 21st-century climate projections.(e.g., Schmidt et al., 2013;Phipps et al., 2013;Fernández-Donado et al., 2013;Coats et al., 2013).With regard to the latter, the abundance of paleoclimate data during this period provides numerous estimates of past climatic conditions against which GCMs can be tested (Jansen et al., 2007;Randall et al., 2007;González-Rouco et al., 2009).While decadal to centennial variability and change can be well characterized during the last 2000 yr by indirect methods, instrumental records rarely exceed 100-200 yr in duration.GCM performance on longer timescales is nevertheless of particular relevance to future climate projections, which require accurate representations of radiatively forced change and internal variability over decades and centuries.
Among the variously available proxy indicators, geothermal measurements from terrestrial borehole temperature profiles (BTPs) provide a unique estimate of past changes in the Earth's surface energy balance (see Pollack and Huang, 2000, Bodri and Cermak, 2007, and Gonzalez-Rouco et al., 2009, for reviews).These records have most widely been used to estimate past temperature changes, but can also be used as estimates of past energy fluxes at the land surface.These latter reconstructions are possibly less vulnerable to processes that could potentially interfere with the coupling of surface air temperatures (SATs) and ground surface temperatures (GSTs), which have been discussed extensively in the literature (e.g., Beltrami, 1996;Schmidt et al., 2001;Beltrami and Kellman, 2003;Stieglitz et al., 2003;Bartlett et al., 2004Bartlett et al., , 2005;;Smerdon et al., 2003Smerdon et al., , 2004Smerdon et al., , 2006Smerdon et al., , 2009;;Demetrescu et al., 2007).In other words, subsurface heat storage Q s is independent of any surface temperature association.Estimating Q does not require the solution of an ill-posed and assumption-laden inverse problem, as is the case for GST history inversion, but can be estimated from observed BTPs, with only minimal assumptions about the volumetric heat capacity.If conductive processes dominate, changes in the heat content of the shallow subsurface do not involve complicated nonlinearities.Previous work that has used BTPs to estimate continental heat storage over the last 500 yr has shown that the continental subsurface was second only to the oceans in terms of the total amount of heat absorbed in the second half of the 20th century (Beltrami, 2002b;Beltrami et al., 2002;Huang, 2006).
Interpretations of borehole temperatures as records of past energy fluxes at the land surface (e.g., Wang and Bras, 1999;Beltrami et al., 2000;Beltrami, 2002a) are important within the context of the quantification of the heat content (particularly its rate of change) of the Earth's major climate subsystems (ocean, atmosphere, cryosphere, and continents) from the perspective of climate system dynamics (Pielke, 2003;Hansen et al., 2005;Trenberth, 2009;Hansen et al., 2011Hansen et al., , 2013)).Numerous studies on the contributions of these distinct components have been published (Levitus et al., 2001(Levitus et al., , 2005(Levitus et al., , 2012;;Beltrami, 2002b;Bindoff et al., 2007;Davin et al., 2007;Murphy et al., 2009;Church et al., 2011;Ortega et al., 2013;Rhein et al., 2013).Moreover, the correct partitioning of heat into the various climate subsystems is important within climate model simulations that seek to project the overall energy imbalance of the planet and how energy fluxes will ultimately impact the character of each subsystem.With regard to the continental subsurface specifically, a large infusion of heat into the subsurface can impact important soil processes, including hydrology (Zhu and Liang, 2005;Bense and Kooi, 2004), biogeochemical processes such as CO 2 production via microbial and root respiration and long-term carbon storage (Risk et al., 2002(Risk et al., , 2008;;Kellman et al., 2006;Bekele et al., 2007;Bindoff et al., 2007;Diochon and Kellman, 2008), and the spatial distribution and depth of permafrost (Sushama et al., 2007;Lawrence et al., 2008).Not only are these processes important for risk assessments associated specifically with soil processes and functions in the future, they are also linked to climate feedback mechanisms and are relevant to simulations of the climate system as a whole.
Given the above-discussed importance of subsurface temperatures and heat storage for the continents, comparisons between paleoclimatic GCM simulations and borehole estimates of past GSTs or energy fluxes comprise a valuable opportunity to evaluate GCM performance.These comparisons are nevertheless dependent on the robust characterization of past changes in the subsurface energy reservoir and the associated uncertainties.Herein we specifically investigate an important but underappreciated uncertainty in the interpretation of BTPs as past temperature and energy flux estimates.Our focus is principally on the impact of the late stage of the last glacial cycle (LGC), including the Last Glacial Maximum (LGM, ≈ 26-21 ka), and the subsequent warming to the current deglacial state on the background steady state of the subsurface thermal regime.Most climatic interpretations of BTPs rely on a principal assumption: downward-propagating signals associated with energy balance changes at the land surface are imprinted on a subsurface steady-state signal associated with the outward flow of heat from the Earth's interior (changes in this outward flow occur on the order of millions of yr).This steady-state signal is approximated below a few hundred meters in BTPs by a linear function of increasing temperatures with depth, assuming that this signal represents the geothermal heat flow from below (see Sect. 2.2).The thermal effects of the LGC have been discussed previously as sources of bias in estimates of this background steady-state signal (e.g., Hotchkiss and Ingersoll, 1934;Jessop, 1971), but recent results have more quantitatively characterized the important impact of the LGC on BTP interpretations (Rath et al., 2012).Though known to be significant for a long time, LGC effects on heat flow density estimations have only recently been characterized in appropriate detail within general geothermal studies (e.g., Slagstad et al., 2009;Majorowicz and Wybraniec, 2011;Westaway and Younger, 2013).Despite important insights gained from earlier work on the LGC and its relation to BTPs (Hartmann and Rath, 2005;Jaupart and Marechal, 2011;Rath et al., 2012), results were based on highly idealized GST histories used to represent the temperature evolution from the LGM to the present day.Ice sheets are, however, inexorably linked to climate in a complex set of relationships; they are a dynamical system responding to climate change on the local scale, but also generate large-scale climate variability.Ice sheets store fresh water, are reflectors of shortwave radiation, alter surface topography, change large-scale atmospheric circulation patterns, and create a subsurface thermal environment that varies in spatial and temporal scale throughout their history (e.g., Marshall et al., 2002;Rolandone et al., 2003).Characterizing the actual impacts on observed BTPs in North America therefore requires modeling of the complexities of surface conditions originating from the last glaciation.This is far from trivial, as these surface conditions do not only include the effects of the ice sheet and its partly unknown dynamics, but also other effects related to the full glacial cycle.These include the influence of oceanic transgressions due to the timing of isotatic effects, as well as proglacial lakes, which could be relevant in stages of progressive deglaciation.

Clim
Herein we study the effect of the buildup and retreat of the Laurentide ice sheet on subsurface temperature profiles using basal temperatures derived from the Memorial University of Newfoundland (MUN) glacial systems model (MUN-GSM).Details on the procedures employed for constructing surface temperatures from the MUN-GSM and the methods relevant to modeling and inversion of the resulting synthetic BTPs are given in Sects.2.1 and 2.2, respectively.In Sect.3, we investigate the perturbation of the background steady-state signal, which in turn perturbs the corresponding heat flux profile and therefore the overall subsurface heat storage.Subsequently, in Sect.4, we quantify the effect of this perturbation on GST inversions from synthetic BTPs, focusing on the depth ranges typical for borehole measurements that have been used in global or hemispheric studies.

Theoretical framework
Energy balance changes at the surface of the Earth propagate into the subsurface and, depending on their magnitude and duration, can leave a detectable signature after periods as long as 100 kyr.This characteristic of heat propagation in the terrestrial subsurface has been exploited to estimate past changes in GST and energy fluxes from BTPs measured across all continental regions.Inversion of BTPs to derive GST histories is an operation that transforms a temperaturedepth profile T (t 0 , z) at a given time t 0 (the time of measurement) into a temperature-time profile T (t, z 0 ) at the ground surface z 0 (e.g., Beltrami and Mareschal, 1991;Pollack and Huang, 2000;Huang et al., 2000;Harris and Chapman, 2001;Beltrami, 2002a;Rath and Mottaghy, 2007).If the rock prop-erties are assumed to be vertically homogeneous, depth and time are linked by the thermal diffusivity κ.This is a bulk thermophysical property of the subsurface soil and rock defined as the ratio between thermal conductivity λ and the volumetric heat capacity C = ρc p .More information on the relevant physics and methodologies can be found in the recent reviews of the borehole paleoclimatology literature by Bodri and Cermak (2007) and González-Rouco et al. (2009).In the following section, we briefly expand on the theoretical basics of BTP modeling and interpretation, as far as they are relevant to our subsequent analysis.Our first focus, however, is a more detailed description of the MUN-GSM that we employ to estimate the ground surface temperature history over North America over the LGC.

The MUN-GSM
We use the MUN-GSM to simulate the thermal history of the ground surface underneath the Laurentide ice sheet during the LGC.This simulated temperature history will then be used to forward model the subsurface propagation of the associated thermal perturbation in BTPs down to several thousand meters in depth.The MUN-GSM is a 3-D thermomechanically coupled ice sheet model that includes a permafrost-resolving thermal model at its base, and is asynchronously coupled to a viscoelastic model of the glacial isostatic adjustment process.Ice dynamics are computed under the shallow ice approximation (SIA) with fast flow due to sliding or subglacial till deformation when the basal temperature approaches the pressure melting point.The SIA, described in detail by Greve and Blatter (2009), is appropriate for long-term simulations of continental-scale ice sheets away from fast-moving ice streams.The model was originally described by Tarasov and Peltier (1999), while subsequent improvements were discussed in follow-up publications (Tarasov and Peltier, 1999, 2002, 2004, 2007;Tarasov et al., 2012).In contrast to other GSMs, it has been calibrated against a comprehensive body of available data (Tarasov et al., 2012) using a Bayesian methodology.The calibration produces a posterior distribution of higher-likelihood ensemble parameter sets given model fits to a diverse set of constraint data.This model is driven by a climate based on interpolation between a modern state from the NCAR/NCEP reanalysis (Kalnay et al., 1996) and an LGM state derived from PMIP II archived 21 ka slice GCM simulations (http: //pmip.lsce.ipsl.fr/).The interpolation uses an index function derived from Greenland ice cores (Tarasov andPeltier, 2003, 2004).The interpolation and definition of the LGC climate state are all subject to a set of calibrated ensemble parameters.
The gridded output from the MUN-GSM is available in the domain defined by the corner points (172.5 • W, 34.75 • N) and (42.5 • W, 84.75 • N), with a resolution of 1 o longitude by 0.5 • latitude, covering most of North America, at constant 1000 a time intervals from 120 ka to the present.The region of interest in the present paper lies between (172.5 • W, 34.75 • N) and (52.5 • W, 75 • N), and is shown in Fig. 1.GSTs for each grid cell are set to the basal temperature of the ice when ice covered, 0 • C when lake covered (simplest approximation for proglacial lakes receiving ice discharge), and to that computed by a temperature at the top of the permafrost (TTOP) correction to the calibrated model climate forcing when sub-aerial.The latter TTOP model-derived correction (Smith andRiseborough, 1996, 2002;Riseborough, 2002) seeks to correct for varying snow cover and vegetation, and for the differing thermal conductivity of thawed versus frozen land (Tarasov and Peltier, 2007).
There is currently no basal hydrology in the GSM used for this study.Though there is some evidence for the existence of subglacial lakes in the area of the Laurentide ice sheet (Livingstone et al., 2013), the lack of basal hydrology has probably limited relevance for the thermodynamics, given all the other uncertainties.On bedrock, it will have no impact except for some horizontal advection of heat (basal water will be associated with basal ice at the pressure melting point), which is likely negligible given the approximately 50 km grid cell scale.On sediment, one can likely assume that all subglacial sediment is saturated, as is done for the permafrost calculation in the GSM.
The uncertainty estimates for the surface temperature forcing are from a Bayesian calibration of the MUN-GSM against a large and diverse set of geophysical and geological constraints (Tarasov et al., 2012).The North American configuration of the model had 39 ensemble parameters, endeavoring to capture uncertainties in the climate forcing, ice dynamics, and ice calving.The calibration of these 39 parameters involved over 50 000 full GSM runs and Markov chain-Monte Carlo sampling of tens of millions of parameter sets using Bayesian artificial neural network emulators of the GSM.As the calibration targets did not involve direct paleoclimatic proxies, the surface temperature chronology confidence intervals are likely too narrow, especially over regions that were predominantly ice free.For the duration of the ice cover, given the constraints and represented physics, the calibrated results offer the best available reconstruction of subglacial temperatures.The 75 % confidence intervals are invoked from application of Chebyshev's inequality to the two standard deviation sample range of the ensemble.In practice, we invoke the Saw et al. (1984) extension that uses sample mean and standard deviation.The assignment of a weaker interpretation of the 2σ confidence interval avoids the assumption of an underlying Gaussian distribution of misfits required for the standard 96 % assignment.

The subsurface thermal model
Given simulated GSTs from the MUN-GSM, the subsurface propagation of the LGC thermal signal into the subsurface at specific locations can be modeled using the onedimensional thermal diffusion equation, assuming a homo- production (Wm −3 ), z is depth (m), and t is time (s).We assume the common convention that z is geneous, isotropic half space (Carslaw and Jaeger, 1959), as where κ is the thermal diffusivity of the rock defined as , and t is time (s).We assume the common convention that z is positive downwards.
A solution of Eq. ( 1) requires the specification of boundary conditions (BCs).Temperatures in the first few hundred meters of the Earth's interior are governed by the outward flow of heat associated with internal energy, and temperature perturbations propagating downward from the land-atmosphere boundary associated with changes in the surface energy balance.Because the changes associated with the outward flow of heat occur on timescales of millions of yr, we assume a constant heat flux (Neumann type) BC at infinity within the half space.
The upper BC imposed at z = 0 is of the Dirichlet type and a function of time.We parametrize this time dependence as a series of step changes in temperature at given points in time based on the output from the MUN-GCM.This BC at the surface was chosen because of both physical and practical reasons.The surface condition in the area and time period influenced by the LGC may be characterized by different regimes.At the ice base, the temperature is continuous, but not the gradient, because of the phase change and energy sinks or sources at the boundary.These result from a variety of processes.They include phase change (melting/freezing), advective heat transport, and heating related to sliding and ice deformation work, which are both source terms in MUN-GSM.Though it is possible to construct a Neumann boundary condition at the ice base, there is much better control on the temperature, as the melting temperature of the ice is well known as a function of pressure, and thus of ice thickness.In the water-covered times and areas, the temperature is again better constrained than the gradient, as the density maximum is well known for water of a different salinity, and can be used to define the BC, though this may be a matter of discussion.For the shallow proglacial lakes (as opposed to oceanic transgressions), it is set to the melting temperature for the current study.For direct contact with the atmosphere, meteorological simulations use a complex parametrization of the boundary layer physics (e.g., Stensrud, 2007).For the long time periods involved here, however, it has been shown that the GST is related to the SAT following a simple relationship (see Smerdon et al., 2003, 2006, and Stieglitz and Smerdon, 2007).
Given the above-described BCs, an analytic solution to Eq. ( 1) gives the temperature anomaly at depth z and time t, for a series of K step changes at the surface, as where T s (z) represents the linear steady-state temperature profile depending only on thermal conductivity and basal heat flow density, while T t (z, t) represents the perturbation by transient changes in the surface temperature boundary condition (Mareschal and Beltrami, 1992).The equilibrium BTP is written as Here, the third term, representing the volumetric heat production h, is often neglected for borehole paleoclimatic studies that use shallow boreholes (< 1000 m) and estimate changes on a timescale of several centuries.Equation (3) reduces to a linear relationship between temperature and depth when h is assumed zero.This linear component can be removed from BTPs, thus reducing Eq. ( 2) to only the transient component of the temperature signal associated with the downwardpropagating surface perturbation.This transient component can be expressed as The subsurface anomaly T t (z, t) described by Eq. ( 4) represents the cumulative effect of the energy balance at the surface.The magnitude of T t (z, t) is directly proportional to the subsurface cumulative heat integral over depth, and it represents the total heat absorbed or released by the ground according to Here, ρc is the volumetric heat capacity, chosen here as 3 × 10 +6 J m −3 K −1 after Čermák and Rybach (1982).
A discrete version of the transient part of Eq. (2) at N depths z i can be derived, leading to a linear matrix equation defined as where T is a vector of borehole temperatures at depths z, and T gst represents the vector of discrete temperature steps T k at times t k .A detailed derivation of this equation was lately given by Brynjarsdottir and Berliner (2011).Equation ( 5) can thus be inverted to estimate T gst given a measurement of T obs .The equation is nevertheless not easily invertible because the problem is ill posed (see Hansen, 1998, andAster et al., 2013, for a detailed introduction).In the case considered herein, the truncated singular value decomposition (SVD; see Lanczos, 1961, Varah, 1973, and Lawson and Hanson, 1974) is used to perform the inversion.This approach was introduced in the context of GST estimation by Beltrami et al. (1992), and has since been a standard tool for analyzing borehole temperatures as paleoclimatic indicators.In this method, stabilization is achieved by keeping only the p largest singular values of M, using an appropriate threshold .A generalized inverse, M g p , can then be formed: In this equation, T obs is a vector of borehole temperature observations, USV T is the SVD of M, and p is the number of singular values λ i λ max > .Many methods exist to determine the optimal number of singular values (see Hansen, 1998Hansen, , 2010)).Herein we select the number of singular values based on a fixed value of = 0.025, which has been defined in previous studies according to the precision of borehole temperature measurements.
In practice, this method requires the determination of the equilibrium surface temperature T 0 , the geothermal heat flow density q 0 , the bottom boundary condition and the timevarying upper boundary condition from the measured T (z) data.Assuming uniform thermophysical properties of the subsurface, a simple estimate of the geothermal gradient can be used for inversions: q 0 is calculated from the the linear trend determined from the deepest part of the temperature profile, which is assumed to be least affected by the recent ground temperature changes.T 0 then follows from upward continuation of this linear trend.

Impact of postglacial warming on the determination of the steady-state heat flux profile
Here, we use basal temperature values from the MUN-GSM simulation to quantify the extent of the perturbation resulting from postglacial warming to steady-state temperature profiles.We then use this simulated perturbation to derive spatial maps of the estimated impacts on GST inversions and estimates of subsurface energy storage over North America that do not account for the LGC perturbation in their analysis.In order to assess the subsurface temperature anomaly field induced by the ensemble of the GSM basal temperature field, we express the temperature time series (120 time steps of 1000 yr with constant temperatures) at each grid point as an anomaly from the local mean over the entire simulation surface temperature variation and its corresponding uncertainty envelope.period.We treat the 75 % probability intervals for each time series in the same manner.The gridded basal temperature anomalies are then used as a forcing function, i.e., the timevarying upper boundary condition, in the forward model expressed in Eq. ( 4) to estimate the perturbations to the subsurface thermal field resulting from the MUN-GSM 120 ka surface temperature variation and its corresponding uncertainty envelope.
Basal temperatures are taken from five sites in a latitudinal transect as represented in Fig. 1 that also approximate the locations of observational measurements in boreholes.The basal temperature time series from these sites are shown in Fig. 2a as departures from the respective site-specific means.Figure 2b plots the modeled subsurface temperature anomalies associated with the selected basal temperatures as functions of depth.A thermal diffusivity of κ = 10 −6 m 2 s −1 was assumed to be representative of common nonporous crustal rocks (e.g., Drury, 1986;Beardsmore and Cull, 2001;Pasquale et al., 2014).Each of the basal temperature histories cause perturbations within the upper 1000 m of the selected profiles that could obscure the correct estimate of the geothermal gradient ( = dT dz ) as a function of depth in BTP studies.Figure 2c further plots the perturbation to the heat flow density (q = −λ dT dz ) as a function of depth under the as- sumption of a homogeneous subsurface and a constant thermal conductivity λ = 3.0 W K −1 m −1 .No efforts were made to take into account any geological differentiation, vertically or horizontally.In this respect, it is even less developed than the subsurface module of MUN-GSM, which includes the effects of a simplified geology (sediment versus crystalline), laterally varying heat flux, the TTOP thermal offset, and a simple permafrost model (Tarasov and Peltier, 2007).This should be included in future investigations, as variations in subsurface porosity, heat capacity, and thermal conductivity have a significant impact on thermal diffusivity and permafrost-related processes.Thus, the derived quantities should be seen as low-order estimates.
The subsurface thermal effect of the LGM is spatially variable because of differences in the advance, extent and retreat of the Laurentide ice sheet.Figure 3a, b shows the spatial distribution of the vertically averaged heat flux perturbation for two depth ranges of 100-200 m and 500-600 m, respectively, while Fig. 3b, d illustrates the spatial distribution of the associated 75 % uncertainty intervals obtained from the MUN-GSM variability.Positive perturbations refer to heat gain by the ground subsurface.Negative ones refer to ground heat loss.According to the MUN-GSM simulations, the heat flux perturbation is nearly constant for the depth range and locations at which the great majority of the temperature versus depth profiles in the borehole climatology database (Pollack et al., 1998)  Figure 4b.For consistency, we used the same synthetic history as described in Beltrami et al. (2011) over the full domain, despite the fact that these histories vary spatially over our experimental domain.
The spatially homogenous last-millennium history nevertheless simplifies our interpretation of LGC impacts on BTP analyses by isolating only the spatial heterogeneity of the LGC basal temperature 320 changes within our analysis.
We restrict ourselves to three examples of 120-ka basal temperature histories as upper boundary conditions in our analysis.These curves are shown in Figure 4a. Figure 4c shows the combined 119 ka of the glacial-deglacial temperature plus the 1000-year artificial history from Beltrami  the heat flux perturbations are not spatially homogeneous: positive heat flux perturbations are persistent and of the same magnitude in eastern Canada and the northeastern US, while negative anomalies are apparent in the eastern regions of the Canadian Arctic.

Impact of postglacial warming on the recent GSTH and subsurface heat content
Section 3 quantified the transient effect of the LGC on subsurface thermal fields using the basal temperatures simulated by the MUN-GSM.In the subsequent section, we investigate how the LGC signal may impact temperature and energy flux histories estimated for the last millennium from BTP analyses that do not account for the LGC impacts.All global analyses of BTPs (Pollack et al., 1996;Harris and Chapman, 2001;Beltrami and Bourlon, 2004) have assumed that the influence of the LGC on BTPs was negligible within the first 600 m of the profiles, the maximum depth that is typically used to estimate last-millennium surface changes.We therefore investigate how studies making these assumptions may or may not include additional biases due to the unaccounted LGC impact.
Our approach assumes a synthetic GST history for the last millennium of the form shown in Fig. 4b.For consistency, we used the same synthetic history as described in Beltrami After constructing the synthetic BTPs as described above, the profiles are processed in exactly the same manner as real-world observations.The lower section of each BTP is fit to a linear function, which is subtracted from the whole BTP to form a temperature anomaly.Figure 6 shows three et al. ( 2011) over the full domain, despite the fact that these histories vary spatially over our experimental domain.The spatially homogenous last-millennium history nevertheless simplifies our interpretation of LGC impacts on BTP analyses by isolating only the spatial heterogeneity of the LGC basal temperature changes within our analysis.We restrict ourselves to three examples of 120 ka basal temperature histories as upper boundary conditions in our analysis.These curves are shown in Fig. 4a. Figure 4c shows the combined 119 kyr of the glacial-deglacial temperature plus the 1000 yr artificial history from Beltrami et al. (2011).The resulting present-day subsurface temperature anomalies generated by the one-dimensional forward model assuming isotropic thermophysical properties are shown in Fig. 5a.A long-term surface temperature of 8 • C and a constant geothermal gradient of 20 K km −1 are added to these anomalies to generate a set of simulated BTPs based on the simulated locations from the MUN-GCM.These background steady-state conditions were chosen after Beltrami et al. (2011) for consistency.An example of one completed BTP is shown in Fig. 5b.Note that similar BTPs were calculated at each grid point across the modeled domain and used in the standard inversion procedure described in Sect.2.2.
After constructing the synthetic BTPs as described above, the profiles are processed in exactly the same manner as realworld observations.The lower section of each BTP is fit to a linear function, which is subtracted from the whole BTP to form a temperature anomaly.Figure 6 shows three examples examples of such anomalies as a function of depth, where the solid and dotted lines of each colour represent the mean subsurface temperature anomaly at eachaux location and 75% uncertainty range of the anomalies, respectively.The solid pink lines in each profile are the true subsurface temperature anomaly for the synthetic 1000-a GSTH used as a reference.The positive biases on the subsurface anomalies that extend to about 400 m are consistent across the profiles, and appear to have less impact on the shape of the temperature perturbation as a function of depth, but more implication for estimates of the subsurface energy storage that is estimated from the integrated area between the profiles and a line of zero anomaly.
The temperature anomaly in Figure 6 is consistent with the overall spatially-averaged subsurface heat content of the domain shown in Figure 7, where the cumulative integrals from the surface to a given depth are shown for the case of LGC plus the last millennium history (black), synthetic last millennium history (red), and their difference (i.e., LGC contribution) (green).These results clearly indicate that heat content estimates over the complete domain in this study can differ by as much as 50% for cases in which the LGC impact is not addressed.of such anomalies as a function of depth, where the solid and dotted lines of each color represent the mean subsurface temperature anomaly at each location and the 75 % uncertainty range of the anomalies, respectively.The solid pink lines in each profile are the true subsurface temperature anomaly for the synthetic 1000 a GSTH used as a reference.The positive biases in the subsurface anomalies that extend to about 400 m are consistent across the profiles, and appear to have less impact on the shape of the temperature perturbation as a function of depth, but more implication for estimates of the subsurface energy storage that is estimated from the integrated area between the profiles and a line of zero anomaly.The temperature anomaly in Fig. 6 is consistent with the overall spatially averaged subsurface heat content of the domain shown in Fig. 7, where the cumulative integrals from the surface to a given depth are shown for the case of LGC plus the last millennium history (black), synthetic last millennium history (red), and their difference (i.e., LGC contribution) (green).These results clearly indicate that heat content estimates over the complete domain in this study can differ by as much as 50 % for cases in which the LGC impact is not addressed.
The subsurface cumulative heat integrals do not change much below 300 m for this specific experiment.The subsurface energy contribution from the LGC according to these estimates for the first 150 m is about 1 ZJ (1 ZJ = 10 21 J), which is about one-tenth of the energy stored in the subsurface for the second half of the 20th century, estimated from the uncorrected global database by Beltrami (2002a) for all continental areas, and about one-fifth of that estimated for the Northern Hemisphere's land (Beltrami et al., 2006); clearly, a non-negligible quantity of energy.The overall cumulative subsurface energy contribution to the subsurface energy content as a single quantity is shown in Fig. 7; however, there are important spatial variations, as expected from the variability of the LGM impact at the surface.The spatial variability is displayed in Fig. 8, where the mean subsurface energy contributions are shown for several depth ranges as stated in the caption.Figure 9 shows the corresponding 75 % uncertainties.
These results represent the spatial distribution of the potential subsurface energy that has remained unaccounted for (Beltrami, 2002a;Beltrami et al., 2006) in previous estimates of the ground energy contribution from geothermal data to the energy balance of the climate system.The magnitude of such contributions could be important, as climate models attempt to include surface processes and soil thermodynamics in their simulations.This is particularly important, as spatial-scale resolution continuously increases in nested regional models; that is, local effects recorded in boreholes may be helpful for ascertaining surface or subsurface processes at locations where these processes can be modeled on small scales.
To date, most borehole climate reconstructions have dealt with estimating the ground surface temperature changes over the last millennia, rather than estimating the subsurface energy content.Therefore, we have also estimated the temporal variability of the upper boundary condition from inversion at Three examples of the results from inversion are shown in Fig. 10.Here, the pink lines represent the inversion results for the true synthetic last millennium history, and the color lines represent the inversion results for the LGC plus the last millennium cases.The color ranges represent the 75 % uncertainty intervals for the LGC and last millennium cases.In Fig. 11, we show the spatial distribution of the resulting deviations of estimated last millennium histories at different time intervals for each grid point.Note again that we have assumed a constant last millennium history at all locations in the domain; thus, the resulting spatial variations are due entirely to differences in the LGC basal histories locally.

Discussion and conclusions
In this paper, we assess the potential magnitude of the thermal energy contributions from the LGC to the shallow subsurface for those regions that had significant ice coverage during the LGC.In this context, the results from our experiments indicate that the subsurface disturbances to the semiequilibrium geothermal gradient are not large at depths from the surface to 600 m, where most of the data in the borehole climatology database used in existing global studies are found.GSTHs from previous studies would show qualitatively similar characteristics in their temporal evolution, although the magnitude of the temperature changes would be different, because the energy contributions from the LGC distort the quasi semi-equilibrium geothermal gradient, and thus alter the reference against which GST changes are estimated.Although some significant changes in the geothermal gradient within this depth range are observed at some locations (e.g., the magenta line in Fig. 2), they usually occur at higher latitudes where differences may be due to permafrost and active layer phenomena, so that the purely conductive air-ground coupling assumption in our forward model does not translate well into a one-dimensional conductive subsurface model.
Thus, according to our experimental estimates of the effect of the LGC on the thermal regime of the shallow subsurface from the MUN-GSM, the global temperature reconstructions (Pollack et al., 1996;Harris and Chapman, 2001;Beltrami and Bourlon, 2004) will not be largely affected.The effect of the LGC on our synthetic GSTH example appears to be small for the last 150 yr.From 1500 to 1850 CE, its effect amounts to a suppression of the warming; thus, the most relevant effect may be a larger magnitude of warming for the recovery from the Little Ice Age (LIA) than what BTP analyses have already estimated.This effect nevertheless would be expected to have the largest influence on only a fraction of BTPs at the northern latitudes, in other words, those areas that have been under the ice sheet, and those that have the subsurface signature of the LIA.
Our results present, for the first time, an estimation of the energy contribution to the ground subsurface from the LGC.The spatial variability of the subsurface energy remaining from the LGC is significant.Over the complete experimental domain, the integral for the cumulative subsurface heat content due to the thermal effects of the LGC is about 1 ZJ for the first 100 m and 2 ZJ for the first 300 m.The magnitude of subsurface thermal perturbation does not change significantly

Discussion and conclusions
In this paper, we assess the potential magnitude of the thermal energy contributions from the LGC to the shallow subsurface for those regions that had significant ice coverage during the LGC.In below this depth, because our experiment is limited to the maximum depth of existing global GSTH analyses (600 m).However, this quantity of heat is on the same order of magnitude as the heat absorbed by the ground in the last 50 yr of the 20th century, which was estimated from geothermal data to be 4.8 ZJ for the Northern Hemisphere, and about 8 ZJ for all continental areas except Antarctica.The total heat absorbed by the continental areas of the NH since the beginning of industrialization is estimated to be 13.2 ZJ. 36 % of this heat gain occurred in the last 50 yr of the 20th century (Beltrami, 2002b;Beltrami et al., 2006).
Recent calls have noted the importance of monitoring the temporal rate of change in the energy stored in climate system components (Hansen et al., 2011(Hansen et al., , 2013)).Although most of the energy changes will be associated with the ocean (Kosaka and Xie, 2013), where the majority of the energy has been stored (Rhein et al., 2013;Goddard, 2014), it remains relevant to monitor the energy stored in other climate subsystems, including the continental subsurface.
To explore the potential of geothermal data in estimating future continental heat storage, we assume a conservative estimate of a 1 K increase in surface temperature, and examine the expected changes in the continental heat content for this situation.The time derivative of the temperature for a step change in the surface temperature is For a 1 K increase at 100 yr BP, this function has a maximum value at 75 m, and over a 10 yr period, the temperature change will be 0.024 K.For the same change at 50 and 20 yr BP, the maximum rates of change will be at 56 and 36 m, with changes of 0.05 and 0.12 K over 10 yr, respectively (Mareschal and Beltrami, 1992).For the warming expected at most locations in central and eastern Canada (Bel- that the purely conductive air-ground coupling assumption in our forward model does not translate 395 well into a one-dimensional conductive subsurface model. Thus, according to our experimental estimates of the effect of the LGC on the thermal regime of the shallow subsurface from the MUN-GSM, the global temperature reconstructions (Pollack et al., 1996;Harris and Chapman, 2001;Beltrami and Bourlon, 2004) will not be largely affected.
The effect of the LGC, on our synthetic GSTH example, appears to be small for the last 150 years.trami and Bourlon, 2004), the maximum temperature change in the subsurface temperature profile should be on the order of 50 mK in a decade, and the detection of such changes is already achievable.
In order to estimate the expected future rate of change of the heat stored in the subsurface due to changes in the surface temperature, we can assume a 1 K step change taking place at present; thus, the subsurface heat gain in the first 600 m will be on the order of 60, 134 and 190 MJ m −2 for 10, 50 and 100 a into the future, respectively.These correspond to mean rates of subsurface heat changes of about 190, 86, and 60 mW m −2 , respectively.Such changes are significantly larger than most of the energy fluxes measured from the interior of the Earth, and should be easily detectable, particularly in areas of the Northern Hemisphere where the projected warming is expected to be much higher than 1 K.
Comparing these numbers with the estimates provided by Beltrami et al. (2006), the mean fluxes for the period between 1780 and 1980 from existing borehole temperature data are about 21 and 30 mW m −2 , with the estimated cumulative heat flux absorbed by the ground since the start of industrialization at about 100 mW m −2 .
Preliminary subsurface heat content estimates from borehole temperature data in the experimental domain for North America are 1 ZJ from the surface to a depth of 100 m obtained from 372 BTPs, and 1.96 ZJ to a depth of 200 m obtained from 260 BTPs.The LGC cumulative contributions to subsurface heat content (green curve in Fig. 9), from the ground surface to the above depths, are 0.6 ZJ and 1.2 ZJ, respectively.These results indicate that the contribution of LGC to the current subsurface heat content, according to the MUN-GSM, is about 60 % of the measured energy, and cannot be ignored.
Overall, our results suggest that it is potentially important to monitor the temporal variation of temperature profiles.In principle, such observations would allow (1) a determination of the rate of change of continental energy storage over time periods as short as one decade, while minimizing the uncertainties, (2) the acquisition of robust GSTH estimates by differentiating between stationary and transient perturbations to the geothermal gradient, as well as increasing the GSTH model resolution (Beltrami and Mareschal, 1995), and (3) assistance in decreasing the uncertainties of the thermal regime of the soil, of which a proper characterization plays an important role in the modeling of near-surface phenomena, such as those governing the spatial distribution of permafrost and the associated active layer changes, as well as those related to soil carbon and its stability in a future climate.
The Supplement related to this article is available online at doi:10.5194/cp-10-1693-2014-supplement.

Fig. 1 .
Fig. 1.Simulation domain (172.5 o W, 34.75 o N) and (52.5 o W, 75 o N) and grid cell locations within the MUN-GSM.Colored stars mark the locations of the five grid cells from which simulated temperatures are used forsubsequent studies (see, e.g.Figure2).

Figure 1 .
Figure 1.Simulation domain (172.5 • W, 34.75 • N and 52.5 • W, 75 • N) and grid cell locations within the MUN-GSM.Colored stars mark the locations of the five grid cells from which simulated temperatures are used for subsequent studies (see, e.g., Fig. 2).

Fig. 2 .
Fig. 2. (a) Basal temperature evolution simulated by the MUN-GCM for the last 120Ka.(b) Resultant presentday subsurface temperature perturbations associated with the basal temperature time series shown in panel a.(c) Resultant present-day subsurface heat flux perturbations associated with the basal time series shown in panel a.

275Figure 2 .
Figure 2. (a) Basal temperature evolution simulated by the MUN-GCM for the last 120 ka.(b) Resultant present-day subsurface temperature perturbations associated with the basal temperature time series shown in panel (a).(c) Resultant present-day subsurface heat flux perturbations associated with the basal time series shown in panel (a).

Figure 3 .
Figure 3. Heat flux anomalies for depth intervals averaged over depth intervals of (a) 100-200 m and (b) 500-600 m.Panels below represent the corresponding distribution of the absolute value of the 75 % uncertainty of each of the above heat flux perturbations.Additional results for the full depth interval down to 2000 m are given in the Supplement to this article at http://editor.copernicus.org/index.php?.

Fig. 4 .
Fig. 4. (a) Surface forcing function: Temperature change from GSM for three of given sites.(Insert b) Theoretical GSTH for the most recent 1000 years to be used as control.(c) Superposition of temperature changes from (a) and (b), i.e., the forcing variation at the surface. 13

Figure 4 .
Figure 4. (a) Surface forcing function: temperature change from the GSM for three of the given sites.(b) Theoretical GSTH for the most recent 1000 yr to be used as a control.(c) Superposition of temperature changes from (a) and (b), i.e., the forcing variation at the surface.

Fig. 5 .
Fig. 5. (a) Subsurface temperature anomalies from the forcing functions in Figure 4c.(b) To simulate the synthetic temperature versus depth profile shown, we added to the anomalies in (a) an assumed steady-state with a surface temperature of 8 • C and a geothermal gradient of 20, K km −1 (black solid line).Red dots denote resulting the steady state as modified by the subsurface anomalies in (a).et al. (2011).The resulting present-day subsurface temperature anomalies generated by the one-325

335 14Figure 5 .
Figure 5. (a) Subsurface temperature anomalies from the forcing functions in Fig. 4c.(b) To simulate the synthetic temperature versus the depth profile shown, we added to the anomalies in (a) an assumed steady state with a surface temperature of 8 • C and a geothermal gradient of 20 K km −1 (black solid line).Red dots denote the steady state as modified by the subsurface anomalies in (a).

Fig. 6 .
Fig. 6.Color lines represent the anomaly profiles from functions in Figure 4c.Dotted lines denote the 75% uncertainty range as generated from the MUN-GCM simulations.Pink lines represent the perturbation from the artificial GSTH in Figure 4b.

15Figure 6 .
Figure 6.Color lines represent the anomaly profiles from functions in Fig. 4c.Dotted lines denote the 75 % uncertainty range as generated from the MUN-GCM simulations.Pink lines represent the perturbation from the artificial GSTH in Fig. 4b.

Fig. 7 .Figure 7 .
Fig. 7. Mean heat content change over the experimental domain as a function of depth.Points represents the heat content integral from the surface to a given depth with respect to geothermal equilibrium.Contributions from the GSTH and (LGC + GSTH) are represented by the red and black curves, respectively.The green curve represents their difference, that is, the contribution of LGC to subsurface heat content.

Figure 8 .
Figure 8.Heat content contribution per unit area from the LGC.Shown are the estimates of the contribution to the subsurface thermal energy from long-term climatic changes.Results are shown for depth ranges of (a) 0-80, (b) 0-160, (c) 0-320, and (d) 0-600 m.

Figure 9 .
Figure 9. Spatial distribution of the 75 % uncertainties in subsurface thermal energy storage for the corresponding panels in Fig. 8.

Fig. 10 .
Fig. 10.Sample mean GSTH from inversion of subsurface temperature anomalies.Dashed lines denote the 75% confidence intervals.Pink lines denote the "true" GSTH obtained using the theoretical function of Figure 3. bution of the resulting deviations of estimated last millennium histories at different time intervals for each grid point.Note again that we have assumed a constant last millennium history at all locations in the domain, thus the resulting spatial variations are due entirely to differences in the LGC basal histories locally.

385Figure 10 .
Figure 10.Sample mean GSTH from inversion of subsurface temperature anomalies.Dashed lines denote the 75 % confidence intervals.Pink lines denote the "true" GSTH obtained using the theoretical function of Fig. 3.

Fig. 11 .
Fig. 11.Spatial distribution of differences in estimated GSTH with and without the influence of the LGC; i.e. blue areas are those where reconstructions are colder than the true solution, and red areas are those where they are warmer.Shown are the GST differences over the time intervals from (a) 0-150 years BP, (b) 150-300 years BP, (c) 300-450 years BP, and (d) 450-600 years BP.

400Figure 11 .
Figure 11.Spatial distribution of differences in estimated GSTH with and without the influence of the LGC; i.e., the blue areas are those where reconstructions are colder than the true solution, and red areas are those where they are warmer.Shown are the GST differences over the time intervals from (a) 0 to 150 yr BP, (b) 150 to 300 yr BP, (c) 300 to 450 yr BP, and (d) 450 to 600 yr BP.