Numerical study of the error sources in the experimental estimation of thermal diffusivity: an application to debris-covered glaciers

A surface debris layer significantly modifies underlying ice melt dependent on the thermal resistance of the debris cover, with thermal resistance being a function of debris thickness and effective thermal conductivity. Thus, these terms are required in models of sub-debris ice melt. The most commonly used method to calculate effective thermal conductivity of supraglacial debris layers applies heat diffusion principles to a vertical array of temperature measurements through the supraglacial debris cover combined with an estimate of volumetric heat capacity of the debris as presented by Conway and Rasmussen (2000). Application of this approach is only appropriate if the temperature data indicate that the system is predominantly conductive and, even in the case of a pure conductive system, the method necessarily introduces numerical errors that can impact the derived values. The sampling strategies used in published applications of this method vary in sensor precision and spatiotemporal temperature sampling strategies, hampering inter-site comparisons of the derived values and their usage at unmeasured sites. To address this, we use synthetic datasets to isolate the numerical errors of the temporal and spatial sampling interval and the precision of sensor temperature and position in recovering known thermal diffusivity values using this method. On the basis of this, we can establish sampling an analytical strategy to minimize the methodological errors. Our results show that increasing temporal and spatial sampling intervals increases (or leads to) truncation errors and systematically underestimates calculated values of thermal diffusivity. The thermistor precision, the shape of the diurnal temperature cycle, the debris thermal diffusivity, and misrepresenting the vertical thermistor position also result in systematic errors that show strong cross-dependencies dependent on signal-to-noise ratio with which spatiotemporal temperature gradients are captured. We provide an interactive analysis tool and best-practice guidelines to help researchers investigate the effect of the sampling interval on calculated sub-debris ice melt and plan future measurement campaigns. These findings can be used to plan optimal field-sampling strategies for future campaigns and as a guide for common reanalysis of existing datasets to allow intercomparison across sites.

1 Introduction

Debris-covered glaciers can be found in tectonically active mountain regions, such as Alaska, the European Alps, High-mountain Asia, or New Zealand (Herreid and Pellicciotti, 2020), where large amounts of debris migrate into the ice via glacial and periglacial processes (Shugar and Clague, 2011; Scherler et al., 2018; Anderson et al., 2018). Debris falling onto the ablation zone contributes directly to any surface debris load, while debris added to the glacier surface in the accumulation zone or sourced subglacially is transported englacially to the ablation area of the glacier, where it melts out and contributes additional debris load (Nicholson and Benn, 2006; Kirkbride and Deline, 2013; Anderson et al., 2018), as shown in Fig. 1a. In comparison to clean ice, thin or patchy debris amplifies ice melt due to its higher absorptivity of short-wave radiation, while thicker debris layers reduce ice melt due to the insulation and attenuation of the diurnal heating signal (Inoue and Yoshida, 1980; Kayastha et al., 2000; Kirkbride and Dugmore, 2003; Mihalcea et al., 2006; Brock et al., 2010; Reznichenko et al., 2010; Fyffe et al., 2014; Minora et al., 2015). The relationship between debris thickness and ablation rate varies for different debris layer compositions and prevailing climatological conditions but retains the same character (Fig. 1b). The critical debris thickness beyond which sub-debris ice ablation is inhibited compared to clean-ice ablation ranges from 15 to 115 mm (Østrem, 1959; Mattson, 1993; Nicholson and Benn, 2006) depending on the optical and thermal properties of the debris and the ambient climate (Inoue and Yoshida, 1980; Nakawo and Takahashi, 1982; Adhikary et al., 1997; Reznichenko et al., 2010). Therefore, in contrast to clean-ice glaciers, where the melt increases towards the glacier tongue in response to typical environmental temperature lapse rates, the spatial pattern of melt of debris-covered glaciers depends more on the debris thickness than on the elevation (e.g. Benn et al., 2012; Rowan et al., 2021; Nicholson et al., 2021). Herreid and Pellicciotti (2020) found that 7.3 % ± 3.3 % of all mountain glacier area is covered by a rock debris cover, which, at a global scale, delays the loss of debris-covered glaciers for the coming decades (Rounce et al., 2023). With continued glacier decline, debris-covered glacier surfaces are expected to increase in absolute and percentage terms in the future (Deline and Orombelli, 2005; Kellerer-Pirklbauer et al., 2008; Quincey and Glasser, 2009; Bhambri et al., 2011; Bolch et al., 2012; Kirkbride and Deline, 2013; Thakuri et al., 2014; Scherler et al., 2018; Tielidze et al., 2020), highlighting the need for accurate modelling of sub-debris ice melt to be included in future glacier projections (Rounce et al., 2015).

Figure 1 (a) Schematic of a debris-covered glacier with debris transport of subglacially sourced rock debris from the release area to the meltout area. The inset shows a classical thermal diffusivity measurement site, consisting of thermistors at several heights between the near surface and the debris–ice interface. (b) Measurements of the so-called Østrem curves for different glaciers show a common pattern of variation in daily melt rate versus the debris depth, with site-specific variations in maximum ablation and the associated debris thickness. Redrawn from Mattson (1993).

Although, under certain circumstances, heat can be transferred through the debris by convection, advection, and radiation, observations (e.g. Conway and Rasmussen, 2000; Nicholson and Benn, 2012) show that the system often, and especially under dry stable meteorological conditions, approximates Fourier's law of conduction, $q=-k\mathrm{\nabla }T$, where q represents the local heat flux density, k represents the thermal conductivity, and T represents the temperature (Fourier, 1955; Cannon, 1984). Consequently, in models of glacier ice melt, the energy supply for ice melt beneath the debris cover is typically treated as if it were heat conduction only (e.g. Reid and Brock, 2010; Fyffe et al., 2014), driven by the surface temperature, the debris thickness, and a value of debris thermal conductivity to be supplied as a model parameter. As a second consequence, Fourier's law of conduction has also been used to derive representative parameter values of effective debris thermal conductivity for horizontally homogeneous debris layers from field observations of spatiotemporal variations in debris temperature. To do this, the one-dimensional heat conduction equation for a homogeneous, isotropic medium (Eq. 1) is used to derive the apparent thermal diffusivity, κ, from the spatiotemporal variation in a vertical profile of temperature measurements by finding the gradient of the regression line between the first derivative of temperature with time and the second derivative of temperature with depth (Conway and Rasmussen, 2000). Effective thermal conductivity k can then be calculated from κ and the volumetric heat capacity of the debris, given by the specific heat capacity c_s and the material density ρ (Eq. 2), including the porosity for a granular material.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial T}{\partial t}}=\mathit{\kappa }{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial x^{\mathrm{2}}}}+\text{const.},\text{(2)}& {\displaystyle }\mathit{\kappa }={\displaystyle \frac{k}{\mathit{\rho }\cdot c_{\mathrm{s}}}}\begin{array}{l}\leftarrow \text{thermal conductivity}\\ \leftarrow \text{heat capacity}\end{array}\end{array}$$

Application of this method therefore requires (1) a vertical array of temperature measurements through the supraglacial debris cover (Fig. 1a) for conditions in which the debris heat transfer closely approximates that of a conductive system from which the apparent κ is derived and (2) an estimate of the volumetric heat capacity of the debris used to convert the apparent κ into effective conductivity.

To meet the first requirement, a sample site must be chosen for which lateral heat transfer can reasonably be expected to be negligible, so a site that is horizontally homogeneous in factors such as slope, debris type, and thickness and without evidence of any hydrological heat transfer. Then, the observed temperatures must be evaluated to find specific time periods and vertical subsets of debris temperature profile data that are identified as being “well-behaved” approximations of a conductive system; data that show evidence of non-conductive processes can be excluded from subsequent analysis (Conway and Rasmussen, 2000). To meet the second requirement, estimates of the debris porosity and the rock density thermal properties must be made. Commonly used values for these terms are porosity of 0.3, rock density of 2700 kg m⁻³, and rock specific heat capacity of 750 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, with a 10 % error applied to the combined terms (Conway and Rasmussen, 2000). Most studies assume that the pore spaces are air-filled when calculating the volumetric heat capacity, but, in principle, if the debris cover is known to be fully saturated, a water-filled case can be used to obtain the volumetric heat capacity of the sampled debris layer (Nicholson and Benn, 2006). In an ideal case, this workflow can yield a reliable estimate of effective thermal conductivity from a homogenous dry portion of the debris with stable meteorological forcing conditions and minimal non-conductive processes. Further use of these effective dry-debris thermal conductivity data in surface energy balance models can allow non-conductive processes and non-uniform debris layers to be included in the model structure by, for example, accounting for stratification in the debris porosity and air flow through the debris (Evatt et al., 2015), stratification of moisture content, and associated phase changes within the debris layer (Collier et al., 2014; Evatt et al., 2015; Giese et al., 2020).

As natural debris covers often show vertical variation in porosity, grain size, and moisture content, recent studies have explored multi-layered applications of the thermal diffusion representation of the debris layer. Laha et al. (2022) perform multiple rather than single regression analysis to account for (i) unknown depth variation in κ in a two-layer model and (ii) non-conductive heat sources/sinks. They apply various methods to synthetic datasets to highlight that applying the original method of Conway and Rasmussen (2000) produces large errors when trying to recover a target κ that varies with depth and that unequally spaced temperature measurements introduce substantial truncation errors. If unequal spacing of measurements cannot be avoided, their new Bayesian method of determining κ outperforms that of Conway and Rasmussen (2000). Petersen et al. (2022) also include a term for depth-varying κ into the heat conduction equation and perform multiple linear regression to solve for its variation with depth in natural debris cover, identifying non-conductive processes as the residual from a comparison of the observed and modelled time-dependent temperature evolution. They find non-negligible heat transfer related to air motion and latent heat fluxes within the debris on Kennicott Glacier. These approaches offer solutions for the potential of vertically varying debris properties and allow quantified assessment of non-conductive processes in measured field sites.

Despite these new developments, the method of Conway and Rasmussen (2000) has been historically widely used (e.g. Nicholson and Benn, 2006; Haidong et al., 2006; Juen et al., 2013; Chand and Kayastha, 2018; Rounce et al., 2015; Rowan et al., 2021) and has provided the majority of published debris thermal conductivity values used in generalized surface energy balance models (Reid and Brock, 2010; Fyffe et al., 2014; Evatt et al., 2015) and for regional intercomparisons of supraglacial debris properties (Fontrodona-Bach et al., 2025). Thus, many studies of debris-covered glaciers rely upon the robustness of debris thermal properties produced following Conway and Rasmussen (2000). The limited number of datasets used to provide generalized values of effective thermal conductivity have deployed very different field and analytical strategies, with temporal and spatial sampling intervals, thermistor placement within the debris, debris depth of the sampled site, and sensor precision all selected ad hoc in different studies and differing from measurement site to measurement site (e.g. Juen et al., 2013; Chand and Kayastha, 2018; Rowan et al., 2021). For example, spatiotemporal sampling intervals range from 2 cm to tens of centimetres and from 5 min to 6 h, sometimes including time-averaged rather than sampled temperatures (Appendix B). The impact of these choices on the derived κ values is not well addressed in the published literature, but, for example, the same data from Imja Glacier in Nepal analysed at 30 min (Rounce et al., 2015) and 60 min (Rowan et al., 2021) intervals yielded thermal conductivity values that differed by almost 40 % despite the same properties being used to derive thermal conductivity from κ. This highlights that baseline literature values that are used in surface energy balance modelling may be differently influenced by sensor, installation, and numerical truncation errors and indicates that care should be taken when comparing across sites for which different instrumental and analytical choices have been made (e.g. Rowan et al., 2021; Miles et al., 2022). Therefore, a deeper exploration of the error sources of this method is warranted, and it would be advantageous to develop standardized field and analytical implementation strategies.

2 Aim of this study

This study explores the effect of measurement setup on κ values derived using the method of Conway and Rasmussen (2000) in order to highlight the potential dependency of published values of thermal conductivity on the spatiotemporal intervals chosen for the analysis and on the sensor precision and locational accuracy. To achieve this, we apply the method of Conway and Rasmussen (2000) to data generated using a forward diffusivity model for a purely conductive system with a specified value of κ and assess how closely the known κ is recovered when varying choices of instrumental and analytical setups. Since the approach recommended by Conway and Rasmussen (2000) is only valid for conductive systems, we focus our study on a purely conductive system to provide a baseline reference for individual method-related error sources, expanding the analysis of the impact of irregular spacings performed in Laha et al. (2022) to include an assessment of a wider range of field measurement choices. By isolating the individual roles of these different error sources, they can be quantified and their tendencies can be understood, thereby making possible a more critical reassessment of the extent to which differences in published effective thermal conductivity values reflect real-world differences in debris properties or instrumental and analytical choices. We provide an interactive tool (https://github.com/calvinbeck/TC-DTD, last access: 30 June 2025) to allow analysis of the combined errors for any given measurement procedure and a best-practice guideline on how to minimize the systematic errors of using this method (Appendix A).

3 Methods

3.1 Artificial data for benchmarking-derived thermal diffusivity

To test the method of Conway and Rasmussen (2000) for different scenarios, we generate synthetic data for debris cover thicknesses of 30 and 100 cm and κ values of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ and $\mathrm{10}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹ to represent a range of values obtained from previous field studies from glaciers across the globe (Laha et al., 2022). The interactive tool allows users to perform analyses for any alternative choice of debris thickness and κ. To generate data for a perfectly conductive system, we force the heat equation with five 10 d (days) surface temperature time series (Fig. 2) and a 0 °C boundary condition for the debris ice interface. The first 2 d of temperature-forcing data is used to initialize the model, and the different debris layer thicknesses are represented by varying the number of vertical grid points in the domain while maintaining equidistant spacing.

Figure 2 Characteristics of the surface temperature forcing for the artificial data generation, which consists of 10 d time series of two analytical sine curves and three experimental temperature measurements within the debris layer. The sine curves have an average temperature of 7.5 °C and the same amplitude. Surface forcing from field data is derived from the uppermost thermistor, which lies 1–5 cm below the surface, as indicated in brackets. Field data 1 and 3 were recorded at Lirung Glacier (Nepal) during September 2013 (5 cm below surface) and April 2014 (1 cm below surface) respectively and were provided by Chand and Kayastha (2018). Field data 2 was recorded at Vernagtferner (Austria) during June 2010 (4 cm below surface) and was provided by Juen et al. (2013). The colour scheme of these forcings is used in subsequent figures.

We use the Crank and Nicolson (1947) method to solve the heat conduction equation for this set of given constraints. This implicit finite-difference method is convergent second-order in time and numerically stable. The method is based on the trapezoidal rule and is a combination of the Euler forward and backward methods in time. For the thermal heat equation, it results in the following equations:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{T_i^{n+\mathrm{1}}-T_i^n}{\mathrm{\Delta }t}}=\\ & {\displaystyle \frac{\mathit{\kappa }}{\mathrm{\Delta }{x}^{\mathrm{2}}}}(T_{i+\mathrm{1}}^{n+\mathrm{1}}-\mathrm{2}{T}_i^{n+\mathrm{1}}+T_{i-\mathrm{1}}^{n+\mathrm{1}})\text{(forward Euler)},\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{T_i^{n+\mathrm{1}}-T_i^n}{\mathrm{\Delta }t}}=\\ & {\displaystyle \frac{\mathit{\kappa }}{\mathrm{\Delta }{x}^{\mathrm{2}}}}(T_{i+\mathrm{1}}^n-\mathrm{2}{T}_i^n+T_{i-\mathrm{1}}^n)\text{(backward Euler)}.\end{array}\end{array}$$

Combining these results in the Crank–Nicolson scheme,

$$\begin{array}{}\text{(5)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{T_i^{n+\mathrm{1}}-T_i^n}{\mathrm{\Delta }t}}=\\ & {\displaystyle \frac{\mathit{\kappa }}{\mathrm{2}\mathrm{\Delta }{x}^{\mathrm{2}}}}\left(\right(T_{i+\mathrm{1}}^{n+\mathrm{1}}-\mathrm{2}{T}_i^{n+\mathrm{1}}+T_{i-\mathrm{1}}^{n+\mathrm{1}})+(T_{i+\mathrm{1}}^n-\mathrm{2}{T}_i^n+T_{i-\mathrm{1}}^n\left)\right).\end{array}\end{array}$$

Because of the implicit nature of the Crank–Nicolson scheme, an algebraic equation or linearizing the equation is necessary to solve the next time step. In our case, we can use the boundary conditions $T(x=\mathrm{0},t)=f\left(t\right)$ and $T(x=D,t)=\mathrm{0}$, where f(t) represents the arbitrary temperature-forcing function (Fig. 2). Although the method is unconditionally numerically stable for the heat equation (Thomas, 2013), unwanted spurious oscillations can occur if the time steps are too long or the spatial resolution is too small. To avoid this, we use the following stability criterion:

$$\begin{array}{}\text{(6)}& \mathit{\kappa }{\displaystyle \frac{\mathrm{d}t}{\mathrm{d}{x}^{\mathrm{2}}}}\le {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

Meeting this criterion (Eq. 6) for both tested values of κ and all five forcing datasets (Fig. 2), the simulated temperatures are produced at 5 min and 2 cm resolution with float-point precision. The resulting generated data (e.g. Fig. 3) provide an ideal reference from which temperatures can be sampled in space and time to replicate field measurements from “well-behaved” portions of vertical temperature profiles within supraglacial debris, meaning subsets of the data that can be shown to closely approximate a conductive system.

Figure 3 Some 5 d examples of the artificially generated debris layer temperature time series data for the skewed sine forcing (a) and the field data 3 forcing (b) for a 30 cm debris layer with κ of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹ using the Crank–Nicolson scheme. (c, d) Daily averaged debris layer temperature profile for the full 10 d time series of the boundary conditions in the upper panels, showing that the often-used steady-state assumption (Evatt et al., 2015) of the daily mean debris layer temperature, shown by a linear temperature gradient, is only fulfilled for periodic daily temperature forcings.

3.2 Experiments performed

We apply the Conway and Rasmussen (2000) method of deriving apparent κ for a selected range of analytical setups as described in the following subsections. When calculating κ from data resampled from the synthetic cases, we calculate a single diffusivity value for the last 8 d of each forcing dataset, although the interactive tool also offers the option to calculate κ at a daily scale for assessment of field datasets. The calculation of the centred spatial derivatives is suitable for unequal grid spacing, but we do not include analysis of unequal vertical thermistor spacings in this study as this was presented in a previous study (Laha et al., 2022). The properties of the analytical setup that are varied are Δt, Δx, varying the precision of the temperature data, and adding Gaussian noise to assess statistical uncertainty. The performance of each experiment at recovering the known κ prescribed in the artificial data is assessed by calculating the relative error:

$$\begin{array}{}\text{(7)}& \text{relative error}={\displaystyle \frac{{\mathit{\kappa }}_{\text{true}}-{\mathit{\kappa }}_{\text{estimated}}}{{\mathit{\kappa }}_{\text{true}}}}.\end{array}$$

Positive relative error values thus correspond to an underestimation of κ compared to the known value. As effects of individual potential sources of error are contingent on other properties of the experimental setup, we present illustrative examples of the error tendencies and their co-dependencies over a range of properties. The full potential parameter space can be explored in the interactive tool. Firstly, the synthetic data are resampled without any added sensor or installation uncertainty to examine the behaviour of numerical truncation errors. Subsequently, the errors associated with the sensor and installation uncertainty are presented.

3.2.1 Quantifying truncation errors in space and time

In theory, the numerical solution to the diffusion problem should be equal to the analytical solution for infinitesimally small spatial and temporal sampling intervals. Truncation errors are expected to scale with the temporal and spatial increment of the analysis with respect to the diurnal forcing cycle (Laha et al., 2022). Higher-order approximations would reduce the truncation error, but errors due to measurement uncertainties would dominate, as described by Zhang and Osterkamp (1995).

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}& \underset{\mathrm{\Delta }t\to \mathrm{0}}{lim}{\displaystyle \frac{T_{t+\mathrm{1}}-T_{t-\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }t}}=\dot{T}\text{and}\\ & \underset{\mathrm{\Delta }x\to \mathrm{0}}{lim}{\displaystyle \frac{T_{x+\mathrm{1}}-T_x+T_{x-\mathrm{1}}}{(\mathrm{\Delta }x{)}^{\mathrm{2}}}}=T^{\prime \prime }\end{array}\end{array}$$

For $\mathrm{\Delta }x,\mathrm{\Delta }t=\mathrm{0}$, the equations are not solvable.

For the temporal truncation error, we resample the artificial data both by skipping and by averaging over an increasing Δt (Fig. 4) from 5 min (the native resolution of the artificial data) to 6 h intervals to encompass the highest- and lowest-resolution temporal sampling of published field data (Appendix B). When skipping, we select every nth value and omit the rest. When averaging, we take the mean temperature over n values. While most studies store samples of the thermistor data at fixed Δt, we include an assessment of this averaging approach, as some published field data collection campaigns are based on measurements of temperatures averaged over Δt (e.g. Rowan et al., 2021). For the spatial truncation error, we resample by skipping data points in space over a range of intervals to decrease the resolution of the 2 cm resolution artificially produced data. For this analysis, we use the highest-resolution temporal forcing with Δt of 5 min and calculate κ for the centre of the debris layer, expanding Δx symmetrically around this point. For assessing truncation errors due to both temporal and spatial resampling, the temperature values are used with their float-point accuracy from the generated data, which implies perfect sensor precision.

Figure 4 Illustrating the two different temporal resampling methods by displaying the temporal grid for different sampling intervals. We compare the method by skipping every nth grid point (a, blue background) or by averaging over n grid points (b, orange background).

3.2.2 Quantifying sensor and installation errors

Thermistors used to record supraglacial debris temperature profiles over time have varying manufacturer-stipulated sensor precision, and there may be uncertainty around their exact location in the debris cover, as this can be challenging to measure with a high degree of accuracy in the field and it can change if the debris moves.

To simulate the effect of temperature measurement precision, we discretize the temperature data to correspond with the measurement precision of 0.1 to 0.4 °C, which is representative of the precision of thermistors typically used in the field. The error properties of these differing sensor precisions are examined for a range of spatiotemporal resampling, in which we ensure symmetrical resampling of Δx by resampling from the centre of the debris layer outwards. Because the observed temperature changes and gradients are smaller at depth, it is expected that a higher precision of temperature measurement is required to capture them. Therefore we also examine how the relative error due to sensor precision varies with the depth in the debris layer at which the analysis is performed. For this, we also consider the potential gain from even higher-precision sensors by including a 0.01 °C temperature discretization, although this is more precise than any of the thermistor properties reported in the literature.

To simulate cases where either the vertical location of the temperature measurement is inaccurate or the thermistor is displaced vertically over time, we use the sampled temperatures at float precision and add a time-invariant vertical offset to each temperature measurement position. Each offset value is randomly sampled from a Gaussian distribution with a standard deviation of 0.5 cm around the true vertical measurement position to represent an inaccurate field measurement of the vertical position. If thermistors move within the debris due to settling or debris migration, the positional inaccuracy could even be larger, but this would likely be discernible from evidence of debris movement or identified when the thermistors were removed from the debris layer, allowing affected data to be excluded from further analysis. For both analyses of the effects of sensor precision and location accuracy, we present only the idealized sinusoidal-forcing data to best isolate the systematic error patterns and how they co-vary with the truncation errors established by the first analysis steps (Sect. 3.2.1).

3.3 Statistical uncertainty estimation

The method of Conway and Rasmussen (2000) is only valid in well-behaved conductive systems; therefore the aim is to only apply the method to a time period and vertical section where this assumption is largely fulfilled. Therefore, our error analysis so far assumes the debris to be a purely conductive, vertically and horizontally homogeneous system, while, in nature, the debris cover will not be perfectly homogeneous and some non-conductive processes are expected to contribute to temperature data even in “well-behaved” sections.

To show that the model-related error sources studied remain relevant despite additional external error sources, we add random statistical noise to the data time series that we perform our analysis on. For this, we use the pure sine curve forcing for a 100 cm thick debris layer, with Δx of 2 cm and Δt of 5 min resolution for a κ of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹. Subsequently, each individual float precision temperature value of the generated temperature time series is modified by a value randomly sampled from a Gaussian distribution with a mean value of 0 °C and a standard deviation of 0.1 °C. This procedure is repeated 20 times to generate a small ensemble of individually perturbed temperature time series. The introduction of this statistical noise of σ_T=0.1 °C does not account for any specific physical processes, since non-conductive processes and effects due to spatial inhomogeneity would produce systematic temperature shifts on a multi-hourly to seasonal timescale, as observed in some field datasets (Conway and Rasmussen, 2000; Nicholson and Benn, 2012; Petersen et al., 2022). The σ_T=0.1 °C is rather selected to statistically perturb the model system and simulate the effect of additional errors. By increasing or decreasing the selected σ_T value, the effect of the perturbation is respectively amplified or attenuated, but the general impact remains the same. The data are analysed as in the previous sections by varying the temporal sampling interval and the vertical position in the debris layer for three selected vertical grid spacings Δx (4, 8, 16 cm) to capture the co-dependencies of the error properties with these measurement choices. The temporal resampling is performed by skipping to preserve the maximum temperature perturbations to illustrate the effects of a maximum perturbation. When resampling by averaging, the perturbed values would equal out for longer temporal averaging periods. For each parameter combination, the mean of κ is calculated from the ensemble with a respective standard deviation to display the value spread.

4 Results

While the interactive tool provided allows a full range of sampling strategies to be explored, here we present results for selected cases within the range of realistic instrumental setups. Our focus is to provide illustrative examples that characterize the error properties of each individual source.

4.1 Error due to temporal truncation

We illustrate the behaviour of the temporal truncation error calculated for a 100 cm thick debris layer with κ of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹ for up to 6 h sampling intervals for both skipping and averaging resampling methods. As few field studies use Δx as small as our 2 cm resolution artificial data, we show an example with Δx of 6 cm to better represent field observations. We show the behaviour at two depths within the debris layer to illustrate the depth dependency of the error behaviour.

The relative error in κ due to temporal truncation error shows a general pattern of monotonic increase with increasing Δt for the skipping method (Fig. 5a and b). Consistently positive relative errors indicate that increasing the temporal sampling interval systematically underestimates κ. At shallow depths, the less sinusoidal the temperature forcing is, the larger the error at all sampling intervals (illustrated by the 8 cm depth cases shown in Fig. 5a and c). At great depths, the error for the sinusoidal forcing remains similar to that in the near surface, while the noisy surface diurnal signals are smoothed at depth and the associated error tends to be more similar to that of the sinusoidal surface forcing (illustrated by the 50 cm depth cases shown in Fig. 5b and d). When data are resampled by averaging, the temporal truncation error is very similar for the sine curve, but, for the noisy-field-forcing data, averaging reduces the error compared to the skipping resampling method (Fig. 5c and d). These patterns of error behaviour are also seen for κ of $\mathrm{10}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹.

Figure 5 Relative temporal truncation error of recovering κ using different temporal sampling intervals: comparison of different temperature forcings for skipping (a, b: blue boundary) and averaging (c, d: orange boundary) resampling methods for two different depths in the 1.0 m debris layer with a target κ of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹.

Considering the maximum relative error produced by typical field installations, we can take the case of calculating diffusivity at a point as close to the surface as is reasonably possible at 4 cm, requiring a thermistor spacing of 2 cm combined with the longer typical time sampling interval of 1 h and calculating over a period with noisy surface forcing. This combination yields a maximum temporal truncation relative error of 25 %. To minimize the error from a truncation perspective, a minimum temporal resolution is desirable, and selecting days with surface temperature forcing that is closer to sinusoidal will decrease errors that may otherwise be significant at shallow depths.

4.2 Error due to spatial truncation

We illustrate the behaviour of the spatial truncation error calculated for a 100 cm thick debris layer with κ of 5 and $\mathrm{10}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹, for Δx up to 50 cm, using a sample of the five surface-forcing datasets at Δt of 5 min.

Spatial truncation error values (Fig. 6) remain quasi-constant for low Δx, up to when the centred differencing scheme spans more than 20 cm, and thereafter increase rapidly with increasing Δx. The spatial truncation error is relatively insensitive to the different surface temperature forcings and, in contrast to the temporal truncation error, does not vary markedly with debris depth. Instead, the κ imposes a strong influence, with higher κ having smaller errors, shifting the respective curves to the right as shown for the case of the sinusoidal forcing in Fig. 6. Given that the diffusivity is the target of sensor installations, this parameter cannot be known in advance, and the results suggest that Δx of below 14 cm is desirable to minimize spatial truncation errors across a range of potential κ. The consistently positive error values mean that the spatial source of truncation error also has the tendency to systematically underestimate κ, increasingly so with more widely spaced temperature measurements.

Figure 6 Comparison of the spatial truncation error for two different κ values and forcing types, calculated for the central position in a 1.0 m debris layer for symmetrically increasing Δx. For clarity, we show only one curve for the higher diffusivity value, as all curves are shifted similarly when varying the target κ. The forcing datasets are at float precision with Δt of 5 min.

4.3 Error due to thermistor precision

To illustrate the role of temperature sensor precision, we firstly focus on the range of sensor spacings that are not affected by the spatial truncation error, i.e. for Δx up to 14 cm (Fig. 7), and show the relative error for a Δt ranging from 5 min to several hours. The error due to temperature discretization is generally less pronounced for smaller temperature discretizations, representing greater thermistor precision. Maximum errors occur for small values of Δx, decreasing to stable relative errors of <25 % for Δx>6 cm, above which the error also decreases systematically with decreasing temporal sampling interval. Values of Δx between the dominant spatial truncation error (Sect. 4.2) and the error due to the sensor precision are desirable, so between ca. 6 and 14 cm for the representative parameter space explored in our analyses.

Figure 7 Relative error of estimated κ due to thermistor temperature discretization of 0.1 and 0.4°C for vertical sampling intervals up to 0.14 m and Δt resampled by skipping (a, b) and averaging (c, d) for the intervals shown in the legend, such that the 5 min dataset is identical for both methods. The case presented is a centred sampling of a 0.3 m thick layer with target κ of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹ forced with a sinusoidal surface forcing.

The depth dependency of the error associated with discretization indicates the importance of high-precision sensors for sampling the debris at depth (Fig. 8). For a Δx of 2 cm, only measurements with a maximum thermistor uncertainty of 0.01 °C would produce correct values and then only for the first 20 cm of debris. Increasing Δx to 6 cm, the relative error decreases for all curves. Still, the thermistors used in most field experiments, which have reported precision ranging from 0.1 to 0.4 °C, would not produce correct values at depth. For the case shown, it would become difficult to obtain reliable values at depths beyond 60 cm even with high-precision thermistors. The error behaviour is dependent on capturing temperature gradients sufficiently well, so the specific error limits are dependent on the amplitude of the surface-forcing fluctuations and diffusivity and the chosen discretization and spatiotemporal sampling. For a given discretization, meaningful values can be obtained at greater depth by enlarging the Δx, but higher-precision sensors are always an advantage. As for both types of truncation error, the sensor precision error systematically underestimates the target κ.

Figure 8 Relative error of estimated κ due to thermistor discretization by depth for a 1.0 m debris layer, with 5 min sinusoidal surface forcing for Δx of (a) 2 cm and (b) 6 cm, with the different coloured lines corresponding to different values of temperature discretization.

4.4 Error due to vertical thermistor position inaccuracy

Conway and Rasmussen (2000) report that a vertical error of 0.5 cm would result in a marginal temperature difference of 0.1 and 0.02 °C for their measurement setups. They and others (e.g. Nicholson and Benn, 2012) interpret this to mean that a vertical thermistor displacement would not affect the results as long as this value does not change in time.

Our analysis, however, shows that low-accuracy knowledge of the temperature measurement location could produce a systematic error for smaller Δx. For example, in the relatively rare case that sensors are installed with a Δx of 2 cm, the resultant error on calculated values of effective thermal diffusivity is so large that the data would become unusable. With increasing Δx, the relative error decreases, such that the mean κ over the depth of the layer recovers the target value. This error source is the only one in this study that has the potential to increase κ values, as shown in Fig. 9 by the spread of κ above the known reference value.

4.5 Statistical uncertainty estimation

In contrast to the noise-free case shown in Fig. 5, with the addition of statistical noise, the relative temporal truncation error now increases with depth, and the predominance of relative errors <25 % in the sub-hourly Δt range can now only be recovered in the near-surface portion of the debris layer (Fig. 10). The standard deviation of the error curves nearer the surface is less than a few percent of the relative error, therefore showing a minimal ensemble spread, while, at depth, the ensemble spread is larger. From this, we can see that, where the random noise introduced is large compared to the spatiotemporal temperature gradients, as is the case at greater depths in the debris layer, the method is essentially no longer applicable. Increasing the Δx decreases the relative error found at depth but has little impact on the smaller errors nearer the surface. At larger Δx, even the near-surface values now have a non-zero relative error for short Δt; this is due to the spatial truncation error of the vertical sampling interval as displayed in Fig. 6 coming into play, while, at greater depth in the debris, the larger Δx decreases the relative error, although this still remains >0.6 with a large relative error and standard deviation values of ∼10 %. In our example, the combination with the most precise recovery of the target κ, with relative error approaching zero, was for Δx=8 cm at an 18 cm depth and at a 5 min temporal sampling interval. For this combination, the relative error due to temporal truncation error increases to ∼10 % and ∼20 % at Δt of 120 and 240 min respectively.

Displaying the noise-induced relative error of κ more explicitly in relation to the depth in the debris over the span of the shared calculation range (0.18–0.82 cm) highlights that there are characteristic transition zones between where the method is still applicable and where it is not, and, as these scale with the relative magnitude of the noise, the transition location is dependent on the Δx used in the analysis (Fig. 11), along with the amplitude of the surface forcing, the diffusivity, and the temperature discretization. For example, for the uppermost section of the artificial debris layer, all curves with 5 and 60 min sampling intervals provide relative errors below 10 %, while, in the data combination we show, the transition to relative error >20 % for these sampling intervals is 0.3, 0.45, and 0.65 m depth for vertical grid spacings of 4, 8, and 16 cm respectively. Therefore, as is the case for the depth dependence of temperature discretization (Fig. 8), increasing the Δx increases the depth at which meaningful values can be recovered when noise is present. However, increasing the grid spacing also results in a Δx truncation error, which is visible in Fig. 11c as a vertical displacement in relative error values additional to the displacement caused by the temporal truncation error.

5 Discussion

In previously published data, most apparent thermal diffusivity derived using the method of Conway and Rasmussen (2000) are below $\mathrm{10}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹, typically ranging from 1 to 30, with some outlier values exceeding $\mathrm{100}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹ (see Table 2 in Laha et al., 2022). The implementation errors that our analysis reveals are often comparable to this range of published values, highlighting how relevant it is to correctly consider the numerical errors in choosing how to apply this method.

While the interactive tool accompanying our analysis allows a wider range of the parameter space to be explored, the cases we present were chosen to characterize the main numerical error sources inherent in the method within the parameter space of published values (Appendix B). The numerical and measurement implementation error sources investigated here all tend to systematically underestimate κ, while the relative error associated with uneven thermistor spacing (tested for a three-thermistor case by Laha et al., 2022) was previously identified to systematically overestimate κ by up to 50 % at thermistor spacing ratios of 1:5.

Figure 9 Illustrating the influence of thermistor displacement on estimated κ by randomly displacing locations by a normal distribution with a standard deviation of 0.5 cm over a range of vertical spacing intervals. The true/target thermal diffusivity is shown by the horizontal black line, showing that, for small temperature sampling intervals, sensor displacement results in large inaccuracies in κ.

In general, the numerical errors associated with applying this method are all related to how well the temperature gradients in space and time within the debris cover can be captured by the instrumental setup. Temporal truncation errors in the absence of statistical noise are typically <25 % in most expected deployment settings at sampling intervals of ≤60 min. Near-surface measurements suffer more error because the diurnal temperature cycle at the surface is the most non-sinusoidal and therefore produces larger temporal truncation errors. Consequently, conditions that more closely approximate sinusoidal conditions (i.e. clear-sky stable atmospheric conditions) reduce the errors in the near-surface layers, but this becomes less relevant at depth, as surface noise introduced by weather is progressively smoothed out at greater depth in the debris. Spatial truncation due to the choice of thermistor spacing is not very sensitive to the non-sinusoidal forcing but becomes ≥25 % at Δx above ∼25 cm for the range of κ reported in the literature, and the error is larger for smaller κ. The Δx range at which errors are small and similar regardless of the forcing and κ is ≤14 cm, providing a conservative upper bound to limit spatial truncation errors. Even though a Δt or Δx→0 would produce a minimal truncation error, sampling intervals that are too small can also produce erroneous results because, for a Δt→0, the linear regression coefficient of determination decreases strongly. In practice, this is not a problem for the temporal sampling, since short temporal sampling intervals can always be resampled afterwards. A more significant problem occurs if low-precision thermistors are positioned too close to each other, especially if the profile comprises only a few thermistors, making it impossible to spatially resample the temperature data. While this effect diminishes to a stable value of relative error ≥25 % for Δx above ∼6 cm, with increasing depth, the thermistors must be further apart, otherwise the thermistor measurement uncertainty dominates the measurement. Therefore, although the highest-precision thermistors should always be chosen if possible, using thermistors with maximum precision becomes even more important at greater depths in the debris layer. The only error source investigated here that has the potential to overestimate κ is that due to inaccurate temperature measurement location. This can happen due to poorly measured positions or due to debris settling after sensor installation if the thermistor profiles are installed on a slope, which is subject to gradual gravitational sliding or reworking. In contrast to Conway and Rasmussen (2000), who showed that a constant error in the thermistor position was not important to the analysis, we find that, at least for very small Δx, the calculated κ does depend on the thermistor positions relative to each other being correctly known and sustained over the measurement period. However, thermistors are typically placed more a than a few centimetres apart; this error source might be expected to have little effect if the κ is calculated at several levels in the debris cover, as the mean value of the location-perturbed cases recovers the target diffusivity. Introducing a statistical noise term highlights the manner in which noise degrades the temperature gradients that the method relies on, particularly at greater depths in the debris, where temperature variations in space and time are small compared to the introduced noise term. Thus, care must also be taken to assess if the method is being applied to portions of the debris layer where the gradients are well captured.

Figure 10 Relative errors of thermal diffusivity of statistically perturbed ensemble data for a 1.0 m debris layer varied by temporal sampling interval for three different depths in the debris layer and three different Δx values (4, 8, 16 cm). The ensemble consists of 20 cases, with each individual temperature value being perturbed by a Gaussian distribution with a standard deviation of σ_T=0.1 °C. The solid line is the mean relative error value, and the shaded background represents the standard deviation of the relative error.

Figure 11 Relative errors of thermal diffusivity of statistically perturbed ensemble data for a 1.0 m debris layer varied by the depth in the debris layer for three different sampling intervals and vertical grid spaces of Δx (4, 8, 16 cm). The ensemble consists of 20 individual runs of each 8 d, with each individual temperature value being perturbed by a Gaussian distribution with a standard deviation of σ_T=0.1 °C. The solid line is the mean relative error value, and the shadowed background represents the standard deviation of the relative error.

In the best-practice guidelines (Appendix A), we address all sources of methodological error discussed in this paper, suggesting optimal implementation strategies for future field studies that wish to deploy these methods of analysing representative thermal conductivity of natural debris layers following the method of Conway and Rasmussen (2000). Our recommendations differ somewhat from those of Laha et al. (2022), as the purpose is different. While Laha et al. (2022) sought to determine the optimal method to determine sub-debris ablation rates directly from temperature sensors using a minimal number of thermistors, we seek to understand the best way to determine a representative κ from which effective thermal conductivity suitable for onward use in generalized surface energy balance models can be derived. For their purpose, they propose to “set the sensor spacing to be one-fifth of the debris thickness at the location”; however, the non-linear nature of the single-error sources presented in this paper indicates that we cannot generalize such statements if the goal is parameter determination rather than direct ablation determination. Furthermore, they stated “the top sensor should be placed approximately at the middle of the debris layer”, as this captures the relevant flux being delivered to the underlying ice. Our analysis indicates that, while it is true that thermistors too close to the surface produce large truncation errors, the same is valid for thermistors that are too deep, as the temperature gradient is too small relative to the thermistor precision. By providing an open-source interactive tool that can be used to explore all the methodological sources of error in implementing the most widely used method of determining κ, we offer a ready-to-use means for determining the field setup that minimizes these numerical methodological errors. The intention is that, prior to a new field deployment, the error response of the expected conditions of debris thickness, surface-forcing amplitude, sensor number, and precision can be explored and the best possible field deployment of sensors can be made.

In addition to the errors related to measurement setup and analysis procedure investigated in this study, non-conductive processes within the debris layer (e.g. rain, phase changes) can also be present (Conway and Rasmussen, 2000; Nicholson and Benn, 2012; Petersen et al., 2022). Unfortunately, it is not always clear in the published literature that the thermal diffusivities and associated thermal conductivity values were derived from optimal conditions sampled within the dataset. The suitability of the sampled debris temperature profiles for determining debris thermal parameters must be carefully evaluated on a case-by-case basis, using meteorological data and closely evaluating the measurements and their gradient functions (Petersen et al., 2022) in order to establish that the data subset represents predominantly conductive conditions, before applying the method of Conway and Rasmussen (2000). Once a suitable effective thermal conductivity is established based on “well-behaved” conditions, these base values can be modified for implementation within a surface energy balance model to account for changes in the pore fluid type to allow simulation of varying wet-/dry-debris conditions (Collier et al., 2014; Giese et al., 2020).

The recently published database of supraglacial debris properties, DebDab v1 (Fontrodona-Bach et al., 2025), reveals that, from the 176 values of debris thermal conductivity, only 33 report an associated uncertainty, and, while 121 include the debris layer thickness, only 23 report on details such as the thermistor depth. To facilitate the intercomparison of these data, it would be valuable to include the temporal sampling used, along with the rock properties and porosity used to convert κ to thermal conductivity. Deeper consideration and potential common reanalyses of these data would require the original thermistor data to be publicly available, which is not always the case. Reanalysing previously published vertical temperature profiles with common resampling strategies, based on the findings of this study, would facilitate intercomparison of κ values, while reanalysis using the methods of Petersen et al. (2022) and/or Laha et al. (2022) might yield more robust and representative global values by providing respectively a more rigorous assessment of non-conductive processes and inclusion of multi-layered thermal properties within the natural debris layers that have been sampled.

6 Conclusion

Conway and Rasmussen (2000) provide a practical method to estimate thermal diffusivity values from a vertical array of thermistors in the supraglacial debris layer, which is applicable for spatially homogenous debris and behaves as a close approximation to a purely conductive system. Although this method has become the standard method for determining effective thermal conductivity to be used in surface energy balance models of sub-debris ice ablation (e.g. Nicholson and Benn, 2006, 2012; Juen et al., 2013; Rounce et al., 2015; Chand and Kayastha, 2018; Rowan et al., 2021), our analysis demonstrates several ways in which the derived κ is sensitive to numerical errors related to instrumental setup and analysis choices, even when solving for a pure-conduction case. The method has regularly been used without considering these error sources, making it difficult to robustly compare published values derived using this method.

To address this, we provide an open-source tool (https://github.com/calvinbeck/TC-DTD, last access: 30 June 2025) where researchers can investigate the combined opportunities and limitations of applying the method by Conway and Rasmussen (2000) to glaciology and beyond. We hope this facilitates more consistent and rigorous experimental design in future field measurements determining debris thermal properties by allowing users to simulate their own artificial data, which most closely approximate their planned field site, and repeat all our analyses presented here with their own artificial or field datasets.

In this paper, we used this tool to provide illustrative examples of the magnitude and tendencies of the systematic errors associated with individual instrumental and analytical choices. Based upon our findings, we provide a set of best-practice guidelines (Appendix A) to minimize systematic errors in applying the method of Conway and Rasmussen (2000). While recent publications highlight limitations of the simplest deployment of the heat diffusion equation in natural debris layers due to the role of non-conductive processes and internal debris stratification (Laha et al., 2022; Petersen et al., 2022), our analysis and best-practice guidelines show the sampling strategies that will yield the best results, provided that the temperatures underpinning the analyses demonstrably sample conditions that closely approximate a homogeneous conductive system. Our analysis also highlights that it is challenging to interpret derived debris thermal properties if the sensor and the analysis system are not reported and accounted for. In the light of this, we encourage more rigorous reporting of implementation strategies and uncertainty in order to facilitate cross-comparison of reported results.

Appendix A: Best-practice guidelines

Our analysis leads us to the following best-practice guidelines to help other researchers to get as much as possible out of their measurements.

Thermistor precision

Use a temperature sensor with the highest possible precision, but not exceeding 0.1 K.

Debris layer thickness

To determine a representative thermal diffusivity from which robust, generally applicable thermal conductivity values can be derived, sampling a minimum of 40 cm but ideally deeper (e.g. 100 cm) debris thickness is advised. The maximum depth that can be meaningfully sampled is limited by the thermistor precision and temperature gradients in the debris layer, which can be simulated beforehand using the tool provided.

Number of thermistors

The method requires at least three thermistors, but more thermistors make it possible to calculate diffusivity values for different depths and therefore make it possible to identify non-conductive processes or other inconsistencies within the debris layer. With only three temperature sensors, it is difficult to assess if the sampled debris meets the requirement of closely approximating a conductive system. A second redundant set of thermistors can also be helpful to rule out measurement errors.

Thermistor installation

A site should be chosen that is not expected to be subject to gravitational reworking or sliding of the debris and where lateral heat fluxes are expected to be minimal. Thermistors should be placed at equal vertical intervals of 8 to 20 cm. Even though the uppermost layer often does not produce ideal results, it can be helpful to place a thermistor at or near the debris surface to provide surface-forcing data. Depending on the depth, the thermal diffusivity, and the temperature gradient of the debris layer, the method produces more significant errors with a greater depth, limiting the depth where it makes sense to place thermistors. The sweet spot can be determined by simulating the debris layer of interest beforehand with model parameters from previous measurements or other estimations.

Thermistor recovery

Thermistors have to be carefully extracted and their vertical positions have to be carefully recorded at the end of the measurement period to make sure that they have not moved in the debris while being deployed. In cases where the thermistors have moved, it might be necessary to discard the dataset. Therefore, mounting thermistors to a thermally insulated rod or set of rods so that their positions are fixed is a valuable approach to eliminate this potential error source.

Temporal sampling interval

One should sample with a temporal resolution as short as possible and then average over a 5 min period. Over such a short period, the temperature is assumed to be nearly constant and therefore not to reduce gradients. By averaging the temperature over a short interval, discretization is reduced.

Measurement duration and conditions

The measurement duration and conditions depend on the scientific objective and seasonality, but at least 1 week of suitable stable meteorological conditions is needed. Therefore, if one has unlucky conditions, a measurement duration of several months could be necessary. A shorter period of predominantly sinusoidal surface temperature forcing, with evidence that non-conductive processes are minimal, is the best way to obtain robust values, so avoiding periods of precipitation, seasonal change, and phase change is advised.

Appendix B: Field measurement overview

Table B1 Overview table of thermal diffusivity field measurement sites. DT: debris layer thickness; SR: sampling rate; AC: thermistor accuracy; TM: number of thermistors. Data from Nicholson (2005), Nicholson and Benn (2012), Juen et al. (2013), Rounce et al. (2015), Chand and Kayastha (2018), and Rowan et al. (2021).