Borehole climatology: a discussion based on contributions from climate modeling

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Borehole climatology: a discussion based on contributions from climate modeling J. F. González-Rouco, H. Beltrami, E. Zorita, M. B. Stevens

in climate reconstructions, as well as on the realism and limitations of model simulations. This paper explores research specifically related to paleoclimate modeling and borehole climatology as a branch of climate reconstruction that has contributed significantly to our knowledge of the low frequency climate evolution during the last five centuries. 10 The text flows around three main issues that group most of the interaction between model and geothermal efforts: the use of models as a validation tool for borehole climate reconstructions; comparison of geothermal information and model simulations as a means of either model validation or inference about past climate; and implications of the degree of realism on simulating subsurface climate on estimations of future climate 15 change.
The use of multi-centennial simulations as a surrogate reality for past climate suggests that within the simplified reality of climate models, methods and assumptions in borehole reconstructions deliver a consistent picture of past climate evolution at long time scales. Comparison of model simulations and borehole profiles indicate that bore- 20 hole temperatures are responding to past external forcing and that more realism in the development of the soil model components in climate models is desirable. Such an improved degree of realism is important for the simulation of subsurface climate and air-ground interaction; results indicate it could also be crucial for simulating the adequate energy balance within climate change scenario experiments.

Introduction
The purpose of this text is to illustrate and discuss advances at the interface of climate reconstruction studies and General Circulation Model simulations. The discourse specifically focuses on the climate of the last millennium and on borehole climatology applications as an example of blending of climate reconstruction and climate modeling 5 efforts. This section furnishes the general context in which these studies are rooted and establishes the structure of the sections to come.

Paleoclimate context
The last decades have witnessed an important growth in research efforts and resulting knowledge regarding climate variability of the last millennium. This has been parallel to 10 an increasing concern about the rising global and hemispherical temperatures and has offered important insights into present and future climate (Houghton, 2005;Solomon et al., 2007). This area of paleoclimate research has helped not only to place the relatively short instrumental record in a broader temporal context (Jones et al., 2001;Mitchell et al., 2001;Jones and Mann, 2004;Mann, 2007;North et al., 2006) but also to 15 explore the mechanisms behind climate events and periods that have been the subject of longstanding scientific interest like the Late Maunder Minimum within the Little Ice Age or the Medieval Warm period (Shindell et al., 2001;Bradley et al., 2003;Rind et al., 2004;Zorita et al., 2004). Such insight into pre-instrumental variability has come both from climate reconstruction exercises using indirect (proxy) sources of information and 20 from simulation of past climatic states with General Circulation Models (GCMs).
Numerical simulation has used models of varying complexity, from energy balance models (EBM, e.g. Crowley, 2000;Osborn et al., 2006) and Earth system Models of Intermediate Complexity (EMIC, e.g. Bauer and Claussen, 2006) to comprehensive atmosphere ocean GCMs (AOGCM, e.g. Zorita et al., 2004; several other research avenues of relevance for the understanding of past climate and future anthropogenic climate change such as: analysis of the response to natural and anthropogenic external forcing and mechanisms involved (e.g. Cubasch et al., 1997;Haigh, 1999;Rind et al., 2004;Zorita et al., 2005;Goosse et al., 2006;Tett et al., 2007;Ammann et al., 2007); comparison of model simulations and climate reconstructions 15 to detect the signal of the various external forcing factors through the last millennium (Crowley, 2000;Bauer et al., 2003;Hegerl et al., 2003Hegerl et al., , 2007a or to constrain estimates of climate sensitivity (Hegerl et al., 2006); use of model simulations as a surrogate reality in which pseudo proxies are built through the deterioration of simulations at the grid-point scale in order to test methods and assumptions relevant in proxy-based cli-Interactive Discussion EGU et al., 2003;Jones and Mann, 2004). Climate reconstruction efforts have targeted a wide range of spatial scales: from the local and regional (e.g. Cook, 1995;Overpeck et al., 1997;Luterbacher et al., 2004) to hemispherical and global scales (e.g. Briffa et al., 1998;Mann et al., 1998Mann et al., , 1999Esper et al., 2002;Moberg et al., 2005;Rutherford et al., 2005;Hegerl et al., 2007a). Though most efforts have been devoted to the 5 reconstruction of past temperature variability, a great deal of research has also concentrated on reconstructing hydrological and atmospheric circulation indices (e.g. Villalba et al., 1998;Diaz et al., 2001;Cook et al., 2002;Luterbacher et al., 2002;Pauling et al., 2006).
Hemispheric scale temperature reconstructions show general consistency in depict- 10 ing broad climate periods through the last millennium (Jansen et al., 2007) like the Medieval Warm Period (MWP, centered around 1000 AD) or the Little Ice age (ca. 1500 to 1850 AD). Additionally hemispheric temperatures also decreased during events of shorter duration like the Late Maunder Minimum (LMM, ca. 1700 AD), the Spoorer (SM, ca. 1450 AD) or the Dalton Minimum (DM, ca. 1800 AD) that are likely related to minima in solar activity. Nevertheless, quantitative agreement has been so far elusive and thus a matter of considerable debate in the last few years (e.g. von Storch et al., 2004Storch et al., , 2006Mann et al., 2005;Buerger and Cubasch, 2005;Buerger et al., 2006;Moberg et al., 2005;McIntyre and McKitrick, 2005;Wahl et al., 2006;Zorita et al., 2007;Mann et al., 2007a,b;Smerdon and Kaplan, 2007;Wahl and Ammann, 2007;Ammann and 20 Wahl, 2007) 1

.2 Borehole climatology
Within the context of millennial climate reconstructions, borehole temperature profiles (BTPs) have been one source of information that has significantly contributed to our understanding of centennial temperature changes (e.g. see Pollack and Huang, 2000; 25 Bodri and Cermak, 2007). Climate reconstruction based on BTPs leans on the assumption that surface air temperature (SAT) changes are coupled to ground surface temperature (GST) changes and propagate to the subsurface by thermal conduction 5 5 Contrary to other climate reconstruction procedures, temperature inversions obtained from borehole data are not calibrated against the instrumental record, thereby providing an independent measurement of past temperature. This approach has offered a rather singular view of the amplitude of global and hemispheric warming along the last five centuries Beltrami, 2002a;Pollack and Smer-10 don, 2004), although recent proxy reconstructions that have sought to preserve lowfrequency variance are consistent with borehole estimates (Esper et al., 2002;Moberg et al., 2005;Hegerl et al., 2007a). Nevertheless, the different magnitude of SAT changes inferred from borehole inversions and those based on other proxy reconstructions (Briffa and Osborn, 2002;Hegerl et al., 2007b) has fostered recent discussion 15 and examination of uncertainties in the borehole approach to climate reconstruction Mann and Schmidt, 2003;González-Rouco et al., 2003aChapman et al., 2004;Schmidt and Mann, 2004;Pollack and Smerdon, 2004) In addition to the recovery of the past temperature evolution BTPs have also been used to evaluate the role of heat storage in the lithospheric crust within the global 20 energy balance (Levitus et al., 2001;Pielke, 2003;Levitus et al., 2005;Hansen et al., 2005) by calculating the amount of energy stored in the ground due to GST warming (Beltrami, 2001b;Beltrami et al., 2002Beltrami et al., , 2006a. This information can be thought of as a byproduct of the temperature reconstruction analysis and, as such, it can be considered potentially subjected to the same advantages and shortcomings. 25 The geothermal approach, as any other method of inferring past climate, is not free from unknowns and limitations such as Majorowicz et al., 2004;Bodri and Cermak, 2007): a progressive inability with depth to resolve past climate variations due to the smearing nature of heat conduction; site specific non Introduction EGU climatic factors contributing to noise in BTPs (e.g., interaction of topography and hydrology, fluid flow and horizontal advection); long term variability in surface climate parameters which potentially disturb the SAT-GST relationship (e.g., freezing, vegetation, snow cover and evaporation changes); inhomogeneous distribution of properties related to sampling such as scarce and/or irregular spatial distribution of boreholes 5 at regional and larger scales, spatial variation in logging dates and depth of profiles, vertical resolution, uncertainties/errors in measurements, etc. These methodological and experimental problems have been addressed in a variety of studies, the balance of which suggest that such shortcomings can be addressed and that ultimately borehole profiles, if treated appropriately, can deliver reliable information about low frequency 10 climate evolution (e.g. Blackwell et al., 1980;Clow, 1992;Fernández and Cabal, 1992;Kukkonen and Clauser, 1994;Shen et al., 1995;Pollack et al., 1996;Kohl, 1998Kohl, , 1999Safanda, 1999;Serban et al., 2001;Bense and Kooi, 2004;Pollack and Smerdon, 2004;Nitoiu and Beltrami, 2005;Hartmann and Rath, 2005;Ferguson et al., 2006;Bense and Beltrami, 2007;Chouinard and Mareschal, 2007). 15 During the last few years climate model simulation studies have targeted questions related to borehole climatology from a variety of perspectives. As in the case of other reconstruction methodologies some studies have addressed aspects of the borehole method using model simulations as a surrogate reality in which the coupling between SAT and GST at long time scales can be addressed as well as some of the uncertain-20 ties mentioned above (e.g. Mann and Schmidt, 2003;González-Rouco et al., 2003b).
Other studies have paid attention to comparison of model simulations and borehole reconstructions or relevant thermo-physical subsurface properties, serving both the purpose of model validation and also the understanding of mechanisms contributing to climate variability (e.g. Beltrami et al., 2006b;Stevens et al., 2008). 25 Other questions involve the expected change of subsurface properties in climate change scenarios and also whether climate models address the physics of the subsurface to account for such processes in a realistic way (e.g. Stevens et al., 2007).
The purpose of this paper is to provide an integrating discussion that incorporates 8 research efforts within the field of borehole climatology and GCM applications oriented towards the simulation of soil thermodynamics, model-observation comparison or methodological evaluation of the borehole reconstruction approach. This research avenue integrating experimental and model work will be illustrated along the next few pages by providing in each section a brief background on the basis of which a discus-5 sion is furnished. For some specific purposes we will complement the ideas discussed in the text using data from the ECHO-g GCM simulations as well as standard geothermal models. Section 2 presents the GCM simulations and geothermal models used herein. Section 3 provides a perspective of past discussion and work concerning the validation of 10 some hypothesis and methodological issues of the borehole approach to climate reconstruction. Section 4 focuses on the realism of GCMs in reproducing the behavior of the geothermal observations and the subsurface thermal regime. Finally, Sect. 5 presents a brief discussion on the implications of the realism of GCMs in simulating subsurface geothermics within the context of climate change simulations.

Models and methods
GCM simulations and geothermal models have been used with the purpose of illustrating some of the discussions presented herein. A brief description for them is provided in the following subsections. 20 Model data were obtained from climate simulations produced with the ECHO-g atmosphere-ocean General Circulation Model (GCM). ECHO-g (Legutke and Voss, 1999) consists of the atmospheric and ocean GCM components ECHAM4 and HOPEg. ECHAM4 (Roeckner et al., 1996) is used with a T30 horizontal resolution (ca. 3.75 • ) and 19 vertical levels. HOPE-g (Wolff et al., 1997)  EGU tal resolution which varies towards the Equator where it reaches a minimum meridional grid point separation of 0.5 • for an improved representation of equatorial and tropical ocean currents. Vertical discretization for the ocean incorporates 20 levels.

GCM simulations
In order to avoid climate drift, heat and freshwater flux adjustments were applied to the ocean. These fluxes were diagnosed in a coupled spin-up integration with restor- 5 ing terms that drive the sea-surface-temperature and sea-surface salinity to their climatological observed values. These flux adjustments are constant in time through the integration and their global contribution is zero.
The surface scheme comprises a soil model, hydrology, snow cover physics and vegetation effects on surface evapotranspiration among others. The soil model, an 10 extension of Warrilow et al. (1986), is a five layer finite-difference approximation of the diffusion equation which operates on the T 30 land-sea-mask grid of ECHAM4. Ground temperatures are simulated at five levels with depths at 0.06 m, 0.32 m, 1.23 m, 4.13 m and 9.83 m. A zero heat flux is prescribed at the lowest layer in order to ensure that no artificial heat sources and sinks may affect the energy balance. This issue will be 15 further visited in Sect. 5 This work makes use of three integrations with the ECHO-g GCM: 1000 year long control simulation (CTRL) in which external forcings were fixed to the values of present climate and two forced simulations (FOR1 and FOR2) covering the period 1000 to 1990 AD produced driving the model with estimates of external forcing factors: at-20 mospheric greenhouse gas concentrations, solar irradiance and an estimate of the radiative effects of stratospheric volcanic aerosols. Sulphate aerosols or vegetation changes are not included in these simulations and some relative cooling effect from these forcing factors (Bauer et al., 2003;Osborn et al., 2006) should be expected to damp the simulated warming trends along the 20th century.  Figure 1a shows the estimation of external forcings used to drive the model. The atmospheric concentrations of CO 2 and CH 4 were obtained from analysis of air bubbles in Antartica ice cores (Etheridge et al., 1996(Etheridge et al., , 1998. Concentrations of N 2 O were used as in previous scenario experiments with this model (Roeckner et al., 1999): fixed 10 The effect of volcanic aerosols is incorporated as global effective variations of the an-5 nual values of the solar constant obtained from short wave radiative forcing changes (Crowley, 2000).
All previous forcings are subject to uncertainties, the largest of them being perhaps those related to the natural factors (volcanic and solar irradiance). CH 4 presents relatively small uncertainties (Etheridge et al., 1998); estimations of CO 2 changes are 10 found to be robust in new measurements (Etheridge et al., 1996;Siegenthaler et al., 2005), and N 2 O shows perhaps larger uncertainties though their impact is smaller. Volcanic forcing estimates show a high degree of variability according to different authors (e.g., Crowley, 2000;Robertson et al., 2001;Ammann et al., 2003); in addition to that, taking into account the impacts of volcanic events as an equivalent change in solar 15 constant does not allow for the specification of the place and moment of eruptions and the possible spatial heterogeneity in the spread of emitted aerosols. Solar irradiance changes are responsible for the largest impacts on centennial variability through the last millennium and the amplitude of their variations is under discussion (Bard et al., 2000;Lean et al., 2002;Foukal et al., 2004;Solanki and Krivova, 20 2004). The 10 Be Lean et al. (1995) rescaled estimations used in this case display changes of total solar irradiance between the LMM and present values of 0.30%, which are equivalent to those in other exercises (Bauer et al., 2003;Mann et al., 2005). Recent comparisons of NCAR CSM simulations using several past solar variability scenarios (Ammann et al., 2007) with NH temperature reconstructions are in good agree-25 ment for the scenario of lowest past solar variability (LMM to present changes of 0.1%) with temperature reconstructions showing a lower range of variability through the last millennium (e.g. Mann and Jones, 2003 EGU present of 0.25%), and agrees better with reconstructions showing larger changes in centennial temperature (e.g. Huang and Pollack, 1998;Moberg et al., 2005). Further description and discussion of uncertainties related to external forcing can be found in Gouiran et al. (2007a) and Zorita et al. (2007). Figure 1b shows the NH land temperature evolution in the three simulations. CTRL 5 presents fairly stable values along the 1000 years and FOR1 and FOR2 show a clear response to solar and volcanic forcing in the preindustrial era and the warming during the last centuries following the increase in solar irradiance and greenhouse gases. FOR 1 presents comparatively warmer values in the first centuries of the simulation.  and Osborn et al. (2006) suggest that FOR 1 is unusually warm 10 in comparison with other model simulations, a feature that can be related to warm initial conditions and the relatively short spin down to the forcing conditions of year 1000. Although the possibility of some initial imbalance can also not be ruled out in FOR 2, the level of medieval warming simulated in FOR2 can be supported with results from other more recent state of the art GCM (Ammann et al., 2007) and EMIC 15 simulations  that trace a very similar climate evolution through the last millennium. Analysis and validation of the ECHO-g performance in reproducing internal climate variability in CTRL in comparison with instrumental variability can be found in Min et al. (2005a,b); further details can be found in Raible et al. (2001Raible et al. ( , 2004Raible et al. ( , 2005 who analyses 20 North Atlantic atmospheric circulation in other control simulations with the same model. Analyses of different aspects of the forced simulations used herein can be found in the literature (von Storch et al., 2004;Zorita et al., 2003Zorita et al., , 2005González-Rouco et al., 2003a,b, 2006Beltrami et al., 2006b;Gouiran et al., 2007b,a;Stevens et al., 2007).
It is also interesting to see how the model simulates the thermal behavior of the 25 subsurface. Figure 1c shows NH land SAT averages in FOR 2 in comparison with temperature at the five soil model levels. Clearly, the highest frequencies have been filtered out and show considerable differences that will be discussed later in the text; at the lowest frequencies however subsurface temperatures and SAT share virtually the NH like Northern Africa (southern North America and southeastern Asia present warmer (colder) subsurface temperatures resulting from the availability of soil moisture and the balance of radiative and turbulent fluxes at the surface. These issues will be further discussed in Sect. 3.

10
Climate reconstruction of GST histories assumes that surface temperature variations propagate into the subsurface following the one dimensional time-dependent heat conduction equation (Carslaw and Jaeger, 1959): where κ is the thermal diffusivity, which we can assume constant for the typical range 15 of depths of interest for this work, z is the depth, and T is the temperature. In the absence of transient non-climatic disturbances at the ground surface or systematic perturbations that violate the one-dimensional heat transfer assumption (heat production, topography, ground water flow,variations in thermal properties, etc.) the temperature at depth z is given, at any time, by the combination of the geothermal temperature 20 gradient and the temperature perturbation T t (z) induced by the time varying surface temperatures: where q 0 is the quasi steady state surface heat flow density and R(z) is the thermal depth given by R(z)= Introduction EGU interval ∆z i . Subsurface temperatures can be modeled using surface temperature as an upper boundary condition of an infinite half space. In turn, temperatures at the surface can be thought of as a series of K step changes in temperature such that the subsurface signatures from each step change are superimposed and the temperature perturbation 5 at depth z is given by Mareschal and Beltrami (1992): where T k is the ground surface temperature, each value representing an average over time (t k -t k−1 ), and erf c is the complementary error function. This forward model can be used to derive grid point BTPs using simulated sur-10 face temperature changes as boundary conditions; as suggested by Fig. 1c and in González-Rouco et al. (2006) the selection of SAT or temperatures at the various soil model levels is of no importance. The forward model is used herein to derive gridpoint perturbation profiles 600 m deep using ECHO-g surface temperature changes as boundary conditions. At each 1 m depth interval a temperature anomaly is evaluated 15 from the contribution of surface temperature changes using the solution of the heat conduction equation in (2). Equation (2) is however sensitive to the reference temperature selected for the calculation of temperature changes (T k anomalies). This issue will be further discussed in Sects. 3 and 4. A singular value decomposition inversion model (SVD; Mareschal and Beltrami 20 (1992) is used to derive GST histories from grid-point BTPs. The thermal diffusivity value used for the forward and the inversion model applications was the same as in the ECHO-g GCM (0.75×10 −6 m 2 s −1 ). The model used for the inversion of each grid point BTP consists of a series of 20-yr step changes in GST history. The eigenvalue cutoff was set in different exercises between 0.025 and 0.3 to select the number of SVD 25 components Beltrami and Bourlon, 2004). The singular value cutoff is dependent in experimental practice on the model parameters, model 14 EGU geometry and data noise level, and it varies among BTPs of the same region. Generally, it is convenient to analyze BTPs with as similar geometry as possible to attempt to reduce the effects of different practical resolutions on the ensemble analysis.

Boreholes in the GCM surrogate reality
One strategy to test methods and assumptions in reconstruction approaches has been 5 to use GCM simulations as surrogates for the real climate evolution. Rather than representing the recent or past climate behavior, the simulations are meant to act as plausible climate realizations compatible with the external forcing imposed to the GCM, and complex enough to incorporate a variety of competing factors that contribute to the realism in the application of the target reconstruction technique (Mann and Rutherford,10 2002; Zorita and González-Rouco, 2002;Zorita et al., 2003;Rutherford et al., 2003;von Storch et al., 2004;Mann et al., 2005). Within the geothermal context, this approach has been used to explore potential caveats in the borehole approach to climate reconstruction through assessing various assumptions and methodological aspects. We will revisit some of this work along the following section and use it to further illustrate some 15 aspects relevant within the context of GST reconstructions from BTPs.
3.1 Discussing biasses on borehole temperature profiles Mann and Schmidt (2003) and Schmidt and Mann (2004) used a 50 year long simulation of the GISS model to explore the effect of snow cover changes on seasonal SAT and GST coupling. They find a good summer link between both variables and report 20 that the winter SAT-GST relationship is disturbed by snow cover. They therefore argue that climate changes in winter SAT might not be as well translated into GST changes as in the warm season. If long term changes in snow cover have occurred during the last centuries these could have caused anomalous changes in GST relative to SAT by altering winter ground surface insulation by snow cover (see also Stieglitz et al., 2003 Zhang, 2005). They argue that if a less intense hydrological cycle would have existed in the LIA due to a more frequent negative NAO phase (Shindell et al., 2001;Luterbacher et al., 2001) this would have been consistent with a reduction in winter precipitation and snow cover. Such scenario would produce anomalous cooling of GST relative to SAT in the coldest phases of the LIA (e.g. LMM). Thus, in the event that climate changes over 5 the past few centuries would have been seasonally different there would be potential for a bias in GST reconstructions that could reflect colder temperatures during LIA than might actually have existed. They suggested that interpretations of past SAT trends from borehole-based GST reconstructions may therefore exaggerate LIA cooling and this could account for the differences of ca. 0.5 K found among existing reconstructions 10 (e.g. Mann et al., 1999;Huang et al., 2000). The Mann and Schmidt (2003) arguments were debated by Chapman et al. (2004) using the same GISS simulation and by Bartlett et al. (2005) using observational snow cover, SAT and GST data. Chapman et al. (2004) argue that seasonal differences lead to seasonal offsets in GST and SAT averages and recommend that annual data 15 should be used in such an analysis. In doing so, they report good inter-annual tracking between GST and SAT through the examined 1951-1958 GISS simulated period. In a follow up study, Bartlett et al. (2005) argue that if the influence of snow on GST is not changing from year to year, such offset would be constant (e.g. Fig. 1d) and with no consequences for SAT and GST tracking. They agree that if snow cover is not 20 stationary its influence on the mean SAT-GST offset varies with time such that long term changes in snow cover may disrupt tracking between the two temperatures by introducing transient and persistent thermal signatures in the coupling. This possibility is examined by considering observations of snow cover changes in North America and their influence on the SAT-GST tracking. They quantified the changes in snow cover 25 onset and duration required to account for the approximately half a degree difference between borehole and multi-proxy climate reconstructions and report that it would be unlikely that changes in snow cover could account for this discrepancy. They find that in the 1970-2002 period, when winter SAT warming was the greatest, snow cover EGU inhibited 0.05 K/decade of seasonal warming from entering the ground, this being the result of a relatively stationary snow season under changing wintertime SAT.

EGU
These issues were also discussed using the millennial control and forced simulations with the ECHO-g model described in Sect. 2.1 (González-Rouco et al., 2003b by assessing the relationship between NH SAT and GST at low frequencies. All 5 model simulations showed differences in both variables at high (seasonal and interannual to decadal) frequencies but demonstrated strong coupling between air and ground temperatures at low frequencies essentially as shown in Fig. 1c. This suggested that, within the limits of the reality simulated by the model, changes in surface conditions that perturb ground-air coupling like snow cover or sensible and latent heat fluxes are 10 not strong enough to decouple NH SAT and GST at long time-scales.
Figures 2 to 4 allow a perspective on the spatial variability of SAT-GST in relation to surface processes like snow and soil moisture changes. Figure 2a shows the average spatial distribution of snow depth in the FOR 2 simulation. The mid-to northern-latitude areas of GST insulation due to snow cover and their correspondence with the negative 15 SAT-GST offsets ( Fig. 1d) are apparent. Figure 2b shows the correlation between annual snow depth changes and differences in annual SAT and GST (SAT minus GST at 0.06 m depth, T −0.06 m ) in the same simulation. Air-ground temperature differences are related with changes in snow depth at inter-annual timescales and broadly dominate in the extra-tropical regions, thus supporting the influence of snow cover at annual and 20 longer time scales: negative correlations indicate increases (decreases) in snow depth associated with larger (smaller) SAT-GST negative offsets. There is general agreement with the finding by Bartlett et al. (2005) that snow influence on the SAT-GST difference is greatest north of 45 • in North America. Figure 2c and d shows linear trends in simulated snow depth in an extended winter (DJFMA) season for the periods 1700-1990 Introduction  Figure 3 allows the discussion to be extended to the possible influence of soil moisture on air-ground temperature differences. The literature reports interesting changes in soil water content through the 20th century (e.g. Robock et al., 2000;Li et al., 2007;Sheffield and Wood, 2007) and, as in the case of snow depth/cover variations, these have the potential to change the ground temperature response by changing surface 5 latent heat fluxes and the heat capacity of the ground (Hu and Feng, 2005;Kueppers et al., 2007;Taylor et al., 2007). The average distribution of wetness ( Fig. 3a) depicts areas of larger water content over Europe, the western Siberian plains, Southeastern Asia and the tropical areas of America and eastern North America. Correlations between annual soil moisture and SAT-GST differences (Fig. 3b) illustrate the importance 10 of soil moisture changes on air-ground temperature coupling in the lower sub-tropical and tropical regions complementary to that of snow depth changes in the northern half of the NH (Fig. 2b). Soil moisture changes show opposite sense to those in absolute temperature differences at the surface, i.e., increases (decreases) in soil moisture tend to reduce (increase) the thermal difference between SAT and GST. In the areas of pos-15 itive (negative) SAT-GST offset in Fig. 1d like eastern USA and eastern Asia (northern Africa and western USA) this translates into negative (positive) correlations in Fig. 3b. Such relation suggests that soil moisture could also be a factor in producing long term changes in ground-surface coupling and thus a bias in borehole temperatures. Fig. 3c and d show linear regression coefficients for the period 1700 to 1990 and 1900 to 20 1990 AD in annual soil moisture as simulated in FOR 2. These values are in the range of observed estimations for the second half of the 20th century (e.g. Robock et al., 2000;Li et al., 2007). The simulated changes indicate for instance increases of soil moisture in northern Europe and decreases in mid-latitudes like the Mediterranean. These changes, as well as those shown above in snow depth would be consistent with 25 a northward shift of the hydrological cycle, thus, also compatible with a transition to a more zonal circulation (positive phase of the NAO) in the North Atlantic as hypothesized by Mann and Schmidt (2003). In fact, FOR 1 and FOR 2 show long term transitions to a more zonal circulation from the LIA to the end of the simulations (e.g. Zorita et al., EGU 2005).
Figures 2 and 3 present a geographically complementary view of the domains of relevance of snow depth and soil moisture for air-ground coupling as well as a non stationary view of long term changes in both variables. However, snow cover and soil wetness changes are not the only two factors with potential to disrupt the SAT- 5 GST relation. Other factors introducing long term changes in soil radiative insulation and surface turbulent and radiative fluxes (e.g., vegetation or land use changes) can potentially bias GST respect to SAT variations. Examples are vegetation and land use changes. These type of variations are not included in the simulations and could potentially affect GST at long time scales (Bonan et al., 1992;Lewis and Wang, 1998;10 Bauer et al., 2003;Nitoiu and Beltrami, 2005;Davin et al., 2007;Xu et al., 2007) and thus constitute issues which deserve attention in the future.
Long-term trends in SAT-GST differences during the simulated 20th century, (Fig. 4) show positive values indicating increases of SAT relative to GST over mid-latitudes and GST warming relative to SAT (negative trends) at northern latitudes. There are 15 no trends simulated in the areas of sensitivity of SAT-GST to soil moisture (Fig. 3b), thus the response suggests a relation with snow depth changes. Actually, there is regional agreement in broad spatial features in Figs. 2d and 4d like increasing snow depth (increasing the SAT-GST offset) over northern Eurasia and North America and decreasing snow depth (decreasing the SAT-GST negative offset) over the mid latitudes 20 and eastern Russia. Rather than long-term continuous changes, these linear trends arise in the simulations due to multi decadal changes along the second half of the 20th century. Figure 4c shows an example of SAT, GST and snow depth evolution at two characteristic points. The larger decreases (increases) in snow depth in the last decades in P 1 (P 2) lead to relative GST cooling (warming) with respect to SAT 25 changes due to increased exposure (insulation) to winter surface fluxes. The largest departures between SAT and GST anomalies go along with the largest changes in snow depth. While in P1 GST would under-estimate changes in SAT, in P 2 the contrary would happen. EGU Therefore, Fig. 4 shows that large changes in snow depth/cover can modify the SAT-GST relationship. The changes simulated by the model are reminiscent of those found in observations (e.g. Brown, 2000) in as much as they reflect larger changes of snow cover for some regions along the last decades of the 20th century. Such decadal trends should have little influence on the estimation of multi-century trends that are imprinted 5 on borehole BTPs to depths of several hundred meters. Thus, as suggested by Pollack and Smerdon (2004), it is not easy to make a case for a systematic bias in subsurface temperatures due to these processes. In the surrogate reality of the model simulations, in spite of the changes in all variables represented by Figs. 2 to 4, there seems to be a negligible impact in the reconstruction of NH GST histories from simulated BTPs 10 (González- Rouco et al., 2006).
It is worth noting that previous discussions or in methodological replications (e.g. González-Rouco et al., 2006) do not demonstrate that the borehole approach to climate reconstruction is correct. The arguments rather lean on a null hypothesis that the methodology is actually competent in its purpose and the evidence either reveals 15 itself consistent with this statement or rejects it if results point in the opposite direction. Also, in this line of argumentation, it is not of relevance whether the model simulations reproduce the exact real past climate trajectories of the involved variables. Though desirable, most likely this will not be the case since the evolution of these variables will be largely subjected to internal, non-forced, climate variability and due to limita-20 tions in model resolution and parametrizations (von Storch, 1995;Raisanen, 2007) to reproduce detailed aspects of snow depth and soil moisture dynamics (Brown, 2000;Déry and Wood, 2006;Robock and Li, 2006;Li et al., 2007). The relevant issue in this approach, however, is whether the magnitude and spatial variability of simulated changes are realistic enough to assess whether changes in SAT-GST coupling within 25 the simulated reality would produce long-term biases in BTPs. In the case discussed herein, the model simulations do not show evidence to distrust the skill of the method in retrieving hemispheric estimates of past temperature changes. At the regional scales, however, the weight of trends in snow depth or other surface perturbations might be 20 Introduction EGU larger and not lead to an overall balance as at hemispheric scales; an example addressing regional scales will be shown in Sect. 3.2.

Other methodological aspects
Beyond the questions raised above concerning SAT-GST coupling,  further address various issues concerning observational and methodological aspects 5 that could account for the differences between borehole-based and other proxy-based climate reconstructions. Two major problems were discussed: (i) the influence of spatial aggregation and gridding of GST reconstructions from BTPs and of latitudinal area weighting to produce a NH temperature reconstruction; (ii) the apparent failure in GST reconstructions to reproduce a spatial structure of changes comparable to that of SAT throughout the instrumental period.  developed a NH temperature reconstruction for the last five centuries using the Huang et al. (2000) set of inverted NH BTPs in an optimal signal detection approach which was more in agreement with previous multi-proxy reconstructions (Mann et al., 1999). These results were criticized by Pollack and Smerdon (2004) who examined the impact of aggregation and area weight- 15 ing for various grid resolutions on the spatial consistency of GST and SAT trends. They concluded that low occupancy cells obscured the relation between SAT and GST and demonstrated the existence of spatial consistency in both sets of variables. In a revision of their initial work Rutherford and Mann (2004) would later agree that gridding and area weighting do not resolve the differences between borehole and other proxy 20 reconstructions.
In addition to gridding and area-weighting, other possible methodological and observational aspects that may degrade borehole reconstructions as representative indicators of SAT changes were discussed by  and Pollack and Smerdon (2004). Two of these are borehole site geography and logging dates. Both of this 25 issues can produce unclear biases in the datasets that could affect GST reconstructions. The geographical distribution of boreholes is centered mostly in mid-latitudes and this could introduce a climatological bias in GST reconstructions. As for logging 21 Introduction EGU dates, the borehole data set  was compiled since the 1960's and different regions display some heterogeneity in logging dates (see Fig. 5c). This can produce biases by under-representing trends in regions where BTPs were logged earlier. While a geographical bias would potentially exaggerate warming by misrepresenting low latitude regions and oceans that would warm less, the bias in the logging 5 dates would presumably under-represent warming by excluding the later 20th century decadal warming in regions where BTPs were logged earlier (Pollack and Smerdon, 2004). The potential bias in the BTP geographical distribution was explored in the ECHO-g model simulations (González-Rouco et al., 2006) by comparing the model response 10 using the whole array of model grid-points and a decimated mask overlapping with areas where actual BTPs are available. Both spatial samples lead to virtually identical results, thus supporting the present distribution of BTPs as a viable sample for estimating terrestrial SAT variations. González-Rouco et al. (2006) further illustrated the potential of the borehole methodology to produce GST reconstructions using a forward 15 model (see Sect. 2.2) to simulate synthetic BTPs within the model world and subsequently inverting and averaging them to produce a hemispheric GST reconstruction that could be compared to the model simulated surface temperature.
This approach will be revisited in the remaining part of this section to illustrate how the method can be used to address the potential influence of logging dates and BTP 20 depth on borehole climate reconstructions. We will focus on North America as the example case where both variables show considerable spatial variability as depicted Fig. 5.
The original model grid (gray dots in Fig. 5a) was decimated to retain a subset of grid-points that were closer to the observational borehole network (red and green in 25 Fig. 5b, respectively).
BTPs were simulated at each location using GST anomalies at the deepest (−9.84 m) soil model level. Anomalies were calculated with respect to the full period of simulation and diffused downwards into the ground using Eq. (2) as described in These simulated profiles represent a perfect world in terms of temperature perturbation profiles for which the SAT is known. These profiles can be inverted in order to test the potential of inversion models to retrieve GST histories in a process that mimics the 5 borehole approach to climate reconstruction. Additionally, the profiles can be further degraded making them more similar to real borehole temperature vs. depth anomalies and test the effect of degrading factors on the final reconstructions. Thus, the smaller grid-point subset was further deteriorated to incorporate a more realistic representation of possible biases related to logging dates and borehole depth. 10 The distribution of logging dates in the borehole dataset is spread over the last 40 to 50 years (Fig. 5c); overall the early measurements where made in USA and the more recent logs are distributed mostly over Canada. The potential impacts of this different distribution have been discussed in the previous section. This geographical spread is replicated within the model world over an equivalent range of 40 years (Fig. 5d), 15 ending in 1990 AD in order to accommodate the length of the model simulations. Model grid point series were trimmed according to this distribution before diffusing them into the ground. In this way, the simulated BTPs mimic within the model world different measurement moments. For the control simulation these dates are arbitrary since there is no correspondence to real years and were assigned in such a way to reproduce an 20 equivalent loss of data at the end of the simulation as it was done for the forced runs.
A further degradation step consisted of shortening the simulated BTPs according to the real depth distribution shown in Fig. 5e and replicated in Fig. 5f. Though this distribution ranges over a depth of 1000 m, in practical terms only changes on the first ca. 400 m (see Fig. 6a) can have some effect, since bellow this depth the simulated 25 perturbation profiles virtually converge to zero. This approach attempts to address the impact of loosing part of the climatic history in some areas where boreholes are shallower.
The resulting sets of profiles are shown in Fig. 6a where the larger size of the original parison to the control run, the depth of warming onset is deeper in the forced than in the model simulations. In the control BTPs this shallower warming can be attributed to trends caused by internal variability at the end of the simulation (see CTRL in Fig. 1b).
The simulated BTPs were inverted using SVD  to obtain local GST histories and the mean was area-weighted over the whole domain. For each 10 local GST history, prior to the obtention of the spatial average, temperatures were kept constant since the logging date ( Fig. 5d ) until 1990 AD in order to mimic published exercises (e.g. Harris and Chapman, 2001). This is a rather conservative scenario since in the real world applications there is availability of instrumental or reanalysis data that can be used to fill the gap after the logging date in a more accurate fashion. The cutoff 15 level used for the SVD-eigenvalues in this exercise was 0.025. Figure 6b compares the low frequency evolution of North American SAT in each of the simulations and the latitude-average of the inverted GST histories, both for the complete grid (solid lines) and the perturbed subset (dashed lines); differences between both GST inversions are highlighted. 20 The three GST histories successfully recover the low-frequency changes in North American SAT. The broad differences between the control and forced simulations are recovered, including the different level of MWP warming in FOR 1 and FOR 2. The differences between the full grid-point net and the perturbed subset are negligible in the context of the low frequency variability reproduced by the inverted histories. Thus, 25 Fig. 6b lends support for the borehole reconstruction approach in recovering low frequency changes of past temperature and suggests that the potential effects of variability in logging dates and depths as simulated here are minor and will hardly explain offsets between borehole based and other proxy reconstructions.

EGU
Nevertheless, this exercise rests on some assumptions that are worth discussing. Since this approach only simulates temperature perturbation profiles, i.e. the T t (z) terms in Eq. (1), all the potential problems related to the separation of the climate transient, T t (z), from the geothermal gradient in Eq. (1), T 0 +q 0 R(z), are overlooked. Additionally, there are at least two relevant issues when separating the geothermal and 5 climatological information in BTPs: the actual depth of the profiles and the existence of noise.
The treatment given to borehole depth in Figs. 5 and 6 addresses the fact that shallower boreholes miss a part of the climate history recorded deeper into the ground but skips the problems of discriminating the geothermal from the climatological signal 10 in shallower boreholes which have not fully met the geothermal equilibrium gradient. Thus, the previous results assume that there is a perfect separation between both types of signals and addresses only potential problems related to the subsequent processing of perturbation profiles.
Borehole noise is another issue that further complicates both the inversion of pertur-15 bation profiles and the discrimination between the geothermal background and the climate transient. The latter problem could be easily addressed by adding some realistic noise to the temperature profiles and continuing with the other steps of the methodology. Therefore, an additional improvement of the approach that can be considered in future work is the simulation of a random and realistic spatial distribution of geother-20 mal gradient values and noise that will allow for a more realistic reproduction of all the steps of the procedure. However, it is not clear that the inclusion of these features will lead to large biases in the recovered GST histories since errors in separating the terms in Eq.
(1) at each BTP should randomly distribute contributing to produce under-and over-estimations of surface temperature. 25 The impact of noise in perturbation profiles is in practice overcome by considering smoother solutions to the heat diffusion equation. In the case of the SVD inversion approach used herein this implies the use of larger eigenvalue cutoff levels, thus retaining less information from the measurements. This avoids unstable solutions produced by EGU over-weighting errors in the eigen modes (see discussion in Beltrami and Mareschal, 1995;Beltrami and Bourlon, 2004). The effects of this can be observed in Fig. 6c where various GST histories corresponding to cutoff levels between 0.025 and 0.3 are shown for the case of the FOR 1 simulation. The 0.025 value is rather low and appropriate for noise-free profiles whereas values between 0.1 and 0.3 are typical of real profiles with 5 the presence of noise (Beltrami and Bourlon, 2004). The provided solutions show a progressive decrease of the differences between the original complete set of GST grid-point BTPs and the date and depth perturbed subset. Also a progressive loss of skill in depicting the warming from the LIA to the MWP is observed in addition to less cooling from the simulated 20th century to the LIA.
Interestingly, this behavior suggests that a more realistic replication of the borehole method leads to under-rather than over-estimation of past climate variability and, in particular, LIA cooling.
A last comment in this discussion concerns the assumption implicit in calculating temperature anomalies with respect to the full period of simulation before producing 15 BTPs. This apparently innocuous step is arbitrary and not free from quite strong assumptions since the selection of a reference period, and thus a mean temperature reference value for the calculation of anomalies produces important changes in the simulation of BTPs. However, the selection of a particular reference is irrelevant if the analysis rests on taking for granted that a perfect discrimination between the geother-20 mal gradient and the climate transient is possible. Under this assumption, the analysis starts by adopting the simulated perturbation BTPs as valid and addressing all subsequent aspects of the methodology. Decisions concerning the selection of a mean temperature reference level to calculate anomalies bear more importance in the comparison of observational and simulated BTPs. Therefore, this will be further discussed 25 within the context of the next section. 26 between GCM simulations and observed BTPs; 2) the level of fidelity with which GCMs can recreate the details of the air-ground temperature coupling at local scales. The first point targets the comparison of model simulations and borehole paleoclimatic reconstructions or BTPs. The second question attempts to address whether the degree of realism of the simulated land surface processes is adequate both in the paleoclimatic 10 and in the future climate change context.

Simulated and observed borehole temperature profiles
Reliable projections of future change require that climate models prove effective in simulating past changes in broad agreement with evidence provided by climate reconstructions. Advances in the convergence between simulation and reconstruction 15 approaches will require reducing various types of uncertainties on both sides (Jansen et al., 2007;Hegerl et al., 2007b). Joint analysis of climate reconstructions, model simulations and estimated past external forcing can illustrate the level of agreement between model results and reconstructions (Jones and Mann, 2004;Zorita et al., 2004;Moberg et al., 2005), but also to try to constrain the range of climate sensitivity (Hegerl 20 et al., 2006) and to attribute past changes registered in millennial climate reconstructions to particular external forcings (Crowley, 2000;Bauer et al., 2003;Hegerl et al., 2003).
In the context of borehole climate reconstructions, first steps have been taken in designing a means of comparison of the information provided by BTPs with model EGU averages of ECHO-g simulated SATs were used to construct the subsurface thermal profiles as described in Sect. 3.2. These were subsequently compared with average BTPs measured at each region. It was found that in all cases subsurface anomalies so calculated from the forced simulations were in better agreement with observed profiles than the synthetic profiles derived from the control simulation. 5 Stevens et al. (2008) arrived at the same conclusion in an extension of the analysis to the whole of North America in which alternative ways of comparing real and simulated profiles are explored. Both studies target a broad evaluation of agreement or disagreement between observed BTPs and model simulations, and in their conclusions they qualitatively support the sensitivity of BTPs to external forcing and thus their use-10 fulness for climate change assessment; however these works cannot be included in the context of climate attribution since the model simulations involved cannot discriminate among the influence of various natural and anthropogenic forcings. Stevens et al. (2008) and Beltrami et al. (2006b) also highlighted a seemingly qualitative agreement in the east-west arrangement of trend in the observed and in the 15 synthetic boreholes from the forced simulation over northern North America. Lower (larger) amplitude of warming has been reported for the western (central and eastern) parts of Canada (Majorowicz et al., 2002) and, in particular, changes from the western to the eastern side of the Cordillera have been noted (Wang et al., 1994;Pollack and Huang, 2000). Figure 7 provides a spatial perspective into the trends simulated by the 20 ECHO-g model for the period 1700 to 1990 AD. In spite of different internal non-forced variability in each simulation, both model integrations produce a remarkably similar pattern with smaller trends over western Canada and larger trends in the central and eastern territories. The large-scale temperature response in these model simulations is shown in Zorita et al. (2005) and can be described as a land-ocean thermal contrast. 25 This pattern is superimposed to regional trends that are associated with changes in the intensity of the main modes of atmospheric circulation.
Over North America most of the main land and the regions usually covered by perennial snow the model simulates larger warming in the continent with diminishing ampli-28 EGU tude towards the ocean and with a regional minimum in western Canada that would be consistent with borehole observations. This match may well be coincidental and it is premature to claim that the observed warming pattern can be attributed entirely to the response to external forcing. However, given that the pattern seems to be robust in both ECHO-g simulations, it is interest-5 ing to further explore some broad features of agreement and disagreement between Fig. 7 and observational evidence. In addition to the east-west arrangement of amplitude trends over Canada, the simulated warming is largest over the northern areas of Alaska, also consistent with borehole observations (Lachenbruch and Marshall, 1986).
The western USA displays however, trends that are comparable to those in central 10 and eastern Canada, a feature that contradicts the experimental evidence for the western USA (e.g. Harris and Chapman, 1995) that has shown less warming than central and eastern Canada. It can be speculated that most of the boreholes over this area correspond to early logs in the dataset mostly from the 1960's and 1970's (see Fig. 5c), that fail to register the posterior warming: Stevens et al. (2008) argue that these smaller 15 trends are consistent with meteorological observations which show some decadal cooling before the average logging dates. Thus, if the pattern in Fig. 7  EGU some assumptions to calculate model SAT or GST anomalies which are subsequently forward-diffused into the ground. Usually this has been done by selecting a temporal reference period so that anomalies are obtained with respect to the mean temperature during this period. Such a selection is not evident both because boreholes integrate the influence of climatic events 5 before the start of the 1000 year ECHO-g millennial simulations and also because these simulations do not necessarily yield a deteailed representation of past temperatures (Raisanen, 2007). In fact, Figs. 6 and 7 show that different model simulations provide climate realizations that are not identical and that ultimately translate into slightly different synthetic BTPs. 10 Beltrami et al. (2006b) took a simplified approach and selected the mean of the whole period of simulation (1000 to 1990 AD) as the reference to calculate anomalies. However, this choice is as arbitrary as any other selection since it turns out that the choice of the reference period, and thus the reference mean temperature, influences the shape of the forward simulated BTPs. Figure 8 shows an example of the impact 15 on the shape of diffused profiles for a set of selected reference periods using the SAT series from FOR 2 at an arbitrary grid point in central North America (50.1 • N, 97.5 • W) as a boundary condition. In the various profiles, anomalies are calculated with respect to the different means of the selected reference period.
The problem illustrated in Fig. 8 can be related to that of comparing instrumental time 20 series and observed BTPs (Harris and Gosnold, 1999;Harris and Chapman, 2001;Harris, 2007). Also in this case the selection of a mean temperature reference is not self-evident because boreholes integrate a much larger climate history that the approximately two-century-long instrumental period or even the 1000 year GCM paleoclimatic simulations considered here. 25 This makes the climate signal recorded in boreholes and that diffused into synthetic BTPs not fully comparable since the existence of an influence corresponding to longer time scales (i.e. multi millennial or glacial to inter-glacial) cannot be ruled out in experimental BTPs of the typical depths involved herein.

EGU
The contribution of this background or pre-observational climate has been taken into account in exercises of comparisons of real BTPs with instrumental data making use of the pre-observational mean concept (see Harris, 2007, and references therein).
From this perspective the selection of an elusive reference period respect to which anomalies should be calculated is avoided and substituted by a search for 5 an average temperature reference level that represents the contribution of the preobservational/instrumental local climate to the observed subsurface thermal regime.
The use of such an approach in the case of instrumental time series and borehole temperature profile comparisons is founded in the assumption that GSTs and SATs perfectly track each other. The search for the appropriate pre-observational reference 10 temperature is performed through minimizing a measure of similarity between the diffused instrumental time series of temperature anomalies and the observed borehole temperature profiles. The selected value not only provides minimal error by definition, but this error is also expected to be very small in magnitude due to the assumed GST and SAT coupling. In the case of the comparison between observed and syn-15 thetic boreholes this rationale cannot be invoked since millennial long GCM simulations cannot be regarded as a detailed representation of true past SAT recorded in real boreholes. Thus, the strategy of searching for a minimal separation between real and simulated borehole temperature profiles cannot be designed under the rationale that the GCM-simulated thermal regimes will necessarily match the observed ones as in 20 the comparisons between real boreholes and instrumental data. Yet, in the absence of other reasonable approaches this is still a useful guide that allows for scanning over a wide range of possible mean temperature reference levels respect to which simulated anomalies can be calculated and forward diffused into the ground.
In this fashion, a family of curves can be generated for each selected region by 25 forcing the forward model by simulated SATs series of anomalies considering a range of mean temperature reference levels. This variability in the synthetic BTPs illustrates the uncertainty in the selection of a temporal reference period, or associated mean temperature reference level, in each model simulation. EGU and observed BTPs can be quantified through various criteria which can serve as a means of model-data comparison. This strategy was used in Stevens et al. (2008) using three measures of distance between simulated and real BTPs: the root mean square error between pairs of synthetic and observed profiles, the depth at which the LIA cooling is recorded in BTPs and the amplitude of temperature change from this 5 depth to the surface. Results showed that the forced simulations presented depths of pre-industrial cooling and rates of surface warming which were in better agreement with observed BTPs than the control simulation. This strategy illustrates a possible probabilistic approach that allows to advance in the comparison of observed and synthetic BTPs in spite of the uncertainties re-10 lated to the choice of a reference thermal state that accounts for the influence of preobservational or pre-modeled climate.

A glimpse at the surface geothermal climate
The comparisons of simulated and observed BTPs presented in the previous section aim at finding consistency between model simulations and past climate as registered in 15 deep borehole temperature logs. Such assessment ultimately targets questions within the long term climate change context since only the lowest frequencies of GST are retained after being diffused into the ground.
A central aspect in the climatological interpretation of borehole profiles and also of time series of subsurface temperatures are the physical and biological mechanisms 20 that may modulate the transfer of heat from the atmosphere to the upper layers of the soil. Therefore, a complementary view at the geothermal regime can be attained by considering the evolution of ground temperatures near the surface as represented by time series of soil temperature monitored at various levels, typically to a maximum depth of about 10 m. Understanding spatial variability in ground temperatures and the 25 relation between GST and SAT evolutions at all time scales is relevant both in the context of borehole climatology and also for climate variability and change simulation purposes. 32 how atmospheric and surface conditions perturb the relation between SAT and GST; ii) to verify that the propagation of surface temperature changes to the subsurface is dominated by heat conduction and to identify and understand potential situations that deviate from this basic assumption. Therefore, a large body of research has been devoted to understand the processes modulating the generation and downward propa-10 gation of a climate signal at the ground surface. Some of those mechanisms, like snow cover and soil moisture changes, have been briefly discussed in Sect. 3; others involve vegetation and land use changes, evapotranspiration, precipitation, albedo, snow melting, etc. Most of the studies have a local focus which allows for the consideration of specific conditions, climatological types, seasons or soil types. Through a large variety of situations and case studies, the broad validity of the conductive regime and the coupling of air-ground temperatures at frequencies beyond the annual cycle is supported both by observational evidence (e.g. Putnam and Chapman, 1996;Beltrami, 2001a;Sokratov and Barry, 2002;Smerdon et al., 2003Smerdon et al., , 2004Hu and Feng, 2005;Mottaghy and Rath, 2006;Chudinova et al., 2006;Smerdon et al., 2006) and land surface model 20 assessments (e.g. Lin et al., 2003;Stieglitz et al., 2003;Bartlett et al., 2004;Pollack et al., 2005;Stieglitz and Smerdon, 2007). In turn, GCM modeling efforts are concerned with the degree of realism with which the surface thermal regime is reproduced since soil temperatures and with a broader perspective, soil moisture and other hydrological variables, are relevant parameters in 25 air ground interaction. As discussed in Sect. 3 the soil thermal and hydrological state are shaped by surface conditions (surface temperature, precipitation, snow cover, vegetation, land use, etc., Lewis and Wang, 1998;Hu and Feng, 2003;Zhang, 2005;Zhang et al., 2005;Davin et al., 2007) though also, the thermal and hydrological con-Introduction EGU ditions in the subsurface have the potential to influence surface climate by modulating turbulent latent and sensible heat fluxes (e. g. Peters-Lidard et al., 1998;Sokratov and Barry, 2002;Schaefer et al., 2007) that feed back to regional climate (e.g. Bhatta et al., 2003;Kueppers et al., 2007;Taylor et al., 2007). In addition, the exchange of matter fluxes (e.g., CO 2 , CH 4 ) regulated by soil temperatures and moisture plays a critical 5 role in global climate through the dynamics of biogeochemical cycles (e.g. Risk et al., 2002a,b;Kellman et al., 2007;Riveros-Iregui et al., 2007). Therefore, a realistic simulation of air-ground interactions is crucial to obtain a credible representation of the energy balance at the surface and related climate feedbacks on temperature, precipitation, evapotranspiration, convection, regional circulation, etc.
10 (e.g. Walker and Rowntree, 1977;Dirmeyer, 2000;Zhu and Liang, 2005;Seneviratne et al., 2006;Davin et al., 2007;Fischer et al., 2007b;Miguez-Macho et al., 2005. The work of Zhu and Liang (2005) constitutes an interesting contribution that illustrates the performance of GCM and regional climate models (RCMs) in reproducing soil temperatures.  Hu and Feng (2003) dataset over the 20 contiguous U.S. Overall, the spatial variability, annual cycles and interannual variability simulated by the R 2 and CMM 5 capture realistically the details of soil temperature variability. However some systematic cold biases in the simulated soil temperatures are found for reasons still unknown. The higher resolution CMM 5 simulation produces more realistic regional details and overall smaller biases than the driving R 2 GCM 25 reanalysis output, thus supporting the added value of the dynamical downscaling.
The coarse resolution of GCMs and their relatively simple parameterizations (von Storch, 1995) hamper an accurate representation of specific soil thermal regimes as in the case of permafrost. At high latitudes, permafrost soils and the variability through downscaling approaches are ussually adopted: the use an off-line soil model forced with boundary conditions provided by either a GCM (e.g. Stendel and Christensen, 2002;Sazonova and Romanovsky, 2004) reanalysis data (Malevsky-Malevich et al., 2001;Nicolsky et al., 2007), or observations Moelders and Romanovsky, 2006); or alternatively, the use a nested RCM with a more sophisticated 10 land surface component (e.g. Sushama et al., 2007). The wealth of studies in this branch of research is considerable and deserves attention by itself in its applications both to current and future climate. However, few pieces of work actually comment on the performance of the forcing climate models to simulate soil temperatures. This is due to the general inability of the forcing GCMs to reproduce permafrost and to the fact 15 that this task is ultimately accomplished by the off-line land surface model or the nested regional model; also, often other permafrost related variables are diagnosed and compared to observations for validation like the thickness of the active later, permafrost area, freezing and thawing indices, snow cover, etc. (Malevsky-Malevich et al., 2001;Sazonova and Romanovsky, 2004). Some interesting case studies exist though that 20 describe comparisons of simulated and observed soil temperatures in the process of their validation assessments. For instance Malevsky-Malevich et al. (2001) report on a fairly good agreement between observed soil temperature at 15 m depth over Siberian permafrost and that simulated by their land surface model driven by reanalysis data (Kalnay, 1996). Sushama et al. (2006) provide an example in which they validate the EGU subsurface in winter and spring and warm biases in summer and autumn, with the annual contributions balancing out and presenting reasonable agreement at the deeper layers of the simulated and observed boreholes.
In the remaining part of this section this text will focus on the broad aspects that characterize the conductive regime in the upper 10 m of the soil as simulated by the 5-5 layer soil submodel embedded in the ECHO-g model. This will help both the purpose of providing a qualitative description of the realism of the model in non-permafrost areas and at the same time to introduce the problem of the boundary condition that is typically prescribed at the bottom of the soil sub-model in GCMs. Figure 9 shows the evolution of daily temperatures at two locations representative 10 of a positive and a negative SAT-GST difference in Fig. 1d. In one of them, air-ground coupling is strongly modulated by soil moisture whereas in the other the presence of snow cover is more important as discussed in Figs. 2 and 3. The grid-point selection is not fully arbitrary though since it is intended to serve for comparison with observational data described in Smerdon et al. (2003Smerdon et al. ( , 2004Smerdon et al. ( , 2006. Two model grid-points were se-15 lected close to Fargo and Cape Hatteras. For both cases, the last 6 years in the FOR 2 model simulation are shown. Fargo is a good example of a snow variability dominated site with SAT values which often range between −30 • C in winter to 30 • C in summer (e.g. Smerdon et al., 2004). In the simulated years shown herein the range is slightly smaller, with minimum (maximum) values of −25 • C (25 • C) reached during some win-20 ters (summers). The presence of snow hinders the cooling of the upper subsurface layers, so that the annual cycle is shifted with depth to above zero temperatures producing the typical negative SAT-GST offset (see Fig. 1): at 1.23 m bellow the ground temperatures seldom reach values bellow −5 • C and at 4.13 m they are always positive. Cape Hatteras is an example of site where soil moisture is important. The SAT annual 25 cycle typically ranges between 0 • C and 30 • C in the observations and 0 • C and 25 • C in the simulations. In this case there is no apparent shift with depth since the temperature series scarcely reaches negative values.
The annual cycle is further damped and phase shifted with depth due to the effects 36 EGU of conduction on the downward propagation of a harmonic temperature signal. This becomes apparent at both sites. For each harmonic component of the temperature variations at the surface with amplitude A 0 , harmonic frequency ω and phase ǫ represented by T (z=0, t)= A 0 cos(ωt+ǫ), the conductive propagation with depth of this signal is given (Carslaw and Jaeger, 1959) for depth z and time t by where the wave vector k is π/P κ and the wavelength λ=(4πP κ) 1/2 , P being the harmonic period and κ the thermal diffusivity. Thus, for each frequency component, the amplitude (A) and phase shift (φ) change with depth according to: The amplitude decay and phase shift due to heat conduction is better illustrated in Fig. 10 by focusing exclusively on the annual cycle. A cosine wave of 12 month period has been fitted to the annual temperature evolution at each soil model level to determine the amplitude and relative phase of the wave. In order to consider a larger 15 number of years and attain more robust estimates, the period 1960 to 1990 AD was used. The natural logarithm of the amplitude and the phase shift has been represented for the grid-points close to Fargo and Cape Hatteras over the complete width of each soil model layer (solid lines). Also shown is the expected profile of a conductive regime using the model thermal 20 diffusivity value for the grid-points close to Fargo and Cape Hatteras (dashed lines). The modeled range of surface temperature at Fargo is larger than at Cape Hatteras as in the observations (Smerdon et al., 2004). The average fitted amplitude of the SAT annual wave at Fargo is 13.9 K whereas a value of 17.8 is found in observations along the period 1980 to 1999 AD; for Cape Hatteras the model simulation provides a value EGU the tendency to underestimate the amplitude of SAT annual range would be consistent with the underestimations reported in Zhu and Liang (2005) for the same area. The vertical gradients of amplitude attenuation and phase shift are somewhat intermediate in the model simulations compared to the observed ones at Fargo and Cape Hatteras. This can be justified by the model thermal diffusivity (κ=0.75×10 −6 m 2 s −1 ) 5 being larger than that deduced from observations at Fargo (0.37×10 −6 m 2 s −1 ) and slightly smaller than that at Cape Hatteras (1.04×10 −6 m 2 s −1 ). The evolution of phase shift with depth basically overlaps at both sites. Both changes in amplitude and phase for each soil model level follow quite closely the conductive regime. However, for the lowest model layer slight differences from the conductive regime 10 should be expected, because of the imposed boundary condition of zero flux at the bottom soil model level. Such deviations have been reported by Lynch-Stieglitz (1994) and Sun and Zhang (2004) for the annual cycle and extended to lower frequencies by Smerdon and Stieglitz (2006). The influence of a zero flux lower boundary condition on the subsurface propagating signals at a given depth depends on the frequency of Zhang (2004) extended for various BBCP depths. Smerdon and Stieglitz (2006) assessed the impact of the BBCP on the propagation of harmonic signals of daily to millennial timescales using analytical solutions of the one-dimensional heat conduction equation. The results indicated that the appropriate BBCP is dependent on the time-scale of interest and this should be taken into consideration in soil models within 25 GCMs. Given that typical depths for soil models in GCMs range between 3 and 10 m, it is likely that significant deviations occur in the subsurface simulated thermal regime due to these shallow configurations. In the case of the ECHO-g model, the BBCP at EGU 9.83 m is deep enough to produce only negligible perturbations on the annual wave. However, perturbations at lower frequencies cannot be ruled out. The amplitude damping of frequencies shorter than the annual cycle is illustrated in Fig. 11 by showing a low pass filtered and extended version of the time series in Fig. 9. The annual cycle has been removed by subtracting from each monthly value 5 the corresponding monthly long term average for the period 1960-1990 AD. The series of anomalies for SAT and soil temperatures at each depth reveal now more clearly the low frequency variability. It is apparent that the high frequency variations in the layers closer to the surface are quickly damped and almost only decadal variability is present at the lowest model level. In the case of Fargo the large negative anomalies in SAT are visibly reduced to smaller changes in temperature at the first model layer, thus showing the effect of snow insulation. In the case of Cape Hatteras this effect is absent and the ground temperature anomalies of the upper layers overlay over those of SAT. A face shift with depth for the envelope of the lowest frequencies is also observable at both sites. 15 The amplitude damping and phase shift for frequencies bellow the annual cycle can also be illustrated using spectral analysis. The amplitude of each component harmonic at the top and bottom layer can be obtained using a spectral analysis for which ordinates have been transformed to represent the amplitude of the wave at each frequency. The ratio of spectral amplitudes for the top and bottom time series is an estimation of 20 the amplitude damping at each frequency. In addition, phase cross spectra can be used to picture the relative phase shift for the component harmonics between the top and the bottom layers. Figure 12 shows this analysis at the two selected grid-points together with the values that would be expected from the conductive regime (grey dashed lines) obtained from Eq. (3). The amplitude attenuation is in good agreement with changes 25 expected from the conductive regime with no zero flux condition for high and intermediate frequencies up to about 0.17 yr −1 (ca. 6 years). Bellow that limit a gradual deviation from the pure conductive regime is observed. The phase shift is in good agreement with the conductive regime up to very low frequencies where some slight deviation EGU is perceptible and zero phase shift is reached around 0.025 yr −1 . Thus, the low frequency time-scales are not attenuated as much as expected and become in phase after a certain threshold frequency. For instance, frequencies of 0.1 yr −1 , while phase shifted according to the theory, are only attenuated to 62% of their surface value when they should be reduced to about 42%; frequencies around 0.05 yr −1 get attenuated to 5 about 90% of their value at the surface when they should be reduced to about 55% of it. This indicates that in Fig. 1d the filtered time series at various depths should present a gradual attenuation with depth rather than an almost parallel evolution. The considerable deviations from the infinite half space conductive regime shown in the amplitude damping in Fig. 12 are due to the influence of the BBCP used in 10 the model. This can be easily shown by considering that the amplitude attenuation produced in such scenario is given by (e. g. Smerdon and Stieglitz, 2006): where d is de depth of the BBCP. Equation (4) can be valued for for each frequency thus obtaining the expected amplitude attenuation for the specific case of the ECHO-g 15 model. This is shown with a blue dashed line in Fig. 12. For the phase shift a similar analysis can be performed (not shown here). Thus, as shown by other authors with one dimensional head conduction models (Lynch-Stieglitz, 1994;Sun and Zhang, 2004;Smerdon and Stieglitz, 2006), the use of zero flux conditions and shallow soil models in GCMs can impact the simulation of 20 temperature at low frequencies with a sensitivity that depends on the frequency of the signal, the value of the subsurface thermal properties and the BBCP. In addition to the temperature simulation, it can be argued that the use of shallow soil models impacts other variables like the simulated amount of heat stored into the subsurface. This can be particularly important in the context of simulations of future climate change. and produce other potential impacts on infrastructures, ecosystem changes, carbon storage, etc (Cherkauer and Lettenmaier, 1999;Oechel et al., 2000;Nelson, 2003;Stieglitz et al., 2003). Therefore, it is not extrange that most efforts concerning the simulation of subsurface temperatures in climate change scenarios have been developed with a focus on permafrost regions. 10 Sushama et al. (200610 Sushama et al. ( ) simulated present (196110 Sushama et al. ( -1990 and changes in future (2041-2070) soil temperatures over a domain covering northeastern Canada using a 1-D heat conduction model (Goodrich, 1982) driven by surface boundary conditions provided by skin surface temperatures and snow cover from transient climate change simulations with the Canadian Regional Climate Model (CRCM, Caya and Laprise, 1999) at a reso-15 lution of 45 km latxlon, in turn, forced by the Canadian Coupled GCM (CCGCM2). Their results indicate future warming under the IPCC scenario ISP2a (Leggett et al., 1992) in annual mean, maximal and minimal soil temperatures near the surface and a widening of the active layer above the permafrost by more than 50%. In Sushama et al. (2007) they use an alternative approach by avoiding the use of the off-line soil model and 20 implementing a dynamical downscaling using an improved version of the CRCM that incorporates sophisticated land surface model (CLASS Verseghy et al., 1993). They indicate increases of 4 to 6 • C for the period 2041-2070 relative to 1961-1990 under climate change scenario A 2 with implications also for soil hydrology.
Though some healthy discussion exists, (Delisle, 2007), comparable results have 25 been recently reported for various permafrost areas which add up to a wealth of existing research (see Nelson, 1996, 1997;Malevsky-Malevich et al., 2001;Nelson and Anisimov, 1993;Lawrence and Slater, 2005, and  A detailed account of the results of studies assessing permafrost sensitivity to climate change with current climate and land-surface models is beyond the possibilities 10 of text. Achievements in this area of research have brought to the forefront the necessity of a more realistic simulation of the subsurface climate to adequately reproduce air-ground interaction; this being also valid for non-permafrost areas (e.g. Liang et al., 2001;Seneviratne et al., 2006). In addition, a complementary perspective that addresses the importance of an improved representation of subsurface climate in GCMs 15 is also pertinent since it could have implications not only in the simulation of regional to large scale soil regimes like in the case of permafrost but also in climate at hemispherical to global scales through changes in the energy balance at the surface.
GCMs and RCMs incorporate relatively shallow (between 3 and 10 m deep) soil models with a flux condition (typically zero flux) at the bottom boundary (see Sect. 4.2).
Since such BBCP affects the simulation of subsurface temperatures near the surface (Smerdon and Stieglitz, 2006), it can be argued that the surface could be possibly influenced by this artificial BBCP, for instance by changes in surface heat fluxes, variations in the permafrost, hydrology, ground-air coupling, etc. Improved land-air interactions have proven to have a strong impact in the simulation of future climate (e.g. Senevi-25 ratne et al., 2006;Fischer et al., 2007a); the question arises whether a more realistic representation of the underground thermal regime can have potential impacts on surface climate. One source of influence could be the heat storage capacity of the ground. The shallow soil model configurations can potentially under-estimate the heat storage EGU capacity of the ground and concentrate the simulated warmth near the surface instead of diffusing it deep into the ground. This additional energy which is not stored into the ground can be compared to a change in external forcing in as much as it represents a surplus of energy that is misplaced into other parts of the climate system. It can be argued that as an additional 5 amount of energy it would contribute to warming in the global energy balance (Stevens et al., 2007).
If that were the case, part of the energy which in reality should be stored into a realistically deep subsurface would be available, in the model, for other parts of the climate system and potentially produce a warmer climate at the surface. The capacity of the 10 subsurface to store heat as a function of BBCP was specifically revised by Stevens et al. (2007) who forced a land surface model (Goodrich, 1982) with SAT provided by the ECHO-g simulations presented herein. One of the simulations, FOR1, was continued to 2100 AD under IPCC A 2 and B 2 scenarios  which allowed for an assessment of the amount of heat stored into the soil under future climate 15 change conditions. This assessment was done by initializing an off-line land surface model with a configuration of 1000 m depth with the 1000 to 1990 AD temperature simulations; after that, several sensitivity experiments were produced by driving the land surface model from the initial conditions attained at year 1990 AD with the temperature changes simulated under the A 2 and B 2 scenarios and for different configurations in 20 which the BBCP was imposed at depths between 1 and 1000 m. Results suggested that shallow BBCPs (10 m) could reduce the capacity of the continental subsurface to store heat by about 1.0×10 23 J during an A 2 scenario run. Deepening the BBCP to about 120 m increased the heat storage capacity by about 5 times. This increment is more than an order of magnitude greater than the heat absorbed by both the whole 25 atmosphere and the continents in the second half of the 20th century (Beltrami et al., 2002;Levitus et al., 2005;Beltrami et al., 2006a;Huang, 2006a,b) Figure 13 shows the heat gained by the ECHO-g soil model along the simulations used in Stevens et al. (2007). The spatial distribution of total heat stored into the ground EGU at the end of the simulations is shown in Fig. 13b. At 2100 AD 1.3×10 8 J (0.8×10 8 J) are gained by the 9.83 m deep subsurface in the A 2 (B 2) scenario. These quantities are compatible with those calculated in Stevens et al. (2007) by forcing the land surface model using a depth of 10 m with FOR 1 NH SAT temperatures. If a 150 m deep soil model would be able to store as much as 5 times the original quantity this would mean 5 that about 6.5×10 8 J would be gained if a more realistic subsurface was used. This quantity is comparatively larger than the actual difference of heat accumulated within the A 2 and B 2 scenarios: 0.5×10 8 J. This suggests that changing to a more realistic (deeper) subsurface could have a larger impact on heat storage and on the related energy balance than using a different forcing scenario.

Conclusions
The last years have witnessed a growth of reconstruction and model simulation studies targeting an improved understanding of climate variability during the last millennium. Both types of approaches offer multiple facies of potential interaction that can contribute to our knowledge of past climate through comparison and validation of various 15 aspects of reconstructions and simulations and can also have implications within the context of future climate change assessment, after all a major motivation for paleoclimate research.
It is reasonable to think that research blending information and knowledge from reconstruction and simulation experiments should be designed from an ad hoc perspec-20 tive that adapts to the specific peculiarities of climate reconstructions regarding limitations in methodologies or biases to represent the climate signature of certain seasons, regions or time scales of variability, as well as those related to climate modeling regarding resolution, parametrizations, uncertainties in forcing estimates, etc.
Borehole climatology has significantly contributed to this context by recovering infor-25 mation of past multi-centennial surface temperature trends from the thermal signature imparted in subsurface rocks. This text has reported and discussed about recent steps 44 EGU and ideas that stem from the interaction of this discipline with climate model exercises.
We have focused mostly on continental boreholes allowing very limited or none attention to other subsurface types for which a more extensive treatment would be desirable elsewhere. One of them is permafrost which receives partial attention in Sects. 4.2 and 5. In addition, there are also important pieces of work which tell interesting lessons 5 about boreholes measured in ice cap areas. These boreholes deliver longer temperature reconstructions and some interesting attempts have been made to understand their thermodynamics with the aid of GCMs (Dahl-Jensen et al., 1998;Jouzel, 1999;Krinner et al., 1997;Werner et al., 2000) The use of GCMs as a test bed for climate reconstruction methods has provided 10 evidence for the robustness of the borehole approach in delivering a consistent picture of low frequency climate variability under various methodological issues. In particular, trends in simulated snow cover and soil moisture do not seem to hamper the ability of inversion models to retrieve the past temperature evolution; the same seems to apply with regard to the potential influence of experimental aspects like the distribution 15 of logging dates and depth. Results suggest that there could be an interesting issue in the effects of noise, which could in extreme cases lead to an under-estimation rather than an over-estimation of past climate variability as it has been suggested in the past literature. Future efforts could consider a more realistic involvement of noise, both in borehole profiles and also with a regional perspective incorporating non-climate 20 perturbations like horizontal heat advection, changes in vegetation, etc. The effect of multi-centennial changes in land surface cover has not been considered so far. Comparison of observed and synthetic BTPs has given its first steps. Results point out a sensitivity of the temperature signature in BTPs to trends in external forcing and that there might be some regional coherence between observed and simulated BTPs 25 which deserves further exploration. A relevant issue that stems from this line of work is that uncertainties in the generation of synthetic BTPs lead to the necessity of developing comparison strategies in a probabilistic framework. An additional issue which is lacking in this context is the incorporation of arguments that involve the effects of past EGU glacial-interglacial changes. Assessments of the realism of models in reproducing the subsurface geothermal regime evidence the need for more sophisticated soil models that help to better reproduce features of the air-ground interaction. This has been shown to be of particular importance in the simulation of permafrost regimes in which local and regional detail 5 is achieved either with nested RCMs or with off-line land surface models. The performance of GCMs in reproducing the basics of heat conduction is hampered by the placement of a zero flux bottom boundary condition in the relatively shallow soil model components. Results of different assessments point out the necessity of incorporating soil components with causally detached bottom boundary conditions either by chang- 10 ing the nature of the bottom boundary condition itself and allowing for the dissipation of heat at this end or preferably through increasing the depth of the soil model component. A further aspect of relevance concerning the depth of soil models within GCMs is the capacity of heat storage which increases with depth and may influence significantly the energy balance at the surface in simulations of future climate change scenarios. 15 This last issue constitutes a distinct example of how studies at the interface of proxy reconstructions and model simulations can provide potentially useful information to improve future climate change scenario simulations or at least to understand better the uncertainties associated to them. 007-9276-x, 2007