Introduction
Measurement of the internal temperature of an ice mass usually requires some form of drilling to place temperature sensors in the ice. This creates a thermal disturbance in the ice, and the sensors register not the original temperature but a decay with time toward the undisturbed temperature. The time-scale of the decay can be approximately parameterized by normalizing the time since emplacement with the time span of the initial thermal disturbance (time is thus non-dimensionalized in this note). Reference ZotikovZotikov (1986) has used the concept of a non-dimensional “stabilization time” for drilled holes, in which the thermal disturbance has decayed to a negligible level after several tens or hundreds. For deep drilling projects that create a long disturbance, this stabilization time can be years in real time.
In practice, data are usually obtained for less than the full stablization time. Several factors, including drift in the sensors, logistics of returning to a distant site and straining of electrical cables from internal deformation in the ice mass, often lead to a short record. Hot-water drills, in particular, create large thermal disturbances extending well beyond the initial drilling, as a result of the refreezing of the borehole. In water-filled holes, the normalized record length can be short since the time span of the drilling disturbance may be many days, or the refreezing time may be long, such as when attempting to measure near-basal temperatures in ice that is close to the freezing temperature. A final practical consideration is that knowledge of the ice-temperature field is often necessary in drilling projects, but to get ice temperatures from recently drilled lióles requires estimation with a short record.
Temperature estimation by observation of the thermal decay of a disturbance has been extensively studied in permafrost by Reference Lachenbruch and BrewerLachenbruch and Brewer (1959). Their technique, or an even simpler model of the decay of an instantaneous line source of heat in an infinite medium (Reference Carslaw and JaegerCarslaw and Jaeger, 1959), has typically been applied in glacial studies. These decay models make simplifying assumptions about the thermal disturbance to obtain analytical solutions and, as a result, are limited to record lengths of greater than 5 or 10.
Refreezing of the borehole, with the attendant moving-boundary problem of the freezing walls, makes accurate modeling of shorter decay records a problem which can only be solved with numerical simulation. Reference Jarvis and ClarkeJarvis and Clarke (1974) outlined the basic numerical technique and used it to correct readings made with apparent, data lengths of 3 to 5; however, they did not elaborate the technique to shorter records or discuss some of the associated problems. In the extreme case where the hole does not freeze shut and the hole can be melted to maintain a constant diameter, the technique of Reference HarrisonHarrison (1972) can be applied which uses the refreezing rate of the borehole to get an estimate of the temperature.
However, there is usable information in temperature data from shortly after freez-up, and this note discusses a technique for estimating the original ice temperature using short data lengths of between 1 and 5 (the case where the real-time data record is short or the total disturbance time is long). The technique is illustrated with data from a borehole drilled through Ice Stream Β in Antarctica. An additional objective is to warn against the case with which short data records can be incorrectly fitted with simpler curves of the Reference Lachenbruch and BrewerLachenbruch and Hrewer (1959) type, where goodness of fit is no indication of the accuracy of projected temperatures.
Analytic Models of the Decay of a Thermal Disturbance
When a thermal record is shorter than the stabilization time, the data may be graphically compared with a thermal-decay model. The basic decay model is an analytical solution for the temperature in a solid that has been heated by an instantaneous line source of heat. The solution (Reference Carslaw and JaegerCarslaw and Jaeger, 1959) shows the temperature decrease is proportional to inverse time
where Δ Τ is the temperature difference between the hole axis and the distant, ice, Q is the total heat released per unit length, Κ is the thermal conductivity of the ice and t is time since the thermal disturbance. Typically, when the data length is of order 10 (normalized time) or more, temperature data plotted against inverse time asymptotically approach a straight line, which can be used to extrapolate to infinite time based on Equation (1). Knowledge of the heat input to the ice, Q, is unnecessary since the graphical technique solves for Q/4πK, as the slope of the line.
When the data length is less than order 10, it is necessary to account for the time span of the thermal disturbance. If the heat input to the ice is uniform in time, a better estimate of the short time behavior of the temperature decay is given in Reference Lachenbruch and BrewerLachenbruch and Brewer (1959)
where s is the time span of the (assumed constant) thermal disturbance, and the origin of time is at the start of the disturbance. Although originally developed for permafrost, Equation (2) with minor modifications has been applied to ice. Lachenbruch and Brewer used the drilling time for S; however, in ice the latent heat of the refreezing water in the hole releases a similar magnitude of heat as drilling, and this extends the effective disturbance time. In this note, the drilling time plus the refreezing time is used for s. However, the value of s for use in Equation (2) is not well defined when the thermal disturbance is non-uniform and this is one of several difficulties in using Equation (2) for short records. Equation (2) converges to Equation (1) as t ≫ s, and as Lachenbruch and Brewer pointed out, Equation (2) only differs from Equation (1) in terms of order (s/t)2 if t + s/2 is substituted for t in Equation (1). In practice, Equation (1) (with or without the s/2 term) has been used in most glaciological work. Both are referred to here as inverse-time models.
A normalized time of 5 is the approximate lower limit of applicability for the inverse-time models. At shorter times, temperature in the ice around a borehole deviates from the simple models in response to non-steady heating during drilling, phase changes in the borehole and non-zero hole diameter. This complex response is distinctive and can be compared to a numerical thermal model. Before discussing the numerical model, an example of the short time response of a temperature sensor in a borehole is introduced.
Temperature Decay in the Short Term, Measured by a Sensorl in a Freezing Borehole
Thermal drilling results in a temperature decay in the borehole with a distinctive shape. This is illustrated with data from a hole in cold ice, at Ice Stream B, Antarctica, that was drilled to the bed with a hot-water drill (Reference Engelhardt, Humphrey, Kamb and FahnestockEngelhardt and others, 1990). The dotted line in Figure 1 shows the data record, which is short (length about 2.1) because the site was occupied for less than 3 weeks after the completion of the hole. The sensor was installed immediately after drilling, 4.8 d after the drill first reached 240 m depth.
The dotted line shows data from a sensor at 240 m depth in a hot-water drilled hole at Up-B station, Ice Stream B, Antarctica. The data, from an AD590 temperature transducer and down-hole multiplexer, have had common mode electrical noise removed. Automatic recording was started an hour before sensor freeze-in. The dashed curve is from Equation (2), using the solution from the numerical model to obtain Q. The dotted curve is the numerical solution to the decay of temperature. The origin of time (t0) is when the drill first reached 240 m depth. Drilling lasted 4.8 d and freeze-up occurred at 6.73d (marked by the arrow).
The temperature record goes through two “kinks”. Initially, the sensor reads the freezing point of water, until the refreezing walls of the hole engulf the sensor (at t 0 + 6.3d). As is typical, the sensor is located somewhat off the hole axis, so that after the sensor is frozen the hole remains partially unfrozen for some time and the heat, liberated by the freezing water keeps the ice around the sensor relatively warm. Once the hole freezes completely (at t 0 + 6.73 d), the ice temperature drops rapidly as it decays towards the undisturbed temperature. The two diagnostic kinks, the first for sensor freeze-in and the second for hole freeze-up, are important features to identify. The second kink, which always precedes the steepest temperature drop, is used as the span of the thermal disturbance,s.
Numerically Extrapolating to Final Temperature from Very Short Records
A numerical model of the heat flow around a melting/freezing borehole was constructed (see Appendix) to investigate whether accurate estimations of ice temperature can be made from short time records. The model assumes that the ice/borehole interface is smooth and circular, heat flow is radially symmetric, and the water and ice are pure and have constant thermal properties. Details of borehole freezing, such as slush-ice formation, non-circularity in the borehole or spatial instability in the freezing/melting process, were neglected under the assumption that their importance to the temperature decay after borehole freeze-up would be small.
The inputs to the model are the original ice temperature and the heating history. A heating history is applied to a modeled region of ice and the melting and refreezing of the borehole are followed through time. Borehole size and the surrounding temperature field are recorded. After the end of heating and after the eventual freeze-up, the model follows the decay of the temperature field. The output is typically obtained in the form of a temperature versus time record at some chosen radial location, which can then be directly compared to recorded data.
Neither of the model inputs are known directly from the recorded temperature data. The original ice temperature is of course the unknown we are seeking, but the history is both complex and unknown, although estimates of the heat input can be made from the drilling records (see Reference Humphrey and EchelmeyerHumphrey and Echelmeyer, 1990). However, modeling shows that temperature decay following freeze-up is insensitive to the details of the heating of the hole during drilling. For example, comparisons were made between simulations in which heat was applied at a constant rate and simulations in which the hole was allowed to refrceze partially, with no direct heat input, and then re-opened with later heating. There was little difference in the post-freeze-up temperature response, as long as the total freeze-up time, s, was unchanged by the heating history. This behavior results from the fact that all heat lost to the ice has to traverse the constant-temperature water/ice interface. Variations in heat input to the borehole are compensated mostly by variations in melting/freezing at the walls, and not by changes in the heat flux into the ice. As a result, the heat lost by the borehole is largely determined by the hole size, which is itself parameterized by the time taken for the hole to rcfreeze after drilling stops.
Thus, for the purposes of modeling the subsequent temperature decay, the heating history can be effectively parameterized by the duration of heating (of unknown intensity) and the length of time from the end of heating to complete freeze-up. Since the freeze-up time can be directly read from the temperature record and the heating time is known from the drilling record, the only unknown is the original ice temperature, and the numerical model is essentially a one-parameter model. The model is sufficiently constrained that only a short temperature record containing the initial curvature of the post-freezing decay is required to obtain reasonable temperature estimates.
There are subtleties that limit the accuracy at very short record lengths. For example, the temperature decay does depend on the location of the sensor relative to the hole axis. The sensor location is included as an input parameter in the numerical model (location is chosen to match the first “kink” in the data). However, it has only a minor effect on the results and, after the hole has frozen, such details of the freezing process rapidly decline in importance.
The numerical model is used to produce a suite of temperature-decay curves based on varying the final ice temperature and a best fit to the data is chosen from the suite of outputs. In the model, the heat loss to the ice, Q, is given as an output. The model provides a good fit to data, as illustrated in Figure 1 and Figure 2, in which the final temperature is predicted to be —23.1° C with an error from the curve-matching of order 0.1° C (~0.5%), which is similar to the accuracy of the sensors.
The data from Figure 1 plotted as a fitnction of the time transform from Equation (2). The upper scale is in non-dimensional time (t/s).
For these calculations, the hole size was treated as an unknown. As a check, it should be noted that the borehole at Ice Stream Β was reamed to a nominal minimum size of 5.1 cm radius, while the model calculated a hole size of 5.7 cm, in good agreement with the nominal size. The kinks in the recorded temperature data are typically less abrupt than in the simulated temperatures; this is believed to be a result of unevenness in the freezing front and the finite size of the sensors. Any non-circularity in the borehole or ice-slush formation as the freezing front closes the hole would cause the recorded temperature to vary more smoothly than that derived from the mathematically precise freezing front in the model. In addition, the sensors were inside 2 cm diameter pressure housings which average the temperatures over a similar length scale.
Discussion and Comparison with Analytical Models
The dashed lines in Figure 1 and 2 show temperatures predicted by Equation (2), using Q found by the numerical modeling and s from the total refreezing time. Note that the dashed lines are not fitted to the data, as would be the case if data were being fitted to Equation (2), and Q (and s) was considered unknown. The measured and the numerically modeled temperatures are always colder than given by Equation (2) but the decay is slower. Thus, in the refreezing borehole, heat in the ice is stored out to a larger radius but at a lower temperature difference than in the model underlying Equation (2). This observational evidence shows that a thermal model of steady heating of a line source does not accurately describe the physical situation in the short time. Although the smooth curves of both the data and the numerical model are asymptotically approaching the dashed line of the inverse time model in Figure 2, it is probable that a line chosen to fit only the data would have a lower slope and thus a slightly warmer final temperature.
Choosing lines to fit the data is the same as choosing values for Q and s in Equation (2). Since the temperature-decay data are smooth and monotonie (after freeze-up), a very good fit to any short span of data can be obtained with the two-parameter inverse-time model, if Q and fl are considered “free”. The “goodness of fit” thus obtained is nonetheless spurious unless the temperatures are close to the asymptotic curve. Extrapolation will lead to an incorrect final temperature. It may appear possible to avoid this problem and extend the range of Equation (1) and Equation (2) by estimating Q independently of the temperature records, such as by using the drilling records. However, the estimation of Q is of similar complexity to the calculation of the temperature decay, and it is probably more accurate to perform the temperature calculation directly.
The potential for error from graphically extrapolating the data from Figure 2 using Equation (2) is small because the data record is long enough to approach closely the asymptotic line. However, when the data length is very short, the potential error can become large and the merit of the numerical technique is particularly apparent. Figure 3 illustrates a difficult case where the data are from a sensor in relatively warm ice at a depth of 870 m (160 m above the bed). The hole took 5d after the end of drilling to freeze-up, and a down-hole electrical failure occurred only 1 d after freeze-up. Despite the noise and the short record (data length of 1.13), the modeling shows that a reasonable extrapolation of the data is possible and provides a final temperature of —9.0° C with an error of about ±0.2° C (~20%).
A short and noisy data record from 870 m depth (160 m above the bed). The dotted line shows data from shortly before to shortly after freeze-up. The data before 6.6 d show a 0.6° C pressure /temperature depression. Data from the deep sensors were more contaminated with electrical noise than the sensor shown in Figure 1, probably as a result of water leakage into the signal lines. The solid curves show numerical solutions for ice temperatures of −8.6° and −9.0° C (upper and lower curves, respectively). The match between data and modeling, although not exact, implies an ice temperature of−9.0° C.
The data (plotted as in Figure 2) are shown in Figure 4, with a dashed line showing the asymptotic behavior as predicted by the inverse-time model (but with Q and s from the numerical modeling). The abscissa at the figure top is in non-dimensional time, or data length. Although cases vary in detail, Figure 4 is typical of the temperature response of a sensor in a refrozen hole. It is difficult to assess directly the error associated with using Equation (2), but Figure 4 shows that the data record ends while the temperature is still almost 1° C from the asymptotic line. In addition, a straight-line projection of the data leads to a predicted temperature up to 1° C (~20%) too warm. Indeed, from the numerical modeling, it appears that the difference in temperature between the asymptote and the numerical curve in Figure 4, at any particular data length, gives an estimate of the potential error in using an inverse time model. In this case of a very short data record, the necessity of using the numerical technique is obvious.
Data from Figure 3 plotted as a function of the time transform from Equation (2). The upper scale is non-dimensional time (t/s). In the short term, the data deviate significantly from the simple analytic models.
Conclusions
The length of record required to use accurately an analytical thermal-decay model ( Equation (1) or Equation (2)) for the extrapolation of temperature data from holes drilled with hot water, or from any thermally disturbed holes, may be longer than logistics or the durability of the sensors allow. The short time behavior of the temperature in the hole is not described by the simple inverse time models, but the short-term temperature decay can be fitted with simple models by using appropriate coefficients. Thus, there is potential for errors in extrapolating data to an undisturbed temperature by fitting to an incorrect model.
Any temperature-sensor record that includes both sensor freeze-in and hole freeze-up in a water-filled hole in cold ice contains enough information to yield estimates of the original ice temperature, based on comparison with a numerical simulation. In practice, it was found that data lengths only marginally greater than 1, say greater than 1.1, gives estimates with errors of only a few tenths of a degree or better. The technique requires little fitting, and can be routinely done in the field to check on progress of a temperature-measurement program.
Acknowledgements
This study was supported by the U.S. National Science Foundation under grant DPP-8519083.
The accuracy of references in the text and in this list is the responsibility of the author, to whom queries should be addressed.
MS received 14 March 1990 and in revised form 27 June 1990
APPENDIX A Numerical Model of the Thermal Decay Around an Initially Water-Filled Borehole
The thermal diffusivity of ice is well constrained outside of a small region of temperatures near freezing (Reference HarrisonHarrison, 1972), and mathematical models can be accurately applied. However, the melting and refreezing of a borehole is complicated by the existence of the moving boundary between the water and the ice walls of the hole. Although the problem is well posed, the “Stefan” problem of the moving boundary in a cylindrical coordinate system precludes any useful analytic description of the melting/freezing of a borehole during and after thermal drilling. Fortunately, the character of the problem does allow a particularly simple application of the finite element method.
Following a development by Reference Jarvis and ClarkeJarvis and Clarke (1974), and using the cylindrical symmetry of the problem, the one-dimensional governing equation for the heat flow in the ice is
where r and t are the space and time coordinates, Τ is temperature and a is the thermal diffusivity of ice. A description of the motion of the ice wall is given by a heat balance
where ρ is the density of ice or water, q is the heat flux to the ice wall from the ice or water and L is the latent heat of the phase transition. The heat-flow Equation (A1) can be transformed, using a log transform on the radial coordinate, into
where ω is the log of the radial distance. The transformed equation is particularly advantageous as a basis for dividing space in a numerical method. Equal-spaced nodes in logarithmic space create a high density of nodes near the freezing wall, where the largest temperature gradients occur. Equation (A2) and Equation (A3), plus boundary and initial conditions, were coded into a finite element formulation, which takes advantage of equal sized elements to avoid much of the overhead of matrix creation and solution. After each time step, the elements are remapped on to the solution space, so that the inner element always bounds the freezing front, while the outer element remains at the outer boundary.
The source of heat is the drilling hose during the melt phase (for more details, see Reference Humphrey and EchelmeyerHumphrey and Echelmeyer, 1990) and the latent heat of fusion during the freezing phase. Modeling is fast and stable until the hole freezes closed. At this point, the boundary condition (A2) is lost, which simplifies the problem except that has a singularity at the origin. After hole closure, the temperature decay may be solved analytically in terms of a Fourier-Bessel expansion of the temperature profile in the ice at freeze-up. However, it was found simplest to have the same code calculate the decay of the temperature after freeze-up. A special element of small radius (5 mm) was used to cover the neighborhood of the hole center to eliminate the mathematical singularity at r = 0.
Solutions were found using 25 elements occupying a radial space of 10 m surrounding the hole. The outermost node had a fixed temperature. The model pauses at freeze-up to integrate the heat in the ice, to obtain Q, and this sum is used to compare the numerical results from Equation (1) and Equation (2). The model is small and runs quickly on personal computers. Progress of a temperature measurement can therefore be modeled while in the field, and estimates of the temperature field can be found shortly after freeze-up.
