Impact of the Last Glacial Cycle on Late-Holocene temperature and energy reconstructions from terrestrial borehole temperatures in North America

H. Beltrami, G. S. Matharoo, L. Tarasov, V. Rath, and J. E. Smerdon Climate & Atmospheric Sciences Institute and Department of Earth Sciences, St. Francis Xavier University, Antigonish, Nova Scotia, Canada Department of Physics and Physical Oceanography, Memorial University of Newfoundland, St. John’s, Newfoundland, Canada Universidad Complutense de Madrid, Facultad CC. Físicas, Departamento de Astrofísica y CC de la Atmósfera, Madrid, Spain Dublin Institute for Advanced Studies, School of Cosmic Physics, Geophysics Section, Dublin, Ireland Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York, USA


Introduction
The past 2000 years are an important time period for data-model comparisons as a means of assessing past variability and change, as well as evaluating the performance of the 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 years by indirect methods, instrumental records rarely exceed 100-200 years 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;González-Rouco et al., 2009, for reviews).These records have most widely been used to estimate past temperature changes, but also can 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 and ill-posed and assumptionladen inverse problem as it is the case for GST history inversion, but can be estimated from observed BTPs with only minimal assumptions on the volumetric heat capacity.conductive processes dominate, changes in 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 years 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 Earth's major climate subsystems (ocean, atmosphere, cryosphere, and continents) from the perspective of climate system dynamics (Pielke, 2003;Hansen et al., 2005Hansen et al., , 2011Hansen et al., , 2013;;Trenberth, 2009).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 process important for risk assessments associated specifically with soil processes and function 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 on the continents, comparisons between paleoclimatic GCM simulations and borehole Figures estimates of past GSTs or energy fluxes comprise a valuable opportunity for evaluating GCM performance.These comparisons are nevertheless dependent on the robust characterization of past changes in the subsurface energy reservoir and 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 the Last Glacial Cycle (LGC), including the Last Glacial Maximum (LGM, ≈ 26-21 ka BP) 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 years).This steady-state signal is approximated below a few hundred meters in BTPs by a linear function of increasing temperatures with depth.The thermal effects of the LGC have been previously discussed as sources of error 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).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 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 at 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 the complexities of Figures 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.
Herein we study the effect of the development and demise 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.1, 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.3.2 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 ka.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 temperature-depth 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., Figures Back Close Full 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 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 basis of BTP modeling and interpretation that are relevant to our subsequent analysis.Our first focus, however, is a more detailed description of the Laurentide ice-sheet model that we employ to estimate the ground surface temperature history over North America over the LGC.

Constructing the upper boundary condition
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 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).This approximation, 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 stated derived from PMIP II archived 21 ka BP slice GCM simulations Figures

Back Close
Full (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 between (172.5 • W, 34.75 • N) and (42.5 • W, 84.75 • N), with a resolution of 1 • longitude by 0.5 • latitude, covering most of North America, at constant 1000 year time intervals from 120 ka BP 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 pro-glacial 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 (Tarasov and Peltier, 2007).
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, endeavouring 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 predominately ice-free.For the duration of ice cover, given the constraints and represented physics, the calibrated results offer the best available reconstruction of sub-glacial 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)  avoids the assumption of an underlying Gaussian distribution of misfits required for the standard 96 % assignment.

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 homogeneous, isotropic half space (Carslaw and Jaeger, 1959) as: where κ is the thermal diffusivity of the rock defined as κ = λ ρc (m 2 s −1 ), λ is the thermal is the volumetric heat production (W m −3 ), z is depth (m), and t is time (s).We assume the common convention that z is positive downwards.An analytic solution for Eq. ( 1 Given the above described boundary conditions, 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 equilibrium, steady-state temperature profile often called the "geothermal gradient", while T t (z, t) represents the perturbation by transient changes in the surface temperature boundary condition (Mareschal and Beltrami, 1992).The geothermal gradient 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 time scale 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 downward propagating 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 Q s = ρc z 0 T (z)dz.Here, ρc is the volumetric heat capacity, chosen here as 3 × 10 −6 J m −3 K −1 after Čermák and Rybach Introduction

Conclusions References
Tables Figures

Back Close
Full (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) thus can be inverted to estimate T gst given measurement of T obs .The equation is nevertheless not easily invertible because the problem is ill-posed (see Hansen, 1998;Aster et al., 2005, for a detailed introduction).In the case considered herein, the truncated singular value decomposition (SVD, see Lanczos, 1961;Hanson and Lawson, 1969;Varah, 1973) 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 > .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 time varying upper boundary condition from the measured T (z) data.Assuming uniform Introduction

Conclusions References
Tables Figures

Back Close
Full 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.
3 Experimental setup, results and discussion

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 year mean temperatures) at each grid point as an anomaly from local mean over the entire simulation 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 functions, i.e. the time-varying 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 120 ka basal temperature time series from these sites are shown in Fig. 2a as departures from the respective site-specific means.Figure 2b  modeled subsurface temperature anomalies associated with the selected basal temperatures as a function of depth.A thermal diffusivity of κ = 10 −6 m 2 s was assumed as 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 assumption of a homogeneous subsurface and a constant thermal conductivity 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 3 shows the spatial distribution of the perturbation of the heat flux over the full domain of Fig. 1.Each sub panel in Fig. 3 shows heat flux perturbations for the depth ranges 30-500, 500-1000, 1000-1500, and 1500-2000 m, while eastern US, while negative anomalies are apparent in eastern regions of the Canadian Arctic.

Impact of postglacial warming on the recent GSTH and subsurface heat content
Section 3.1 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. 6b.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 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. 6a. Figure 6c shows the combined 119 ka of the glacial-deglacial temperature plus the 1000 year 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. 7a.A long-term surface temperature of 8 • C and a constant geothermal gradient of 20 K km (2011) for consistency.An example of one completed BTP is shown in Fig. 7b.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 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 8 shows three 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 each 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 year 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 Fig. 8 is consistent with the overall spatially-averaged subsurface heat content of the domain shown in Fig. 9, 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.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. 9, 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. 10, where the mean subsurface energy contributions are shown for several depth ranges as stated in the caption.Figure 11 shows the corresponding 75 % uncertainties.
These results represent the spatial distribution of the potential subsurface energy that has remained unaccounted (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 into 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 at 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 each grid point using the synthetic borehole temperature profiles generated for two sets of artificial data; one of them generated from the forward model of the synthetic last millennium history, and another generated from the LGC merged with this time series.millennium cases.In Fig. 13 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 semi-equilibrium 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 on 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.magenta line in Fig. 3), 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 year.Between CE 1500 to 1850, its effect amounts to a suppression of the warming, thus the most relevant effect may be a larger Introduction

Conclusions References
Tables Figures

Back Close
Full 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 are 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 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 of the same order of magnitude as the heat absorbed by the ground in the last 50 year 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 begining of industrialization is estimated to be 13.2 ZJ. 36 % of this heat gain occurred in the last 50 years of the 20th century (Beltrami, 2002a;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 1 K increase in surface temperature and examine the expected changes in continental heat content for this situation.The time derivative of the temperature for a step change in surface temperature is: Introduction

Conclusions References
Tables Figures

Back Close
Full For a 1 K increase 100 ka BP, this function has a maximum value at 75 m and over a 10 year period, the temperature change will be 0.024 K.For the same change 50 and 20 year BP, the maximum rate of change will be at 56 and 36 m, with changes of 0.05 and 0.12 K over 10 year (Mareschal and Beltrami, 1992).For the warming expected at most locations in central and eastern Canada (Beltrami and Bourlon, 2004), the maximum temperature change in the subsurface temperature profile should be of the order of 50 mK in a decade, and 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 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 of the order of 60, 134 and 190 MJ m −2 for 10, 50 and 100 year 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 flux for the period between 1780 and 1980 from existing borehole temperature data is about 21 and 30 mW m −2 with the estimated cumulative heat flux absorbed by the ground since the start of industrialization to 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's cumulative contribution to subsurface heat content (green curve in Fig. 9), from the ground surface to the above depths are 0.6 and 1.2 ZJ, respectively.These results Introduction

Conclusions References
Tables Figures

Back Close
Full 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 modelling of near surface phenomena, such as those governing the spatial distribution of permafrost and associated active layer changes, as well as those related to soil carbon and its stability in a future climate.Introduction

Conclusions References
Tables Figures

Back Close
Full       Figure 6b.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.calculated at each grid point across the modeled domain and used in the standard inversion procedure described in Section 2.2.
After constructing the synthetic BTPs as described above, the profiles are processed in exactly the 300 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 8 shows three examples of such anomalies as a function of depth, where the solid and dotted lines of each colour 50% for cases in which the LGC impact is not addressed.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 at small scales.
To date, most borehole climate reconstructions have dealt with estimating the ground surface Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | extension that uses sample mean and standard deviation.The assignment of a weaker interpretation of the 2σ confidence interval Discussion Paper | Discussion Paper | Discussion Paper | ) requires the specification of boundary conditions at the surface and at depth within the half space.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 time scales of millions of years, we assume a constant heat flux (Neumann type) boundary condition at infinity within the half space.The upper boundary condition 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.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | A discrete version of the transient part of Eq.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig. 4 plots the widths of the 75 % confidence intervals for the heat flux perturbations.Most of the mean heat flux perturbations are restricted to the first kilometer of the subsurface.This depth range spans the typical depth of BTPs used for borehole paleoclimatic studies.It is therefore important to estimate the impact of the LGC perturbations on GST and heat flow estimates from boreholes.Figures 5a and b show the spatial distribution of the vertically averaged heat flux perturbation for two depth ranges of 100-200 and 500-600 m, respectively, while Fig. 5b and 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 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 vs. depth profiles in the borehole climatology database (Pollack et al., 1998) are available.However, for these depth ranges, the heat flux perturbations are not spatially homogeneous: positive heat flux perturbations are persistent and of the same magnitude in eastern Canada and north Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |generate a set of simulated BTPs based on the simulated locations from the MUN-GCM.These background steady-state conditions were chosen afterBeltrami et al.
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 Discussion Paper | Discussion Paper | Discussion Paper | Three examples of the results from inversion are shown in Fig. 12.Here the pink lines represents the inversion results for the true synthetic last millennium history, and the colour lines represent the inversion results for the LGC plus the last millennium cases.The colour ranges represent the 75 % uncertainty intervals for the LGC and last Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |magnitude of warming for the recovery from the Little Ice Age (LIA) than what BTPs analyses have already estimated.This effect nevertheless would be expected to have the largest influence on only a fraction of BTPs in the northern latitudes.In other words, those areas that have been under the ice sheet, and those which have the subsurface signature of the LIA.
Discussion Paper | Discussion Paper | Discussion Paper |

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 for subsequent studies (see, e.g.Figure2).

Figure 1 .Fig. 2 .
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).

Figure 2 .Fig. 3 .
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 (a).(c) Resultant present-day subsurface heat flux perturbations associated with the basal time series shown in (a).

3. 2
Impact of postglacial warming on the recent GSTH and subsurface heat content.Section 3.1 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

Figure 5 .Fig. 6 .
Figure 5. 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 the absolute value of the 75 % uncertainty of each go the above heat flux perturbations.

285Fig. 7 .
Figure 6.(a) Surface forcing function: Temperature change from GSM for three of given sites.(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.

Figure 7 .Fig. 8 .Figure 8 .Fig. 8 .Fig. 9 .
Figure 7. (a) Subsurface temperature anomalies from the forcing functions in Fig. 6c.(b) To simulate the synthetic temperature vs. 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).

Figure 9 .Fig. 10 .Figure 10 .Fig. 11 .
Figure 9. 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. 335

Figure 11 .Fig. 12 .
Figure 11.Spatial distribution of the 75 % uncertainties in subsurface thermal energy storage for the corresponding panels in Fig. 10.

Figure 12 .Fig. 13 .Figure 13 .
Figure 12.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. 5.