Measured temperatures

We use the temperature-depth profile measured by Reference Alley and KociAlley and Koci (1990) in a 217 m deep borehole at the GISP2 site near the summit of the Greenland ice sheet. Alley and Koci found the shape of the measured profile to be reproducible to ±0.01°C, although the absolute value of the temperatures was consistently ofiset by 0.035°C for two different factory-calibrated thermistors. Temperature gradients in the borehole are ≤0.04°Cm^∡1, less than the adiabatic lapse-rate of air at these temperatures. In addition, the borehole is narrow, limiting the potential size of convection cells. Convective disturbance of the borehole-temperature profile is therefore unlikely (Reference DimentDiment, 1967). A repeat measurement of the borehole temperatures 1 year later by Alley and others confirmed that the trend was undisturbed. For both years, we observed that the thermistors’ resistance values were stable during measurement. Finally, the error resulting from distortion of the temperature distribution in the ice due to the presence of the air-filled borehole is insignificant for these low gradients (Reference SandersonSanderson, 1977). For all these reasons, we consider the borehole-temperature profile to reflect accurately the temperatures in the surrounding ice.

The model

Our thermal model is a one-dimensional finite-difference equation, similar to that used by Reference JohnsenJohnsen (1977), Reference Paterson and ClarkePaterson and Clarke (1978), Reference BolzanBolzan (1985) and Reference Alley and KociAlley and Koci (1990). Thermal energy per unit volume changes with time as (Reference MalvernMalvern, 1969, p. 228–30)

where ρ is ice/firn density, c is heat capacity, Τ is temperature, t is time, t _i and t ₀ correspond respectively to the years AD 649 and AD 1989, z is depth (0 at the ice-sheet surface, increasing downward), k is thermal conductivity, w is vertical velocity, σ is the stress tensor, is the strain-rate tensor, and: indicates the inner product. is the rate of energy production due to non-mechanical sources. The first three terms on the righthand side of Equation (4) account for heat transfer by conduction and advection and heat generation due to mechanical work during strain. Values for heat capacity are taken from Reference YenYen (1981, p. 13, table 2). For thermal conductivity, we use empirical data for firn densities ρ 750 kgm^−s (Reference YenYen, 1981, p. 16, equation 35). For ρ750 kg m⁻³, we use empirical k values for pure ice (Reference YenYen, 1981, p. 15, equation 33), and a theoretical extrapolation for dense firn (Reference YenYen, 1981, p. 17, equation 37).

The initial temperature, T_i(z), is a steady-state profile for an accumulation rate b = 0.24 m a⁻¹ and surface temperature T_i(0). Equation (6) expresses T_i(0) as an offset, Δ, from the average temperature during the model run, T (a, b), which is a function of the coefficients in Equation (1). Note that Δ 0 indicates the average temperature over the last 1340 year was less than that during the preceding years. We have shown that Δ probably falls in the range 0.65–0.90°C, with a most probable value of 0.79°C (Reference Cuffey, Alley, Grootes and AnandakrishnanCuffey and others, 1992). This is similar to the result of Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen (1986), suggesting a Medieval Warm Period of this magnitude at Dye 3. Because the shape of the steady-state profile is fixed, Teoo is uniquely determined by T_i(0) (Equation (7)).

The remaining terms in Equation (4) are given by:

Here, σ_z is vertical stress and σ_x is horizontal stress normal to the ice divide (i.e. east—west) and is the horizontal longitudinal strain rate in this direction. We partition the vertical longitudinal strain rate into components due to firn densification, and flow divergence, . We assume plane strain in vertical planes normal to the ice divide, so no terms appear for the north-south horizontal direction. Because our model includes only the top 20% of the ice thickness and because GISP2 is only 30 km west of the ice divide, we neglect the work terms involving shear deformation, which are orders of magnitude less at these depths and low surface slopes. Note that the product in Equation (9) gives the rate of heat dissipation due to work done moving grains closer together in the densification process, while the sum of the other two terms in Equation (9) gives the rate of heat dissipation due to the flow-divergence strain.

Incomprcssibility requires that the strain rates due to flow divergence sum to zero (Equation (10)). We assume that has a constant value equal to the mean annual ice-equivalent accumulation rate, b, divided by the ice-equivalent thickness, h_i = 3100m (Reference Hodge, Wright, Bradley, Jacobel, Skou and VaughnHodge and others, 1990). We also assume power-law flow (n = 3) to write in terms of the vertical longitudinal deviatoric stress, σ_d , and the usual flow-law pre-factor, A, taken from Reference PatersonPaterson (1981, p. 39) for a temperature of-30°C.

In Equation (11), wc write the vertical strain rate due to firn densification in terms of the material derivative of density (Dρ/Dt) and assume that the density at a given depth docs not change through time. Equation (12) expresses the balance of forces; the vertical normal stress equals the weight of overburden and the sum of the longitudinal deviatoric stresses vanishes. The vertical velocity equals the snowfall rate at the surface and decreases with depth as the integral of the vertical strain rate (Equation (13)).

Finally, we assume that thermal energy released by non-mechanical sources is negligible (Equation (14)). At locations with frequent surface melt, the latent heat released upon refreezing of meltwater at depth can dominate the temperature structure (Reference Paterson and ClarkePaterson and Clarke, 1978). The paucity of inferred melt layers in the GISP2 ice cores (Reference Alley, Grootes, Meese, Gow, Taylor and CuffeyAlley and others, 1991) indicates we can neglect this term. Another source of thermal energy is the elimination of surface area during ice-grain growth We neglect this term, too, based on the following argument. The surface area of grain boundaries per volume of ice, A_grains , is approximately half the area of each grain (because each surface is largely shared by neighbouring grains) times the number of grains per volume. For spherical grains of radius r, A _grains ≈ 1.5/r.The rate of energy release is the product of surface energy associated with the ice—ice boundaries (γ_i = 0.065 J m^-2;, Reference HobbsHobbs, 1974, p. 440) and rate of elimination of grain-boundary area, or, for a steady distribution of grain-size with depth,

will be a maximum for shallow ice, where r ≈ 10^-3 m, and r increases by approximately 1 mm per hundred meters, giving a temperature change at a rate (∂T/∂t) , which is more than an order of magnitude less than strain heating (see below and Fig. 1).

Fig. 1.

Contribution of different heat-transfer and generation mechanisms to rate of temperature change (Ka⁻¹) in the ice sheet, calculated for AD 1989. Absolute values for temperature-change rate were taken to make possible the logarithmic scale. This causes the sharp downward points where the sign of temperature change switches. We omit the lower 200 m to facilitate the view of the more interesting upper section.

Our model extends from the ice-sheet surface to 600 m depth, with uniform 6 m node spacing. Time step is 1 x 10⁶s (11.6 d). The density profile is measured and assumed constant, which is a reasonable assumption because long-term averages of temperature and accumulation rate do not vary significantly in the Holocene. Conductive-and advective-heat transport dominate the model, with heat generation due to densification playing a significant role in the 20–150 m depth range (see Fig. 1).

Forcings

Through time, the model is forced with surface temperatures, T_s , calculated from isotopic ratios (Equation (8)) and accumulation rates, b(t) calculated from seasonal indicators in the ice core (Reference Meese, Gow, Alley, Taylor and RamMeese and others, 1992) (see Fig. 2a and b).The age scale is accurate to a year or two for the most recent times and to a few per cent for older times (Reference Meese, Gow, Alley, Taylor and RamMeese and others, 1992). The isotopic ratios are measured using mass spectrometry. The uncertainty in a single δ ¹⁸O value is less than 0.1‰. We measured ten samples per year of ice and interpolated between these to get our temperature forcing. Mean annual δ ¹⁸O values have an analytic uncertainty of less than 0.03‰. Note that, if our δ measurements are systematically offset from true values, the intercept (b) value in Equation (1) will be offset by the same amount in our calibration.

Fig. 2.

The isotope and accumulation-rate data used in our calibration, smoothed using 10 year running average.

Inversion

Comparison of model equation solutions with measured borehole-temperature data enables us to invert for optimal values of the coefficients a and b Δ is also treated as a free parameter in the inversion. We use a multi-dimensional form of Newton’s method (Reference MenkeMenke, 1984) to find values of model parameters that give best agreement between model and measured temperatures, as defined by Equation (3). The inversion uses model grid points only from 36 to 216m depth, the bottom of the measured temperature record. The top 30 m are truncated to avoid extreme sensitivity of the calibration to the most recent several years. This is desirable because the natural variability of mean annual δ ¹⁸O values is large; because the shallowest part of our borehole may be affected by surface winds; and because the amplitude of shallow temperature variations is huge compared to those at greater depths.

An initial guess at values for the model parameters, m_q ⁰ = [a, b, Δ], is changed by an amount , calculated from

where p = 1,... Ν and q = 1,... Μ (Ν = number of grid points = 31, M = number of free-model parameters≤3). Here is the temperature residual, T_p (m_q ⁰) - ϴ_p , and H_pq is a matrix of the derivatives of T_p with respect to each model parameter, calculated numerically as

which we invert using singular-value decomposition to solve Equation (16). The new values are then modified in the same fashion with decreasing until E _ms ceases to decrease. The inversion to solve Equation (16) is overdetermined. We neglect no eigenvalues and we find that the model resolution matrix is approximately equal to the identity matrix, with off-diagonal elements smaller than 10⁻⁹.

We interpolate measured temperatures to each grid point, because the measured values have 5 m or less spacing, versus 6 m spacing for grid points. Because temperature varies gradually, this interpolation introduces insignificant error. Inversions of synthetic borehole temperatures produced by known values of a, b and Δ indicate that the numerical error of our technique is approximately ±0.001‰°C⁻¹ for b, ±0.05‰ for b and ±0.005°C for Δ. In addition, the solutions do not change with variations in initial m_q0 or with reasonable changes in the size of . We therefore consider our solutions to be accurate numerically.

In this paper, we map the error E_ms as a surface over the a-b plane. The low point of this surface, which has error E_ms ^min overlies the optimal values of the coefficients a a and b, a* and b*. This minimum is found by using m_q = [a, b, Δ] in the inversion. Then we assign fixed values to a and use only m_q = [6, Δ]. This defines a trough in the error surface that gives the best possible error for given a. It is important to define this trough because a is the parameter in our inversion of most interest to paleoclimate studies, since it relates the magnitude of variations in δ ¹⁸O to variations in temperature. Finally, we assign fixed values to both a and b, let m_q = Δ, find the minimum E_ms for these (a, b) coordinates and contour the resulting surface. This surface shows how confident we can be in our calibration. If the optimal values, a* and b*, are well-constrained, the error surface will rise rapidly along trajectories away from the coordinate (a*, b*). For us to claim that (a*, b*) makes a significantly better paleothermometer than (a, b) at a given level of statistical confidence, the ratio E_ms(a, b)/E_ms ^min should exceed the maximum ratio of squared errors expected from random fluctuations in the data alone, which is given by the F-distribution (see Reference MenkeMenke, 1984, p. 96–97). Thus we use an F-test that compares the surface E_ms(a, b)/E_ms ^min to outline regions in the a-b plane that contain the values of the true paleothermometer at various confidence levels. For optimal values of Δ, see Reference Cuffey, Alley, Grootes and AnandakrishnanCuffey and others (1992).

Time-weighting of the calibration

If the slope and intercept of the isotope-temperature relation (Equation (1)) have changed over time, then a and b from our calibration are time-averaged values for the most recent 1340 years. Furthermore, these are weighted averages, with some times contributing more than others. For example, very recent times are weighted lightly in the calibration because we disregard the uppermost part of the borehole-temperature profile where the effects of recent temperature change are largest. We calculate the weight function, W(τ), as a semi-quantitative measure of the influence of events in each of the 1340 years on the modern borehole-temperature profile; W(τ) therefore approximates the time-weighting of a and b in our calibration.

The surface-temperature history given by Equation (1) from the 1340 year δ ¹⁸O record is simplified to a series of annual average perturbations about the mean, δΤ_av(τ), for the rth year before present. In the modern ice sheet, each δΤ_av(τ), will have a temperature response at depth z

Here, is the temperature perturbation as a function of elapsed time that results from a unit forcing (1°C for 1 year) at time . These are Green’s functions, which we compute numerically. Note that, because τ is given as years before present, τ equals .

The modern temperature at depth z depends on the sum over all years of the responses R_z(τ). Comparing the magnitude of R_z(τ) for different years τ therefore gives a measure of how important each year is for determining the modern temperature structure, compared to other years, at that depth. To present this comparison graphically, incorporating all depths used in the calibration, we define the depth-averaged weight function W(T), as

, where z(p) is the depth of the pth grid point used in the calibration.

We normalize W(τ) by setting the integral of W(τ) over the entire range 0 < τ < 1340 years equal to unity. Then gives the relative contribution of the interval τ₁ < τ < τ₂ to the modern temperature structure and hence to the calibration of the isotopic paleothermometer.

Modern local spatial gradient of δ ¹⁸O with temperature

One purpose of our paper is to argue that calibrating isotopic paleothermometers using borehole temperatures is a profitable undertaking. To this end, it is useful to compare results of our borehole calibration with results of the most frequently used calibration method, that based on spatial gradients (as defined in the Introduction).

We calculate a spatial gradient using data from locations within 100 km of the GISP2 site. One of us (unpublished manuscript by J. Bolzan) has measured 10 m temperatures in the GISP2 area and collected a number of shallow cores that span approximately 2545 years, and oxygen-isotope profiles were measured on these cores by the Geophysical Laboratory in Copenhagen (Reference Bolzan and StrobelBolzan and Strobel, 1994). Unfortunately, for logistical reasons, the isotope and temperature measurements were not made at the same sites. As a result, we interpolated the temperature field using a generalized Kriging algorithm (Reference OleaOlea, 1974) to generate mean annual temperatures at the shallow coring sites, except for two sites which we considered too distant from the nearest measured 10 m temperature. We then regressed isotopic values averaged over the most recent 30 year interval in the shallow cores against the interpolated mean annual temperatures. The one core that spanned only 25 years was not used.

Reference Clausen, Gundestrup, Johnsen, Bindschadler and ZwallyClausen and others (1988) measured 10 m borehole temperatures and δ ¹⁸O of ice cores south and east of the GISP2 site. We include these data in our correlation to widen the temperature range and increase the sample size. The isotopic data of Reference Clausen, Gundestrup, Johnsen, Bindschadler and ZwallyClausen and others (1988, fig. 4) are 30 year averages; we use the most recent value for each core. We do not use the data for Crête, because this core was drilled almost a decade before the others and therefore does not include the late 1970s and early 1980s in the 30 year average.

Methodology for the sensitivity tests

The values we use for the borehole temperatures and for some fixed parameters in our model may differ from the true values. It is important to know how this uncertainty may affect our calibration. Therefore, we test the sensitivity of the optimal values, a* and b*, to these uncertain model variables. We alter each model variationsable in question, V, by an amount or replace it by an alternative variable, V_alt , and repeat the inversion (m_q ⁰ = [a,b,Δ]) using or V_alt . If the result has E _ms that is statistically significantly worse than E _ms ^min, then the altered value is probably physically unreasonable. Otherwise, the inversion gives new, acceptable values for a* and b*.

We examine the sensitivity of the calibration to:

• 1. The basal boundary condition. We replace the constant basal temperature T₆₀₀ , with a constant basal heat (lux implied by the steady-state initial temperature profile.

• 2. The initial temperature profile. In Equation (5), we replace the steady-state, nearly isothermal, profile with linear profiles where T _i is the temperature difference between the 600 m depth and the surface. It is unlikely that because there was a cooling trend during the late Holocene (Reference LambLamb, 1982; Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen, 1986). It is also probable that is less than 1°C because the difference in mean annual temperature between the climatic optimum and the first millennium AD is at most several °C (Reference LambLamb, 1982; Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen, 1986). The magnitude of the corresponding temperature maximum at depth will be considerably reduced by diffusion and the maximum will be deeper than 600 m due to advection alone. Therefore, we use and-1.0°C as estimates of this uncertainty.

• 3. The thermal properties of snow. For snow of density ρ < 750 kg m^-3, we increase and decrease the dependence of thermal conductivity on density and temperature by the maximum amount that still approximates the data trends in figures 16 and 17 of Reference YenYen (1981). We alter the constants in Yen’s equation (35), so that thermal conductivity is 0.048 exp(0.0103T + 5.9289 ρ) for the increased dependence case; for the decreased dependence case it is 0.100 exp(0.0073T + 4.1456ρ) (units as in Reference YenYen (1981)).

• 4. The ice-equivalent thickness, h_i The thickness of ice beneath Summit is uncertain by 100 m or more because of the uncertainties in velocity of radar waves through ice and the moderately coarse grid for radar surveys relative to the rough bedrock topography (Reference Hodge, Wright, Bradley, Jacobel, Skou and VaughnHodge and others, 1990). We have used h_i = 3150m (Reference Hodge, Wright, Bradley, Jacobel, Skou and VaughnHodge and others, 1990). Here, we use with .

• 5. The borehole-temperature measurements.

Results of the sensitivity tests

Results of the sensitivity tests show that realistic uncertainties in measured temperatures and model variables do not significantly affect our calibration (see Table 1). a* is insensitive to the basal boundary condition, the ice thickness and the absolute value of the measured borehole temperatures. It is weakly dependent on the initial temperature profile and on snow thermal properties, and more strongly dependent on systematic errors in the measured borehole-temperature trend. Because the temperature trend was highly reproducible using different measuring devices, we consider such systematic errors to be very unlikely. Assuming no such systematic error, the summed uncertainty in a* ranges from +0.040 (which obtains if ) to-0.038 (which obtains if , snow-thermal properties have minimum dependence on density and temperature, ). Thus, the sensitivity of a* is considerably smaller than the 90% confidence interval for a* (Fig. 4).

Table 1.

Tabulated results of sensitivity tests. Changing the model variables listed returned the new optimal values a* + da and b* + db, where a* = 0.5305, b* =-18.18. Recall that the numerical error in a* is of order ±0.001 and in b* is of order ± 0.05. Cases where da or db exceed these values are shown in bold type. Note that increasing the dependence of thermal conductivity on snow density and temperature leads to model results that are statistically unacceptable at the 95% confidence level

b* is insensitive to the basal boundary condition and ice thickness, moderately sensitive to the absolute value of borehole temperatures and more strongly sensitive to deviations of the initial temperature profile from steady-state and to snow-thermal properties. It is most sensitive to systematic errors in measured temperatures. The summed uncertainty in b* ranges from + 1.84 to-0.72, assuming no systematic errors, which is much larger than the 90% confidence interval for b* (Fig. 4).