Late Pleistocene–holocene Ground Surface Heat Flux Changes Reconstructed from Borehole Temperature Data (the Urals, Russia)

We use geothermal reconstruction of the ground surface temperature (GST) history early obtained in the Middle Urals to determine the surface heat flux (SHF) history over the past 35 kyr. A new algorithm of GST–SHF transformation was applied to solve this problem. The timescale of geothermal reconstructions has been corrected by comparing the estimated heat flux and annual insolation at the latitude of 60 • N. The consistency of SHF and insolation changes on the interval 35–6 kyr BP with the linear correlation coefficient R = 0.99 points to orbital factors as the main cause of climatic changes during the Pleistocene–Holocene transition. The amplitude of SHF variations is about 1.3 % of the in-solation change amplitude. The increase of carbon dioxide concentrations lagged by 2–3 kyr from the SHF increase and occurred synchronously with GST changes.


Introduction
The role of orbital factors in Pleistocene climatic variations has been studied for more than 100 years since the works of Joseph Adhemar, James Croll and Milutin Milankovitch.A popular approach is comparing palaeotemperatures reconstructed from proxy data (oxygen isotopes, palynological or others) with theoretically calculated insolation.Some investigators (Peixóto and Oort, 1984;Pielke, 2003;Douglass and Knox, 2012) criticized this approach.They noted that temperature field is not an optimal parameter for climate attribution, particularly for evaluation of climatic reaction on the external radiative forcing.There is a lag between external radiative flux and temperature changes, which disappears if we consider the heat content or the surface heat flux changes.The advantage of heat flux estimation over a temperature estimation was not realized in full until now.Wang and Bras (1999) proposed the integral relation to estimate surface heat flux (SHF) changes from ground surface temperature (GST) variations.A finite-difference approximation of the relation between the GST (represented by a piecewise linear function of temperature) and the SHF was proposed by Beltrami et al. (2002).SHF history reconstructions based on borehole temperature data were made on timescales from several centuries to a millennium (Beltrami et al., 2002(Beltrami et al., , 2006;;Huang, 2006).Another approach was used in Majorowicz et al. (2012).Subsurface temperatures were calculated from solar irradiance change using information about climate sensitivity.
In the paper we use geothermal reconstruction of GST history early obtained in the Middle Urals (Demezhko and Shchapov, 2001) to determine the SHF history over the past 35 kyr.The recently developed improved algorithm of GST-SHF transformation (Gornostaeva, 2014) was applied to estimate the SHF history.Palaeoclimatic interpretation of GST and SHF histories is based on a comparison of them with orbitally driven solar insolation and atmospheric carbon dioxide changes.

The method
The GST-SHF transformation algorithm is based on the relation between surface heat flux and surface temperature changes according to Fourier's equation in one dimension: Published by Copernicus Publications on behalf of the European Geosciences Union.648 D. Y. Demezhko and A. A. Gornostaeva: Late Pleistocene-Holocene ground surface heat flux changes where q is SHF (W × m −2 ), t is time, λ is thermal conductivity (W m −1 K −1 ), and T (z, t) is temperature anomaly ( • C) at a depth z.
If GST is represented by the expression (Carslaw and Jaeger, 1959;Lachenbruch et al., 1982) where D is a constant, n is positive integer (or 0) determining the shape of temperature changes, the transient temperature anomaly at any depth is where a = λ ρC is thermal diffusivity (m 2 s −1 ), ρ is density (kg m −3 ), C is specific heat capacity (J kg −1 • C −1 ), i n erfc(α) is the nth repeated integral of the error function of α, and (β) is gamma function of argument β.Differentiation of Eq. (3) yields SHF q(0, t) = Note that the ratio E = λ/(a) −1/2 represents the rock's thermal effusivity (J m −2 K −1 s −1/2 ), or thermal inertia, characterizing the rate of heat exchange at the surface.We approximate GST history by a sum of temperature changes corresponding to Eq. (2): where i and j are positive integers related to real time by t = i × t, t = j × t , for a constant time interval t.For each addend of this sum, Using a recurrence relation the D i can be estimated for each interval of the temperature curve.The corresponding instantaneous heat flux values can then be calculated as The so obtained SHF history reconstruction will be more accurate if we calculate the average value of heat flux on the interval and refer it to the midpoint of the interval (i − 0.5) The GST-SHF transformation algorithm was tested by applying it to a harmonic function of surface temperature change with amplitude A, frequency ω, and initial phase φ: The propagation of temperature waves in a homogeneous half-space with thermal diffusivity a is described by the expression Differentiating Eq. ( 11) with respect to z, we find the ground surface heat flux change q(0,t): The relationship between the amplitudes of GST and SHF changes is determined by thermal effusivity E and frequency ω.The heat flux changes are ahead of temperature changes by π/4ω, i.e. one-eighth of the oscillation period.
The relative error of SHF estimation was calculated as the ratio of the standard error of the SHF estimation to the real amplitude of SHF variations.The test showed that approximation of temperature history by Eq. ( 5) with n = 2, 3 provides the most accurate results (Fig. 1).When GST discretization is 6 points per period, we obtain the relative error of SHF history estimation equal to 3 %, and given 10 points per period the relative error is less than 1 % (Gornostaeva, 2014).For comparison, the algorithm proposed by Beltrami et al. (2002) under the same discretization conditions provides relative errors equal to 8 and 3.5 %, respectively.

GST data and SHF estimation
As an input to our SHF reconstruction, we used the GSTH obtained from the Urals superdeep borehole SG-4 (58 • 24 N, 59 • 44 E, Middle Urals, Russia) which was presented in an earlier contribution by Demezhko and Shchapov (2001), Fig. 2. We analysed only the last 35 kyr of the GST history for the SHF reconstruction, while the paper mentioned above presents 80 kyr temperature history.Because of the decrease of the GSTH resolution with time, the interval from 35 to 80 kyr BP does not contain any noticeable GST variations.The SHF may be considered as a constant on this time interval.The reconstruction of the surface heat flux history was conducted using the algorithm described above with n = 3 (see Fig. 2).GST and SHF curves are different in shape.The temperature increase started about 15 kyr BP and after a short break it continued to 1 kyr BP, while the heat flux increase began about 3 kyr earlier.The heat flux reached its maximum of 0.08 W m −2 about 8 kyr BP and then it began to decline.

The comparison of the SHF with solar insolation
The reconstructed SHF changes are similar to the Northern Hemisphere solar insolation changes that are determined by the variations of the Earth's orbital parameters like eccentricity, inclination, and the Earth's axis precession (Fig. 3).It is admissible to assume that insolation changes cause the surface heat flux changes.This assumption for the Middle Urals is also supported by the absence of Late Pleistocene ice sheets here (see Velichko et al., 1997;Svendsen et al., 2004, and references therein).However, there is some shift between insolation and SHF changes.The observed shift can be explained by several reasons.The first one is the influence of internal climatic factors and feedbacks translating the external heat flux on the Earth's surface with a certain delay and amplitude attenuation.The second reason is an overestimation of the effective thermal diffusivity that determines the rate of climatic signal propagation into the depth and therefore the timescale of geothermal reconstructions.The effective thermal diffusivity may be strongly affected by different factors like hydrogeological processes, thermophysical inhomogeneities of rocks, thawing-freezing processes, and so on.1)-( 5) SHF history q(t) (E = 2500 J m −2 K −1 s −1/2 , blue line).
To synchronize SHF and insolation ( I ) time series it is necessary to correct the initial value of thermal diffusivity (and timescale, respectively) to maximize the correlation between them.The correlation maximum points to an optimal value of thermal diffusivity, which determines the degree of expansion/compression of SHF and GST timescale.Note that the direct comparison of these series is not so correct.The insolation temporal resolution is constant while SHF resolution power decreases back in time.The minimum event length which can be resolved in the GSTH reconstruction has been estimated by Demezhko and Shchapov (2001) to be 2 × t * /3 where t * is time before present.The procedure of averaging in uneven running windows was proposed (Demezhko and Solomina, 2009) to modify the curve to a form comparable with the geothermal one.The insolation curve for the latitude of 60 • N smoothed according to the resolution power of geothermal method is presented in Fig. 3b.A maximum correlation (R = 0.99) between SHF history and smoothed insolation is achieved by increasing SHF dates by 1.4 times.It corresponds to the thermal diffusivity decrease from initial value of a = 1.0×10 −6 to 0.71×10 −6 m 2 s −1 .The interval of uncertainty for the optimal diffusivity value under R ≥ 0.95 is (0.71 ± 0.06) × 10 −6 m 2 s −1 .
Linear regression analysis of q and I from 35 to 6 kyr BP showed that change of insolation on 1 W m −2 produces an additional surface heat flux equal to 0.013 W m −2 .At that, only a small portion of insolation changes (about 1.3 %) was spent to the increase of the lithosphere heat content.The ratio q/ I may be considered as a dimensionless measure of climate sensitivity of the region under study to long-term orbital forcing variations.
Taking the climatically caused SHF before 35 kyr BP equal to 0 W m −2 and integrating it with respect to time, we estimate changes in heat content.This value characterizes the additional amount of heat adsorbed in a rock column having a cross-sectional area of 1 m 2 and limited by the depth of thermal anomaly penetration (i.e. by a few kilometres).Until 15 kyr BP a total heat balance was negative.A minimum value of heat content of −3.5 TJ m −2 with respect to the reference value at 35 kyr BP was found about 20 kyr BP.From this moment the heat flux became positive.For the next 14 kyr (20-6 kyr BP) the heat content increased to 22.0 TJ m −2 .For comparison, during the period of modern warming , heat content of the continental lithosphere increased by 0.1 TJ m −2 (calculated using data from Beltrami, 2002).

The comparison of the SHF with CO 2 changes
Another source of the additional radiative forcing during the Pleistocene-Holocene transition could be greenhouse effect caused by the increase of carbon dioxide concentration in the atmosphere (see Shakun et al., 2012, and references therein).An additional downward heat flux necessarily would contribute to SHF changes.Figure 4 shows geothermal reconstructions of surface temperatures and heat fluxes from the borehole SG-4 (on the timescale corrected after SHF-insolation synchronization) and carbon dioxide concen- tration changes in Antarctic ice cores (Blunier et al., 1998;Indermühle et al., 1999a, b;Smith, 1999;Barnola et al., 2003;Pedro et al., 2012).Despite the substantial dispersion of CO 2 estimations, a character and a chronology of CO 2 concentration changes are much closer to temperature changes rather than to heat flux variations.It may means no significant contribution of CO 2 forcing to climatically caused heat flux and thus to the temperature increase during Pleistocene-Holocene warming.
About 10 kyr BP the increase of carbon dioxide concentration was replaced by its fall which ended about 8 kyr BP.This local minimum is not consistent with either GST or SHF histories.It is possible that the CO 2 decrease was associated with a sharp increase of vegetation absorbing its excess.

Discussion and conclusions
The reconstruction of the surface heat flux history using data on the past surface temperature changes represents a new instrument for climate analysis.The reconstructed SHF variations and radiative forcing changes may be compared directly because they are expressed in the same units of energy flux (W m −2 ).
-Since the concentration of δ 18 O, δD in the ice cores or marine sediments associated with palaeotemperature fluctuations, the time shift between the orbital insolation and temperature reaction can be estimated only from independent absolute markers.Because of the rarity of such markers it is generally considered that the shift is a constant (Parrenin et al., 2007).Unlike conventional approach we tune another palaeoclimatic characteristics, the surface heat flux, which provides a physically reasonable shift.In Waelbroeck et al. (1995) the phasing between the precession band of mid-June insolation at 65 • N and δD was found about 3 kyr (with the uncertainty ±3 kyr).A reliable estimation of the phase in the obliquity band was not obtained and therefore it was not accounted for.Considering the period of precession 23 kyr and using Eq. ( 12) we obtain the close estimate 23/8 = 2.9 kyr.For the obliquity band the phase shift is equal to 41/8 = 5.1 kyr.
-For correct comparison with geothermal reconstruction the insolation curve must have the same resolution.A procedure of running averages with changing window lengths was applied to the insolation in order to obtain a curve comparable to the geothermal one.Such a procedure limits the tuning interval within the last cycle of precession.
-The reliability of the new timescale after synchronization with the orbital insolation also depends on how much the thermal diffusivity changed from the initial value 1 × 10 −6 m 2 s −1 , which is typical of common crustal rocks (e.g.Pasquale et al., 2014).The SG-4 borehole encountered 3.5 km of basalts, andesitebasalts, and their tuffs (Shakhtorina, 1992).Thermal diffusivity of these rocks may decrease to (0.76-0.70) × 10 −6 m 2 s −1 (Durham et al., 1987, Demezhko andSolomina, 2009), and even to 0.53 × 10 −6 m 2 s −1 (Eppelbaum et al., 2014).In our study the best coincidence of insolation and heat flux (R = 0.99) in the most part of the reconstructed interval is achieved under the thermal diffusivity value 0.71 × 10 −6 m 2 s −1 , which is quite realistic.
-Using the reconstructed surface heat flux instead of the surface temperature does not exclude the existence of residual time shift because the relation between insolation changes and the heat flux may be indirect.Such a shift can be caused by the climate delayed feedbacks.For example, orbital variations of insolation could change the extent of continental and sea ice cover in the Northern Hemisphere, albedo, and North Atlantic warm currents.The secondary heat source distributed in the atmosphere arose, which could significantly affect spatial distribution of the SHF change.However, this is beyond the scope of our study.
Assuming the surface heat flux varies proportionally to the external forcing one can consider the ratio q/ I as an alternative measure of the Earth's climatic sensitivity.The ratio of two heat fluxes is a non-dimensional parameter, and additionally depends less on radiative forcing duration by contrast to traditional index of climatic sensitivity represent-ing temperature reaction on the external radiative forcing ( T / I ).
The reconstructed surface heat flux reflects impact of all possible sources of radiative forcing.In addition to solar insolation, greenhouse gases (such as CO 2 ) can be a source of additional forcing.On the other hand the increase of carbon dioxide may be a consequence of temperature increasing.Comparing the chronology of surface flux, temperature and carbon dioxide concentration changes, we can draw some conclusions about the causes of climate changes.
The described algorithm of GST-SHF transformation is quite easy to realization and allows the estimation of SHF history with high precision.Using this algorithm, we have first estimated long-term surface heat flux changes in the Urals for the past 35 kyr.The amplitude of heat flux variations was about 1.3 percent of the insolation changes range at the latitude of 60 • N. The increase of carbon dioxide concentrations occurred 2-3 thousands of years later than the heat flux increase and synchronously with temperature response.

Figure 1 .
Figure 1.Testing the algorithms of GST-SHF transformation by applying it to a harmonic function of GST change.Relative error of SHF estimation (the ratio of the standard error of the SHF estimation to the real amplitude of SHF variations) versus the GST discretization frequency (points per period).

Figure 3 .
Figure 3.The comparison of SHF history with solar insolation changes in the Northern Hemisphere caused by changes in Earth's orbital parameters and timescale correcting.(a) Annual insolation changes I (t) at the latitudes of 40-70 • N (Berger and Loutre, 1991); (b) annual solar insolation at the latitude of 60 • N smoothed in uneven running windows (green line), SHF history in the initial timescale (a = 1.0 × 10 −6 m 2 s −1 , blue dashed line) and SHF history in the corrected timescale (a = 0.71 × 10 −6 m 2 s −1 , blue solid line).

Figure 4 .
Figure 4.The comparison of GST history T (t) (brown line), SHF history q(t) (blue line) and CO 2 concentration in the Antarctic ice cores (multicoloured markers).