1. Introduction

During summertime, the interaction of Arctic glaciers with the atmosphere results in melting and subsequent freezing of meltwater, followed by a release of heat in the upper layers of snow–firn mass. As a consequence of intense summer melting, precipitation spreads inside the current annual layer and may also penetrate into other annual layers as a result of filtration of meltwater. This kind of ‘stirring’ of annual layers should be taken into account in paleo-reconstructions and in the analysis of ice cores. Methods for determining the past temperature of glacier surfaces, based on the data of measurements in boreholes, have been worked out mainly for cold glaciers with low summer temperatures. In this case, the processes of melting and penetration have almost no effect on the temperature field of the glacier and do not result in the stirring of annual layers (MacAyeal and others, 1991; Cuffey and others, 1994; Dahl-Jensen and others, 1999). To study temperature fields in Arctic glaciers, a model has been developed (Paterson and Clarke, 1978) that takes into account energy accumulation in the process of meltwater freezing. One of the input parameters of that model is the melt feature index (MFI), that is the relative concentration of refrozen ice per unit volume (Koerner and Paterson, 1974; Koerner, 1977; Watanabe and others, 2001). This parameter characterizes the intensity of melting and freezing, and determines the additional heat source in the subsurface layer of a glacier. Another input parameter of the model is the oxygen-isotope ratio (OIR) determined from the analysis of extracted ice cores. Linear dependence is known to exist between OIR and the near-surface temperature of air measured at the moment of formation of precipitation that constitutes annual layers of the glacier. The coefficients of linear dependence are influenced by many factors and can be considered as parameters of the site of location of the glacier. Both OIR and MFI are climatic proxies. If the temperature profile is known together with OIR and MFI parameters, it is possible to calibrate the paleothermometer for glaciers subjected to intense summer melting by applying a method analogous to that worked out for cold glaciers. In the present paper, a paleothermometer calibration method is proposed for two Arctic glaciers, the temperature fields of which are substantially influenced by the formation, penetration and freezing of meltwater. The boreholes in these glaciers were drilled at the summits of Akademii Nauk (Severnaya Zemlya, Russian High Arctic; 80.50^˚ N, 94.83^˚ E) and Austfonna (Svalbard; 79.85^˚ N, 24.14^˚ E) ice caps, respectively, by V. Zagorodnov in 1986/87 (Zagorodnov, 1988; Arkhipov, 1999).

2. Method of Calibrating Oxygen-Isotope Paleothermometer that Accounts for the Formation and Freezing of Meltwater

The well-known calibration of the oxygen-isotope paleothermometer for Greenland is based on an inversion of the heat-conduction equation (Cuffey and others, 1994). Due to low summer temperatures, this approach does not take into account the heat accumulation due to meltwater refreezing. We introduce this heat source using the Paterson and Clarke (1978) model. The other difference from the mentioned calibration consists in finding the solution in an explicit form.

Mathematical formulation of the problem of calculating the temperature field in an ice cap involves a one-dimensional heat-conduction equation that takes into account vertical advection of ice, and initial and boundary conditions. Let axis z be the origin at the surface of a glacier and be directed towards its base (z = H). The distribution of temperature in the bulk of the glacier is the solution of the following problem for the heat-conduction equation:

with the following boundary and initial conditions:

where T = T(z, t) is the glacier temperature, ρ and k are density and the thermal conductivity of firn, respectively, c is the heat capacity of ice, w is the rate of vertical advection, f (z, t) is the glacier heat-source density, ψ (t) = δ¹⁸O(t) is the oxygen-isotope ratio, a and b are constants, Q is the geothermal flux, T ₀ is the original surface temperature (considered as a parameter to be determined by the solution of the inverse problem), T _st(z) is the steady-state temperature corresponding to the geothermal flux Q and the surface temperature T _st(0) = 0°C, and t _f is terminal time which corresponds to the borehole temperature measurements. In this study, the vertical advection is assumed to obey the linear–quadratic model by Dansgaard and Johnsen (1969). Equations (1–4) form the so-called forward problem.

We approximate the firn density in ice caps by the following equation:

where ρ _ice is the density of bulk ice. On the basis of the density measurements in ice caps (Arkhipov, 1999), we choose the following values of parameters in Equation (5): γ ≈ 0.1 m^–1, c ₀ ≈ 0.58 for Austfonna ice cap and γ ≈ 0.28 m^–1, c ₀ ≈ 0.61 for Akademii Nauk ice cap. The thermal conductivity of firn k is determined by the relationship:

where k _f(ρ) = 0.021+ 0.00042ρ + (2.2 ×10^-9 ρ ³)Wm^–1K^–1 (Paterson and Clarke, 1978) and ρ is the density of firn. We introduced the temperature multiplier in Equation (6) according to Yen and others (1991). In order to linearize the heat-conduction equation, we use the following expression for the coefficient k:

where T _st(z) is the steady-state temperature. Hence, the problem (1–4) becomes linear with respect to temperature.

The function f (z, t) in the heat-conduction equation (1) determines the rate of temperature changes due both to deformation processes (Paterson and Clarke, 1978; Cuffey and others, 1994) and to refreezing of meltwater in the subsurface layer. Because the drilling sites of glaciers are located at the summits, the horizontal velocities are close to zero, and the internal heating of the glacier caused by ice shear has a negligible effect on the temperature distribution. According to experimental observations at the Severnaya Zemlya surface, the borehole site is located in a region where the horizontal velocities are zero (Dowdeswell and others, 2002). The surface velocity contours for Austfonna show that the horizontal velocity in the drilling location is zero (Dowdeswell and Drewry, 1989). Investigations of

Devon Ice Cap, Nunavut, Canada, show that during relatively long intervals of time (∼10³–10⁴ years) the refreezing of meltwater in the subsurface layer significantly affects the formation of the temperature profile in the ice cap (Paterson and Clarke, 1978). The specific heat (per unit of time and per unit of length) produced by refreezing of meltwater in the subsurface layer is (Paterson and Clarke, 1978):

where L is the latent heat, c is the heat capacity, w ₀ is the near-surface advection velocity equal to the rate of accumulation (w ₀ ≈ a ₀), ρ ₀ is the density of firn in the subsurface layer, 〈P(t)〉 is the mean value of the relative concentration of refrozen ice by volume unit, the MFI, and g(z) is the distribution function of meltwater in the subsurface layer. The function g(z) related to the area of localization of meltwater in the subsurface layer is not unambiguously defined. For example, in order to take into account the amount of heat of refreezing in Devon Ice Cap, Paterson and Clarke (1978) used for g(z) the ‘triangular’ form

where d is the depth of localization of meltwater in the subsurface layer, and l ₀ = 0.2m is the effective length of the region where the water–ice transition takes place. The value ρ ₀ in Equation (˚) is equal to the density of firn ρ(z) satisfying the relationship z є [z ₁, z ₂], where [z ₁, z ₂] is the interval of non-zero values of the function g(z), and ρ ₀ ≈ ρ(d) for the relatively narrow triangular area.

The measurements of relative concentration of recrystallized ice carried out on 2 m long cores (Arkhipov, 1999) show that the function P(t) significantly fluctuates on time intervals of 1–10 years. Taking this into account, in Equation (˚) we use the mean value of the function where τ is the time-step of the numerical scheme. For Akademii Nauk and Austfonna ice caps, we use for 〈P(t)〉 the ‘smoothed’ values of the function P(t):

where j is an ice-core section number corresponding to time interval [tj , t _j-1].

One can simplify the source function f (z,t) in Equation (˚) by considering possible values of temperature on Akademii Nauk and Austfonna glaciers and the fact that Accordingly, as a zero approximation the amount of heat resulting from refreezing may be taken as f (z, t) = La ₀ ρ ₀〈P(t)〉g(z). An additional parameter P ₀ should be introduced into this equation. This parameter is unknown and must be determined by the solution of the inverse problem. The necessity to introduce P ₀ stems, on the one hand, from possible deviation between the estimated and real values of the amount of meltwater measured by the content of infiltrated ice in the ice cores (Paterson and Clarke, 1978). On the other hand, as already mentioned, the parameters l ₀ and d of the refreezing zone in Equation (9) and the corresponding relationship ρ ₀ ≈ ρ(d) are not unambiguously determined. The parameter P ₀ has to be around 1, and its variation allows estimation of the possible range of the depth of localization of meltwater. Thus, the expression for function f (z, t) takes the form:

The system of Equations (1–4) is linear for the unknown temperature. Therefore, the linear superposition is applied. Let us look for the solution of Equations (1–4) in the form

where T*(z, t) is the solution of Equations (1–4) with a = b = T ₀ = P ₀ = 0; T ₁(z, t) the solution of Equations (1–4) with φ(z) = φ(t) = Q = T ₀ = P ₀ = 0, a = 1; T ₂(z, t) the solution of Equations (1–4) with φ(z) = a = Q = T ₀ = P ₀ = 0, b = 1; T ₃(z, t) the solution of Equations (1–4) with φ(z) = a = b = Q = P ₀ = 0, T ₀ = 1; and T ₄(z, t) the solution of Equations (1–4) with φ(z) = a = b = Q = T ₀ = 0, P ₀ = 1. All conditions are known to find the functions T(z, t), T ₁(z, t), T ₂(z, t), T ₃(z, t), T ₄(z, t). That is done numerically by finite-difference schemes.

In the case of the forward problem (1–4), a set of parameters P ₀, a, b and T ₀ is known, and temperature can be calculated at any given space point z and time t. In our case this set is unknown. Its determination is a solution of the inverse problem that includes Equations (1–4) and an additional condition

where T _bh(z) is the measured temperature–depth profile. Let us introduce the discrepancy functional

Note that S 0. It becomes zero only if the set a, b, T ₀ and P ₀ is the solution of the inverse problem (1–4), (13). In all other cases S ≤ 0. Thus, the required parameter set provides a minimum of the functional (13), and then

This leads to the system of equations

where Thus, we obtain explicit formulas for a, b, T ₀, P ₀, as well as for the surface temperature T(0, t).

3. Results of the Reconstruction of Past Temperatures for Akademii Nauk and Austfonna Ice Caps

The measured OIR δ¹⁸O and the MFI in the borehole of Akademii Nauk ice cap are shown in Figure 1. Both proxies exhibit elevated values in the upper part of the glacier that can be interpreted as warming in recent years. Similar behavior is observed in the measured temperature–depth profile (triangles in Fig. 2): the upper-part temperature is higher than the temperature corresponding to the steady state (dashed curve). The steady-state temperature is calculated here by the measured bottom temperature and the geothermal flux calculated by the measured temperature gradient. One can see the difference between the measured and the steady-state fitted temperatures for a specified part of the temperature profile. The maximal absolute value of the difference does not exceed 0.2˚C, which agrees well with the accuracy of the measurements.

Fig. 1.

The climatic proxies for Akademii Nauk ice cap.

Fig. 2.

Temperatures in Akademii Nauk ice cap (1. measured; 2. fitted steady-state; 3. misfit).

The results of calculations of the coefficients of paleothermometer and temperature profiles are depicted in Figures 3 and 4. The rate of accumulation on the surface of Akademii Nauk ice cap lies in the range 0.3 ≤ a ₀ ≤ 0.5ma^–1 (Arkhipov, 1999). Assuming the value a ₀ = 0.3ma^–1, the relative error in coefficient a of the paleothermometer, related to the uncertainty in the data for δ¹⁸O, does not exceed ˚%. Given a ₀ = 0.3ma^–1, the mean-square deviation S between the calculated and the measured temperature profiles is minimal. The corresponding values of thecoefficient α of inverse dependence of δ¹⁸O on temperature lie in the range α ≈0.57–0.63 ppt ˚C^–1. This agrees with the results of calibration of the paleothermometer for Greenland (Cuffey and others, 1994).

Fig. 3.

Measured and calculated temperature profiles in Akademii Nauk ice cap for different accumulation rate (1. a ₀ = 0.3ma^–1; 2. a ₀ = 0.5ma^–1; 3. measured).

Fig. 4.

Reconstructed temperatures, Akademii Nauk ice cap (1 and 2. surface and 10m temperatures based on the isotope calibration; 3. the 10 m temperatures reconstructed by the Tikhonov regularization method (Nagornov and others, 2001)).

The variations of the depth d of the area of localization of meltwater (i.e. the area of the heat source related to refreezing of meltwater) in the range 0.1–0.8m result in the corresponding variations of parameter P ₀ in the range 0.1–1.43. This variation of parameter d assumes that the square of the triangle g(z) remains constant. In the range of variations of d between 0.1 and 0.8m the coefficients of the paleothermometer are practically constant (vary no more than 0.2%). The minimum of the functional S is achieved for d = 0.3 m. Taking into consideration the annual accumulation rate a ₀ = 0.3ma^–1, we conclude that the meltwater is localized predominantly in the current annual layer, and only its minor part penetrates into the preceding annual layer. The reconstructed surface temperature (Fig. 4, curve 1) agrees well with the temperature reconstructed previously using another method (Nagornov and others, 2001). In that study, the temperature was reconstructed at the depth of attenuation of seasonal signals (at the ‘10 m’ depth (Blatter, 1987)) and then recalculated for the free surface (Fig. 4, curve 2). Thus we conclude that despite the presence of melting processes and meltwater filtration, the approach of the oxygen-isotope paleothermometer is also justified for Arctic glaciers subject to intense summer melting. The parameters a and b of the oxygen-isotope paleothermometer (T = aδ¹⁸O+b) are found for Akademii Nauk ice cap (a = 1.6–1.73˚Cppt^–1, b = 19.8–23.6˚C). The variations of a and b are related to temperature measurement errors in the boreholes and to the uncertainty in the values of accumulation rate and geothermal flux.

The reconstruction of the surface temperature of Austfonna ice cap (Arkhipov and others, 1987; Arkhipov, 1999) is performed in a similar manner. The results of calculations of the parameters of the paleothermometer are obtained for accumulation rate values lying in the range 0.2–0.8ma^–1 (Arkhipov, 1999). The minimum of the deviation between the calculated and the measured temperature profiles is achieved at a ₀ ≈ 0.5–0.6ma^–1. Assuming the value of accumulation rate a ₀ = 0.5ma^–1, the relative error in the coefficient a of the paleothermometer, related to the uncertainty in the data for δ¹⁸O, does not exceed 12%. The corresponding values of the coefficient α of inverse dependence of δ¹⁸O on temperature lie in the range α ≈0.3–0.4 ppt ˚C^–1.

The measured values of temperature at the bottom of the glacier allow for assuming the melting of ice at its lower boundary (Arkhipov, 1999). Let us investigate the influence of the variations of the geothermal flux on the coefficients of the paleothermometer in the case a ₀ = 0.5ma^–1. For values of the geothermal flux in the range 10–50mWm^–2, the coefficient a varies from 2.6˚C^–1 to 2.9˚C ^–1, while the coefficient b varies from 30˚C to 33˚C. The values of the coefficient α of inverse dependence of δ¹⁸O on temperature fall in the range 0.34–0.4 ppt ˚C^–1. Therefore, if one assumes the value Q ≈ 30mWm^–2 for the geothermal heat flux at z = H, the relative error in determination of the coefficient a (related to the uncertainty of the value of Q at z = H) amounts to 6–7%. The maximal error in determination of the coefficients of the paleothermometer resulting from errors in the measurements of the temperature profile in the lower part of the glacier (z ≥ 200 m) does not exceed 7%.

Variations of the depth d of the area of localization of heat source associated with refreezing in the range 0.5–2m lead to corresponding changes of the parameter P ₀ in the range 0.54–3.3. In this case, the coefficients a and b of the paleothermometer vary in the ranges 2.6–2.77˚C ^–1 and 28.5–32.7˚C, respectively, and deviation between the calculated and the measured temperatures does not exceed 0.01˚C. Based on the obtained values of P ₀ and considering the rate of accumulation a ₀ = 0.5ma^–1, we deduce that the depth of the heat source associated with refreezing may lie in the interval 1 < d < 2m. Therefore, meltwater from the surface of Austfonna ice cap is likely to penetrate into three to four annual layers. The variations of the parameter l ₀ have almost no influence on the results of the paleothermometer calibration. Indeed, given a ₀ = 0.5ma^–1 and d = 1.5m, the coefficients of the paleothermometer corresponding to l ₀ = 0.2m and l ₀ = 3 m, respectively, differ by about 1%. For Austfonna ice cap, the best agreement between the calculated and the measured temperature profiles is achieved when the values of the coefficients of the paleothermometer lie in the ranges a = 2.6–2.9˚Cppt^–1 and b = 30–33˚C.

Temperature variations at the depth of attenuation of seasonal signals agree with temperature changes at the surface of Austfonna ice cap reconstructed using the measured δ¹⁸O data (Fig. 5, curves 1 and 2). The temperature is higher at 10 m depth than at the surface of the glacier. This effect is associated with the existence of an additional heat source due to refreezing of meltwater and advective heat transfer towards the bottom of the ice cap. Higher mean annual temperatures at the surface of Austfonna ice cap result in stronger influence of the heat source due to refreezing of meltwater than at Akademii Nauk ice cap. The additional heat source due to the ‘meltwater–ice’ phase transition at the surface of the glaciers under consideration leads to significantly more pronounced warming than in the Arctic region and Northern Hemisphere over the last 150 years (Overpeck and others, 1997).

Fig. 5.

Same as Figure 4, for Austfonna ice cap.

4. Conclusions

A method has been developed for determining past temperatures at the surface of Arctic glaciers subject to intense summer melting. The method is based on the measured oxygen-isotope ratio, the MFI in the ice core, and the borehole temperature profile. It was found that meltwater spreads inside two annual layers for Akademii Nauk ice cap and inside three to four annual layers for Austfonna ice cap. The values of the past temperature at the depth of attenuation of seasonal signals determined in the present paper agree well with those found by using another method (Nagornov and others, 2001). This justifies the extension of the approximation of the oxygen-isotope paleothermometer, previously used only for cold glaciers (Cuffey and others, 1994), to glaciers with intense summer melting. The following parameters of the isotope paleothermometer have been found: a = 2.6–2.9˚C ppt^–1, b = 30–33˚C for Austfonna ice cap; and a = 1.6–1.73˚C ppt^–1, b = 19.8–23.6˚C for Akademii Nauk ice cap. Refreezing of meltwater in subsurface layers of Arctic glaciers results in additional temperature rise (about 6˚C) as compared to the average warming in the Arctic region and Northern Hemisphere over the last 150 years.