Estimation of ice ablation on a debris-covered glacier from vertical debris-temperature profiles

Abstract

A supraglacial debris layer controls energy transfer to the ice surface and moderates ice ablation on debris-covered glaciers. Measurements of vertical temperature profiles within the debris enables the estimation of thermal diffusivities and sub-debris ablation rates. We have measured the debris-layer temperature profiles at 16 locations on Satopanth Glacier (central Himalaya) during the ablation seasons of 2016 and 2017. Debris temperature profile data are typically analysed using a finite-difference method, assuming that the debris layer is a homogeneous one-dimensional thermal conductor. We introduce three more methods for analysing such data that approximate the debris layer as either a single or a two-layered conductor. We analyse the performance of all four methods using synthetic experiments and by comparing the estimated ablation rates with in situ glaciological observations. Our analysis shows that the temperature measurements obtained at equispaced sensors and analysed with a two-layered model improve the accuracy of the estimated thermal diffusivity and sub-debris ablation rate. The accuracy of the ablation rate estimates is comparable to that of the in situ observations. We argue that measuring the temperature profile is a convenient and reliable method to estimate seasonal to sub-seasonal variations of ablation rates in the thickly debris-covered parts of glaciers.

1. Introduction

Glaciers in the Hindu-Kush Karakoram Himalaya (HKKH) constitute an important freshwater reserve for downstream populations (Azam and others, 2021). The discharge of meltwater is affected by ice ablation rates. About 11% of the glacierised area in HKKH has been estimated to be covered by a layer of rock debris (Herreid and Pellicciotti, 2020). Observations and physical theory suggest that a thin cover of debris can enhance ice ablation rates, while a cover of more than a few centimetres debris thickness can reduce ablation by insulating the ice from solar and atmospheric energy fluxes (Østrem, 1959). Therefore, understanding the recent past or predicting the future of HKKH glaciers requires an accurate understanding of the effects of dynamical supraglacial debris cover on ablation (e.g., Banerjee and Shankar, 2013; Anderson and Anderson, 2016; Banerjee, 2017; Ferguson and Vieli, 2020).

The influence of supraglacial debris on ablation depends on its thickness and thermal properties (Mihalcea and others, 2006; Reid and others, 2012; Fyffe and others, 2020). While field observations of debris properties in HKKH are scarce (Conway and Rasmussen, 2000; Nicholson and Benn, 2013; Rounce and others, 2015; Chand and Kayastha, 2018; Rowan and others, 2021), they indicate that debris properties vary within and between glaciers at a range of scales, and through time (Mihalcea and others, 2008; Nicholson and others, 2018; Shah and others, 2019). Although the physical theory is understood, the actual variability of the debris layer properties is not currently well characterised, and the effects on ablation are poorly constrained (Nicholson and others, 2018). This results in considerable uncertainty in many glaciological applications, such as satellite-based debris-thickness estimation (e.g., Mihalcea and others, 2008; Foster and others, 2012; Schauwecker and others, 2015) and glacio-hydrological modelling (e.g., Fujita and Sakai, 2014; Hagg and others, 2018; Zhang and others, 2019; Steiner and others, 2021).

It is generally accepted that the global population of glaciers is undersampled (Mernild and others, 2013) due to the logistical challenges associated with fieldwork. HKKH has become known colloquially as the ‘third pole’ due to its significance in the global cryosphere (Wester and others, 2019), yet, direct measurements of surface ablation have been reported for fewer than 20 debris-covered glaciers in this region (Winter-Billington and others, 2020). Expanding the network of on-site monitoring stations to include a greater number of debris-covered glaciers in HKKH, and sampling from a wider range of geographies (covering a wider range of elevations, e.g., Wang and others, 2019) is necessary to reduce the uncertainties associated with predictions of glacier change.

The glaciological method of monitoring a network of stakes is generally considered to be the most accurate method to measure the surface ablation on debris-covered glaciers (Cogley and others, 2010). This is a labour-intensive method (Kaser and others, 2003). From our direct experience (Shah and others, 2019), the logistical challenges and human resource requirements for performing sub-seasonal glaciological mass-balance measurement on a debris-covered Himalayan glacier are considerable. Such an exercise requires bimonthly field visits to obtain only ablation and surface-displacement data. The alternate method of estimating sub-debris ablation from debris-temperature profiles recorded at auto-logging temperature sensors, which we focus on in this work, is much less labour-intensive. Only a couple of field visits per year are required to obtain both debris thermal properties and sub-debris ablation at up to daily temporal resolution (Nicholson and Benn, 2013).

Conway and Rasmussen (2000) pioneered a method to compute the effective thermal diffusivity of a debris layer and estimate the sub-debris ablation using observed vertical temperature profiles. This approach, hereinafter referred to as the CR method (Conway and Rasmussen, 2000), has been adapted in several studies in HKKH (e.g., Haidong and others, 2006; Nicholson and Benn, 2013; Chand and Kayastha, 2018; Rowan and others, 2021) and elsewhere (e.g., Nicholson and Benn, 2006; Anderson and others, 2021).

The objective of this study is to improve the accuracy of sub-debris ablation rates estimated using the CR method, so that it is comparable to that of the glaciological method. We use the temperature profiles reported in this paper and previously reported in situ observations of ablation rates (Shah and others, 2019) on Satopanth Glacier. We introduce three new methods, including a Bayesian inversion procedure to analyse the temperature data. To achieve a higher accuracy of the estimated sub-debris ablation, we focus on (a) identifying an optimal sensor spacing to minimise discretisation errors, and (b) incorporating the effect of vertical inhomogeneity in the inversion methods.

2. Study area and field data

Satopanth Glacier (30.75°N, 79.36°E, central Himalaya) is in the Garhwal region of India (Fig. 1). The glacier is more than 12 km long with a total surface area of ~19  km², and it spans an elevation range of 3900–6200 m a.s.l (Shah and others, 2019). The glacier valley has steep and high walls, from which large volumes of loose debris are transported onto the glacier by frequent avalanches and rockfalls (Banerjee and Bilal, 2018). Avalanches contribute a major fraction of the total accumulation of the glacier (Laha and others, 2017). Rock debris covers the surface where the average slope is around 7° or less, which is ~60% of the total surface area and most of the ablation zone. The debris cover starts around 4500–4700 m a.s.l elevation and thickens downglacier. The greatest measured debris thickness is 1.27 m, near the terminus. However, it is not an upper bound; for example, boulders with several metres in size are commonly observed near the terminus. It has been observed that locally the debris thickness can change substantially. At a given elevation, the standard deviation of the debris thickness can be $\sim 50\percnt$ of the corresponding mean value (Shah and others, 2019).

Fig. 1. Map of Satopanth Glacier (central Himalaya) showing the locations of debris pits (solid yellow circles) and ablation stakes (solid blue circles), where the data used in this study were collected. The blue and red solid lines denote the boundaries of the glacier and the debris-covered area, respectively. The inset map is the political boundary of India (solid black line) as per the Survey of India, with a solid red circle indicating the location of Satopanth Glacier.

The bedrock and supraglacial debris are calc-silicate, predominantly schist with lesser amounts of biotite-gneiss and granite with tourmaline (Valdiya, 1999). The debris-covered ablation zone is punctuated with ephemeral supraglacial ponds and ice cliffs of various sizes, which is a common characteristic of such low-sloping and extensively debris-covered glaciers (Sakai and Fujita, 2010).

2.1. Vertical temperature profile data

The field measurements of debris temperature profiles were performed during the ablation season (May–October) in 2016 and 2017. The temperature profiles were measured in pits dug into the debris at 11 locations in the lower ablation zone, from 3800 to 4000 m a.s.l., and five locations in the middle ablation zone, from 4100 to 4400 m a.s.l (Figs 1, S1). Debris temperatures were recorded at 15–60 min intervals using HOBO Onset TMC6-HD Water/Soil thermistors (accuracy 0.2°C) at 3–8 depths. The thermistors were connected to HOBO Onset U12-4 External Channel Data Loggers at the debris surface (Supplementary Fig. S2). The thermistors were inserted into the wall of each pit, and the pit was back-filled. The length of the records varied from 7 d to more than a year depending on location, with some data gaps. An example of a year-long temperature record with some data gaps is shown in Figure 2. All the other temperature records from different pits are presented in Supplementary Figure S3. The debris thickness at these 16 pits ranged from 0.22 to 0.77 m. Depending on the debris thickness, 2–3 d were required for the temperature profile to recover from the initial thermal perturbation during the installation procedure. The locations of the pits, the measurement periods and the sensor depths are given in Table 1, and the pit locations are mapped in Figures 1 and S1. There were no ice-cliffs and/or ponds within ~50 m of any of the debris pits. The thermistors were installed with different spacing in each of the pits, at irregular and regular intervals, so that the effect of sensor spacing on the discretisation errors could be evaluated.

Fig. 2. The temperature time series recorded at the depth of 4 cm (solid purple line), 16 cm (solid blue line), and 28 cm (solid red line) at pit SBP1 (see Table 1 for details). The light grey shading denotes the availability of the ablation stake data that is used for ablation comparison.

Table 1. Details of the debris temperature measurements at Satopanth Glacier are given here

Pits with data gaps are marked with a ‘*’.

The temperature profile data from each pit were split into several 7–15 d long time series to facilitate a comparison between the sub-debris ablation estimated from the temperature records and the available ablation stake data during 2016–2017 (Shah and others, 2019). In total, we had 64 temperature time series from 16 pits, each of which was analysed using the methods discussed in Section 4.1.

2.2. Glaciological ablation data

The glaciological ablation data for the ablation seasons of 2015–2017 used in this study werecontributed by Shah and others (2019). The measurements were performed using a network of up to 83 bamboo stakes spanning an elevation range of 3900–4700 m a.s.l (Fig. 1). The debris thickness at the stake locations varied from 0.05 to 1.27 m. The ablation values at the stakes were recorded once in about every 15 d, with some data gaps. Shah and others (2019) showed that sub-debris ablation on Satopanth Glacier is more sensitive to debris thickness than elevation.

3. Theory of one-dimensional vertical heat conduction through a debris layer

On debris-free glaciers, atmospheric energy fluxes are directly incident at the ice surface. In the presence of a debris layer, atmospheric energy fluxes are incident at the debris surface. A positive surface energy balance results in warming of the debris layer. Once the temperature of the debris has been raised above 0°C, a temperature gradient from the debris surface to the ice surface produces a heat flux towards the ice that can result in ablation. Assuming that this is the dominant mechanism of energy transfer, it is possible to estimate the ablation rate by measuring the temperature gradient within the debris layer.

Supraglacial debris is manifestly an inhomogeneous medium, and the interface between the debris and the atmosphere is uneven. Consequently, there could be temperature gradients and heat flow in the horizontal directions (Evatt and others, 2015) due to the lateral variation of debris-layer thermal properties, as well as the debris thickness and the debris-surface temperature. However, past work (Conway and Rasmussen, 2000) suggests that horizontal heat fluxes are not large relative to the vertical heat flux. Hence, a one-dimensional heat equation seems to be sufficient to represent the flux of heat used in ice ablation with reasonable accuracy.

The conceptual mathematical model we use is the one-dimensional heat equation in a vertically inhomogeneous medium. We include an inhomogeneous source/sink term to account for physical sources/sinks as well as the effects of lateral heat fluxes.

(1)$$\eqalign{\partial_{\rm t}{T( z,\; t) } & = {1\over \rho_{\rm d} ( z) C( z) ( 1-\phi( z) ) } \partial_{\rm z} \left(K( z) \, \partial_{\rm z} T( z,\; t) \right)\cr & \quad + s( z) .}$$

Here, z is the vertical distance, measured from a convenient depth within the debris layer (Fig. 3). T(z,  t) is the debris temperature at vertical distance z and time t. ρ _d(z), C(z) and ϕ(z) are the density, heat capacity and the porosity of the debris, respectively, at vertical distance z. K(z) is the thermal conductivity of the debris. The term s(z) accounts for any source or sink of heat that could arise from processes such as condensation/evaporation, convection and horizontal conduction. Therefore, the strength of the term s(z) is a measure of the accuracy of the one-dimensional model (Eqn (1)).

Fig. 3. Schematic of a debris pit with debris thickness d. The vertical positions of the temperature sensors (z = −dz ₁, z = 0 and z = dz ₂) are indicated with solid black arrows. Figures (a) and (b) show the thickness l ₁ (l ₂) of the top (bottom) layer as used in the CR_i, and MC_i methods, respectively (see Section 4.1 for details).

We simplify the model by approximating the inhomogenous debris layer as a layered thermal conductor. In particular, we assume (a) ρ _d = 2700 kg m⁻³, C = 750 J kg⁻¹ K⁻¹ and ϕ=0.3 are constant in z (Conway and Rasmussen, 2000) and (b) K(z) and s(z) are constant in each layer. Then, for each layer, Eqn (1) can be written as,

(2)$$\partial_{\rm t}T( z,\; t) = \kappa\partial^2_{\rm z}T( z,\; t) + s.$$

Here, the thermal diffusivity, κ, is defined as ${K\over \rho _{\rm d} C ( 1-\phi ) }$. Equation (2) can be solved for κ(z) and s(z) using debris temperature profile time series, with the boundary conditions between layers being specified by the continuity of the temperature and heat flux across the interface. The sub-debris ablation rate is calculated using the estimated values of κ, and the following equation:

(3)$$M = {\rho_{\rm d} C ( 1-\phi) \over L_{\rm f} \rho_{\rm w}} \kappa {dT\over dz} .$$

Here, the temperature gradient at the debris-ice interface is ${dT\over dz}$, the latent heat of fusion L _f = 3.34 × 10⁵ J kg⁻¹ and the density of water ρ _w = 1000 kg m⁻³.

4. Methods

In this section, we review the CR method and describe the three new methods proposed by us to infer the values of κ and s from the temperature profile data. Next, we describe the synthetic experiments that we used to assess the accuracy of the values of κ and s inferred by four methods. We then describe the procedures we used to estimate sub-debris ablation rates on Satopanth Glacier from the observed temperature profiles for these four methods. Finally, the methods used to evaluate the accuracy of the estimated ablation rates using the glaciological ablation data are detailed, along with the methods used to estimate the corresponding uncertainties.

4.1. Analysing vertical temperature profiles using the one-dimensional heat equation

Our first new method is a simple generalisation of the CR method (Conway and Rasmussen, 2000) to a two-layered conductor. We denote the original one-layered model as CR_h, and the generalisation to the two-layered model as CR_i. The other two methods are based on a Bayesian Monte-Carlo approach, using the heat equation detailed above as the so-called forward model. We denote the method applied to the one-layered and two-layered models as MC_h and MC_i, respectively.

Since we were mainly interested in estimating the sub-debris ablation rates, we concentrated on the lower part of the debris layer, namely, data from the three bottom-most sensors in each pit (Table I). Hereinafter the top sensor refers to the sensor position z = –dz ₁ (Fig. 3). Thus, the computational domain of our model was from the debris-ice interface to the top sensor rather than the debris surface (Fig. 3). The temperature data from deeper sensors are typically less noisy compared to the data from the shallower sensors (Collier and others, 2014; Giese and others, 2020), which helps to constrain the uncertainty in the estimates.

4.1.1. Finite-difference method

Assuming that the temperature time series at three vertical positions (z = −dz ₁, z = 0 and z = dz ₂ as shown in Fig. 3) are available, the second-derivative term in the right-hand side of Eqn (2) can be approximated as follows,

(4)$$\eqalign{\partial ^2_{\rm z} T ( 0,\; t) & \approx { {T( -dz_1,\; t) -T( 0,\; t) \over dz_1} - {T( 0,\; t) -T( dz_2,\; t) \over dz_2} \over {dz_1 + dz_2 \over 2}} \cr & \quad + {dz_1-dz_2\over 3} \partial^3_{\rm z} T( 0,\; t) + O( dz^3) + \ldots}$$

The lowest order error term is linear in (dz ₁ − dz ₂) when dz ₁ ≠ dz ₂, and O(dz ³) otherwise. This suggests that in this finite-difference scheme, the errors are larger in the case of non-uniform sensor spacing, e.g., dz ₁ ≠ dz ₂.

4.1.1.1. Homogeneous CR method (CR_h)

Given the discrete time series of temperature measured at the three different sensors (as shown in Fig. 3), this method (Zhang and Osterkamp, 1995; Conway and Rasmussen, 2000) employs the following finite-difference approximations to evaluate the derivative terms in Eqn (2) assuming higher-order correction terms (Eqn (4)) to be negligible:

(5)$$\partial_{\rm t}T( 0,\; t) \approx {T( 0,\; t + \Delta t) -T( 0,\; t) \over \Delta t}$$

(6)$$\partial ^2_{\rm z} T ( 0,\; t) \approx { {T( -dz_1,\; t) -T( 0,\; t) \over dz_1} - {T( 0,\; t) -T( dz_2,\; t) \over dz_2} \over {dz_1 + dz_2 \over 2}}.$$

It follows from Eqns (2) that the values of ∂_tT and $\partial _{\rm z}^2\, T$ computed using Eqn (5) and (6) should be linearly related to the measurement uncertainties and the discretisation errors. Therefore, the slope and intercept of the best-fit straight line to these data points were used to calculate κ and s, respectively (Conway and Rasmussen, 2000). The standard errors of the fitted parameters were taken as the corresponding uncertainties.

4.1.1.2. Inhomogeneous CR method (CR_i)

To account for inhomogeneities in the debris thermal properties in a simple way, we propose the following generalisation of the discretisation scheme used in the CR_h method (Eqn (6)) to a two-layered model:

(7)$$\eqalign{{T( 0,\; t + \Delta t) -T( 0,\; t) \over \Delta t} & \approx {{ \kappa_1 {T( -dz_1,\; t) -T( 0,\; t) \over dz_1} - \kappa_2 {T( 0,\; t) -T( dz_2,\; t) \over dz_{2}} }\over {dz_1 + dz_2 \over 2}} \cr & \quad + s,\; }$$

where κ ₁ and κ ₂ are the thermal diffusivities of the top and bottom layers with thicknesses dz ₁ and dz ₂, respectively (Fig. 3a). The finite-difference terms in Eqn (7) were computed using the observed temperature time series from Satopanth Glacier, and we used a three-parameter linear regression to obtain the best-fit values of κ ₁,  κ ₂ and s. The uncertainties in these parameters were computed in the same manner as described in the CR_h method.

4.1.2. Bayesian inversion method

The Bayesian approach to infer the values of κ and s, which we collectively denote as m, from the observed temperature profiles, which we denote by y, is based on Bayes's theorem. The theorem states that the probability of the parameters being m, given the observations y, is given by

(8)$$p( m\vert y) = {\,p( y\vert m) p ( m) \over p( y) },\; $$

where p(y|m) is the probability of observing y given that the parameters are m. We estimated this probability using our forward model as detailed below. p(m) is the apriori probability of the parameters being m. We assume p(m) to be a uniform distribution over a specified range of the parameters. p(y) is an unimportant normalisation constant (Gelman and others, 2014). We describe our implementation of this method for the one-layer case (MC_h) and the two-layer case (MC_i) below.

4.1.2.1. Homogeneous Bayesian method (MC_h)

Vertical heat conduction through a homogeneous layer between the upper sensor and the debris-ice interface was simulated by solving Eqn (2) numerically with an explicit Forward-in-Time-Central-in-Space finite-difference scheme (Slingerland and Lee, 2011). The upper boundary condition was fixed using the observed temperature at z = −dz ₁ (Fig. 3). The bottom boundary at the debris-ice interface was assumed to be at 0°C. The spatial and temporal grid sizes were 0.01 m and 1.0 s, respectively. As the observed temperature data had an hourly to sub-hourly temporal resolution (Table 1), linearly interpolated values of the temperature data from the top sensor were used to set the upper boundary condition. The temperature data for the first day of the simulation were repeated seven times for model spin-up. The modelled temperatures (T ^mod) at z = 0 and z = dz ₂ were used to compute the sum of squared errors (δ²) relative to the observed temperatures (T ^obs) as follows.

(9)$$\eqalign{\delta^2 & = {1\over N} \sum_{t} [ \{ T^{{\rm obs}} ( 0,\; t) - T^{{\rm mod}} ( 0,\; t) \} ^2 \cr & \quad + \{ T^{{\rm obs}} ( dz_2,\; t) - T^{{\rm mod}} ( dz_2,\; t) \} ^2 ] ,\; }$$

where t denotes the time step, and N is a normalisation factor that counts the total number of data points being fitted.

We assumed that p(y|m) ~ exp( − δ²), and that p(m) is a uniform distribution over a given range of the parameters. The corresponding ranges for the parameters κ and s were 0.01 to 10 mm² s⁻¹, and −6 × 10⁻⁴ to 6 × 10⁻⁴ Ks⁻¹, respectively. We searched the two-dimensional space of the parameters (κ,  s) using a Monte Carlo procedure in order to minimise δ²(Eqn (9)). The steps involved in this stochastic minimisation procedure are presented in more detail in the Supplementary Material (Section S1.1 ).

4.1.2.2. Inhomogeneous Bayesian method (MC_i)

This method used two layers, with the boundary being equidistant from the sensors at z = 0 and z = dz ₂ (Fig. 3b). Here, the forward model was solved imposing the continuity of the temperature and the heat flux at the interface between the two layers. In this case m consisted of four parameters, κ ₁,  κ ₂,  s ₁ and s ₂, where the subscripts 1 and 2 refer to the top and bottom layers, respectively. Apart from these differences, the same procedures and the apriori distributions as used in the MC_h were used for the MC_i.

4.2. Uncertainty in thermal diffusivity due to temperature measurement error

The uncertainties in the estimated values of κ due to the measurement errors in the debris temperature profiles were computed using a separate set of Monte Carlo simulations for both the finite-difference and Bayesian methods. We added a zero-mean random Gaussian noise to the temperature time series and to the depths of the temperature sensors. The temperature noise had a standard deviation 0.2°C, which was the specified accuracy of the thermistors. The standard deviation of the noise in the position of the temperature sensors was assumed to be 2 cm. The data with added noise were used as inputs to the four methods to obtain debris thermal diffusivities and ablation rates. This procedure was repeated 100 times, for each method, and 100 copies of the noisy dataset were used to compute the uncertainties of the best fit values of κ and s.

The above analysis was done for two randomly chosen pits, SBP4 and SBP7. We found that uncertainties in both κ and s were similar for the two datasets. Further, we checked that the uncertainties in each parameter (κ and s) did not change substantially if we increased the number of noisy datasets from 100 to 150. The mean percentage error in both κ and s obtained from these computations were assumed to be the same for all the debris pits. Finally, the uncertainties in both κ and s due to measurement errors estimated by the above procedure (σ ₀) were combined with the fit uncertainties (σ ₁) from each method discussed above, to compute total uncertainties ($\sqrt {\sigma _0^2 + \sigma _1^2}$).

4.3. Synthetic experiments to check the accuracy of the methods

We performed two synthetic experiments to check the affect of (1) spacing between sensors, and (2) vertical inhomogeneity in κ, on the accuracy of four methods used to analyse the temperature profiles. Equation (2) was applied in a forward simulation as described in Section 4.1, using specific values of κ. For simplicity, we assumed s = 0 in both the experiments. The time series of temperature recorded by the top sensor in pit SBP6 was used to simulate hourly time series of temperature at vertical distances z = 0 and z = dz ₂ (Fig. 3).

4.3.1. Experiment 1

This experiment was used to evaluate the accuracy of the values of κ estimated using the four different methods described in Section 4.1 for different spacings between the sensors. The heat equation (Eqn (2)) was run as a forward model to simulate the debris temperatures with known values of κ, s, dz ₁ and dz ₂. The value of κ was held constant at 1.0  mm² s⁻¹, and the value of dz ₂/dz ₁ was varied within the range 1.0–5.7, as reported in the relevant literature (Conway and Rasmussen, 2000; Nicholson and Benn, 2006; Brock and others, 2010; Reid and others, 2012; Nicholson and Benn, 2013; Rounce and others, 2015; Chand and Kayastha, 2018; Rowan and others, 2021). In our field measurements, the values of dz ₂/dz ₁ varied in the range 1.0–1.5 (Table 1).

The four methods were applied to the synthetic temperature time series to infer the value of κ. For the two-layered methods (CR_i and MC_i), we computed the effective thermal diffusivity (κ _eff) using

(10)$${d\over \kappa_{\rm{eff}}} = {l_1\over \kappa_1} + {l_2\over \kappa_2},\; $$

where l ₁ and l ₂ are the thickness, and κ ₁ and κ ₂ are the thermal diffusivity of the top and bottom layers, respectively. To estimate the influence of the relative separation between sensors on the accuracy of the corresponding method, we compared the known value of κ with the inferred values of κ _eff for different sets of dz ₁ and dz ₂, and assessed the trend of κ _eff with dz ₂/dz ₁ by a visual inspection.

For an independent validation of the results of Experiment 1, we applied the above procedure to debris temperature data that we copied from a study by Rowan and others (2021). Rowan and others (2021) data were collected on Khumbu Glacier, eastern Himalaya, Nepal, using the same filed method that we used to collect data at Satopanth Glacier. From Rowan and others (2021) data, we considered the pits where the debris temperature had been measured at more than seven different depths during summer, and split the time series from each pit into 15 d windows for consistency with our methodology. For each record, we applied the CR_h method to estimate κ for different values of dz ₂/dz ₁ by selecting three different sensors.

4.3.2. Experiment 2

This experiment was designed to check the accuracy of the methods (Section 4.1) in the case of a vertically inhomogeneous distribution of κ. We kept dz ₂/dz ₁ = 1, s = 0, and used two horizontally homogeneous layers with different κ (κ ₁ for the top layer and κ ₂ for the bottom layer) to generate the temperature data at two bottom sensors in a forward model simulation. Then, the assumed values of κ _eff were compared with the estimated values of κ _eff to assess the accuracy of the inferred κ. The above process was repeated for different sets of κ ₁ and κ ₂ in the forward model simulations.

To ensure that the input data did not bias the results, we repeated both experiments using input data from a different pit, SBP3. The results were the same in substance, therefore the results from pit SBP3 are not discussed further.

4.4. Estimation of sub-debris ablation

The value of the thermal diffusivity of the debris for each temperature record was estimated using each of the four methods described in Section 4.1. Given the estimated values of κ, sub-debris ablation rates were calculated using Eqn (3). In line with previous studies, the assumed values of ρ _d, C and ϕ were assumed to have 10% uncertainties (Macfarlane and others, 1992; Nicholson and Benn, 2013). In both CR_h and CR_i methods, following Conway and Rasmussen (2000), the value of ${dT\over dz}$ was obtained by a linear fit to the mean temperature at all three sensors. In the MC_h and MC_i methods, ${dT\over dz}$ was obtained by a linear fit to the mean simulated temperature at each of the grid points between z = dz ₂ and the debris-ice interface (Fig. 3). In the CR_i and MC_i methods, the thermal diffusivity estimated for the bottom-most layer, κ ₂, was used to compute the sub-debris ablation rates.

The accuracy of the estimated sub-debris ablation rates obtained from the observed temperature profiles was evaluated by a comparison with the glaciological ablation data. To do the comparison, we needed observed ablation for the thickness of the debris at each of the pits. Because of random local variations of debris thickness, we typically did not have ablation stakes installed with exactly the same debris thickness as that of the nearest pit. Following Shah and others (2019), we considered all the observations in the glaciological ablation dataset from across the debris-covered ablation zone for each observation period and interpolated the observed ablation rates as a function of debris thickness to estimate the glaciological ablation rates corresponding to the thickness of the debris at each of the pits. The number of available stakes for any given period varied from 12 to 65. For each of the glaciological ablation rate values, the corresponding uncertainty was computed by the procedure described in Shah and others (2019).

5. Results and discussion

Below, we first present the results of our synthetic experiments and discuss them. Subsequently, we present and discuss the results of our analysis of the temperature profiles observed on Satopanth Glacier, to infer the thermal properties of the debris layer. Next, our estimates of sub-debris ablation rates inferred from the debris-temperature profiles are reported and compared with the glaciological ablation measurements. Finally, we provide a set of recommendations for the experimental design and the numerical method to improve the accuracy of estimated thermal diffusivity and sub-debris ablation rates.

5.1. Factors affecting the accuracy of the methods

5.1.1. Optimal sensor spacing

Synthetic experiment 1 showed that all four methods inferred the input value of κ _eff accurately for dz ₂/dz ₁ = 1 (Fig. 4a). In fact, we found systematic biases in the values of κ _eff estimated using all four methods, which increased as values of dz ₂/dz ₁ increased from 1. Also, using the observed temperature data from Satopanth Glacier, the corresponding values of κ _eff from all four methods were biased for dz ₂/dz ₁ > 1 (Figs 5a–5d). The trend of increasing κ _eff with increasing dz ₂/dz ₁ > 1 from pit data was consistent with the results from synthetic experiment 1. This confirmed the robustness of the finding from synthetic experiment 1, that equispaced sensors produced the most accurate estimates of thermal diffusivity. The trend of increasing κ _eff with increasing dz ₂/dz ₁ > 1 was also observed in the debris temperature data from Khumbu Glacier, which were analysed using the CR_h method (Figs 5f, S5). The true values of κ _eff at Satopanth and Khumbu glaciers are, of course, not known, so the accuracy of κ _eff estimated for those sites cannot be assessed explicitly.

Fig. 4. (a) Synthetic experiment 1: The inferred κ _eff values from four methods are plotted for different values of dz ₂/dz ₁ (see the text for details). The grey horizontal line denotes the true value of κ _eff (1 mm² s⁻¹) used in the synthetic experiments. The vertical grey band is the range of dz ₂/dz ₁ in the field experiments reported in this study (Table 1). (b) Synthetic experiment 2: the fractional errors in the inferred κ _eff relative to the forward model values are plotted against κ ₂/κ ₁ used in the forward model. Here we kept dz ₂/dz ₁ = 1.

Fig. 5. (a)–(d) Estimated values of κ _eff from four methods are plotted as a function of dz ₂/dz ₁ for all the records. The colours of the symbols denote the month of the temperature measurement. (e) The CR_h method estimates of κ _eff plotted for different dz ₂/dz ₁ in a pit (KH1) from Khumbu Glacier during ablation season of 2014 (Rowan and others, 2021).

While the above dependence of the biases in inferred κ _eff on dz ₂/dz ₁ can be qualitatively inferred from Eqn (4), in the sense that the larger the truncation errors, the larger the error in κ _eff, the sign of the bias cannot be determined without knowing that of $\partial ^3_z T( 0,\; \, t)$. With only three sensors used in the methods of analysing debris temperature profiles (Section 4.1), one could not determine the sign of $\partial ^3_z T( 0,\; \, t)$ using finite difference methods. However, if we consider the known analytical solution for the case of an infinite slab with a sinusoidal temperature variation applied on the upper boundary (Anderson and Anderson, 2010), it is seen that $\partial ^3_z T( 0,\; \, t)$ is positive for all t. This is consistent with the positive bias found in synthetic experiment 1 (Fig. 4a).

5.1.2. The choice of method

Synthetic experiment 2 showed that the accuracy of inferred κ _eff depended on the magnitude of the ratio κ ₂/κ ₁, which is a measure of the inhomogeneity of the two layers (Fig. 4b); again, this result held true for all the methods we tested. The root-mean-squared-error (RMSE) of forward model κ _eff relative to that of the CR_h, CR_i, MC_h and MC_i methods were 0.62, 0.08, 0.07 and 0.03 mm² s⁻¹, respectively. As these values show, the accuracy of the CR_h method was substantially lower than that of the other methods. The performance of the CR_i, MC_h and MC_i methods was comparable, but ultimately, the MC_i method outperformed the others, with the smallest RMSE. The accuracy of both the CR and MC methods highly improved in the two-layered model compared to the one-layered one. These results highlight the importance of using a method that accounts for inhomogeneity in the thermal properties of the debris layer when estimating sub-debris ablation rates.

5.2. Thermal diffusivity of the debris on Satopanth Glacier

Based on the results presented in Section 5.1, hereinafter we only consider the temperature records where dz ₁ = dz ₂. This criteria leaves 38 temperature records out of 64. In the field experiments, we found that the debris thickness in four records from the pit SBP1 (for ablation season 2016) had changed after installation, which may be due to debris movement. These were discarded, leaving 34 records. In the CR_h and CR_i methods, four records resulted in unphysical negative values of κ _eff; therefore, these were also discarded. Hence, in the rest of the paper, we discuss only 34 records with respect to the Bayesian methods (MC_h and MC_i), and 30 records with respect to the finite-difference methods (CR_h and CR_i).

In the MC_i method, the mean, standard deviation and range of κ ₁ (κ ₂) were 1.2 (2.0), 0.6 (1.0), 0.3–0.32 (0.7–4.3) mm² s⁻¹, respectively. In the other three methods, the corresponding values and fit qualities are given in Table 3. The uncertainty in the fitted κ values reported in previous studies ranged from 10 to 34% (Conway and Rasmussen, 2000; Nicholson and Benn, 2013; Rounce and others, 2015; Chand and Kayastha, 2018). In this study, all four deployed methods had a median uncertainty of 11%. The magnitude of the difference between the values of κ estimated for each layer using the two inhomogeneous methods (CR_i and MC_i) was similar (Table 3). The range of κ _eff was similar for all the methods, with somewhat larger values compared to the previous reports from glaciers in and outside the Himalaya (Table 2). However, the corresponding median κ _eff in all the methods were similar to that reported in the literature. We speculate that the larger values of thermal diffusivity in a few pits at Satopanth Glacier may be related to changes in moisture content, latent heat fluxes, ϕ, C or ρ _d during the study period. However, additional in situ measurements of K, C, ϕ, ρ _d, moisture content, horizontal and vertical heat fluxes are needed to ascertain that, which will be taken up in the future.

Table 2. Comparison of estimated κ _eff from this study with the values from the literature, within and outside the Himalaya

Here, the range of κ _eff obtained from the pits at Khumbu Glacier is marked with a ‘*’.

Table 3. Details of the thermal properties and heat source/sink terms obtained using the four different methods for the selected pits where dz ₂ = dz ₁, and the goodness-of-fit metrics corresponding to each method

5.2.1. The spatial variability of the estimated thermal diffusivities

We found substantial spatial variability in κ within the uncertainty, and using all methods (Supplementary Figs S6, S7). The mean (standard deviation) of κ _eff estimated using the CR_i and MC_i methods, and that from the literature were 1.3 (0.8), 1.5 (0.9) and 0.9 (0.5) mm² s⁻¹, respectively. Thus, there was a greater within-glacier variability of κ _eff at Satopanth Glacier than reported in the literature. This spread in κ, as well as the vertical inhomogeneity, are important because they imply that debris thickness estimates made using remote-sensing data using the average value of κ may be more uncertain than they have been reported to be (Mihalcea and others, 2008; Foster and others, 2012; Schauwecker and others, 2015). The same may be true of the predictions of sub-debris ablation or meltwater runoff that have been made using energy-balance models (Nicholson and Benn, 2006; Reid and Brock, 2010; Lejeune and others, 2013) or glacio-hydrological models (Fujita and Sakai, 2014; Hagg and others, 2018; Zhang and others, 2019), where a spatially uniform thermal diffusivity has been assumed. The results presented here may help improve the estimates of the prediction uncertainties in future studies.

5.2.2. The temporal variability of thermal diffusivities

A seasonal trend was evident in the estimates made using data of at least a few months in length, in which κ increased after the onset of the monsoon (Supplementary Figs S6, S7). We found the same result using the data from Khumbu Glacier (Rowan and others, 2021) (Figs 5f, S5). However, the observed seasonal amplitude was smaller than the estimated uncertainties in all cases. This implies that seasonal variability in κ is not likely to be a major source of uncertainty in models of debris thickness or ablation.

5.3. Heat sources within the debris layer

The mean, standard deviation and the corresponding ranges of s obtained using all four methods are given in Table 3. The ranges of s ₁ and s ₂ obtained using the MC_i method were close to the range of s used in the corresponding apriori distributions (Section 4.1). One can see in Supplementary Figures S8a and S8b that the δ² of the fits were not sensitive to s ₁ and s ₂, while it had a sharp minima in the κ ₁ and κ ₂ plane. We also found that the MC_i method is not sensitive enough to estimate vertical inhomogeneities in s. This was indicated by an anti-correlation between the observed values of s ₁ and s ₂, so that s ₁ + s ₂ ≈ 0 at almost all the points. This implied that a similar δ ² can be obtained by simply exchanging the sign of s ₁ and s ₂ values.

To understand the thermal impact of the obtained s values, we compared the net heat contributed by the sources with the corresponding melt energy flux. Averaging over all the selected records for the CR_i method, we found that the net estimated heating of the layers due to the sources, 6 ± 9 W m⁻², was not substantial relative to the uncertainty, by comparison to the mean conductive flux of 36 ± 6 W m⁻²; similar results were found using both the CR_h and MC_h methods. The same was evident using the data from Khumbu Glacier (Rowan and others, 2021), where the net estimated heating of the layers due to the sources varied from 4 ± 12 to 16 ± 17 W m⁻².

The above findings suggest that on average, the debris layer can be well approximated by a horizontally homogeneous one-dimensional conductor, where the internal heat sources (e.g., those due to the latent heat exchange or the lateral inhomogeneity) play a relatively minor role. It may be worth exploring if this result is a general feature of debris-covered glaciers in the Himalaya or elsewhere.

5.4. Sub-debris ablation rates on Satopanth Glacier

Considering the selected pits with dz ₂ = dz ₁ that were used for the analysis, the glaciological ablation rates varied between 0.2 and 1.3 cm w.e. d⁻¹ (Fig. 6). The debris thickness in these pits ranged from 22 to 77 cm. In general, for any given period, higher ablation rates were found at pits with thinner debris and vice versa, which is consistent with the results from Shah and others (2019) and other studies (e.g., Winter-Billington and others, 2020). The CR_h, CR_i, MC_h and MC_i estimates of ablation rates varied between 0.02–1.9, 0.2–1.9, 0.04–2.3 and 0.2–2.1 cm w.e. d⁻¹, respectively.

Fig. 6. Comparison of ablation rate estimates from the CR_h (a), CR_i (b), MC_h (c), MC_i (d) methods, with that obtained from the observed glaciological method using ablation stakes. Each point is coloured by dz ₂/dz ₁ of the corresponding pit. The asserted RMSE and $R^2_{\rm {adj}}$ were estimated using the selected temperature records with dz ₂/dz ₁ = 1 (see text for details). The solid grey line is a guide to the eye that denotes perfect match.

The comparisons between the observed glaciological ablation rates with those obtained using the debris temperature profiles (Fig. 6) yielded RMSEs of 0.48, 0.32, 0.36 and 0.30 cm w.e. d⁻¹ using the CR_h, CR_i, MC_h and MC_i methods, respectively (Table 4). The number of fitted parameters in these four methods was 2, 3, 2 and 4, respectively; therefore, we also used the adjusted R ² ($R^2_{\rm {adj}}$) as a measure of the goodness-of-fit of the models. The $R^2_{\rm {adj}}$ values were 0.06, 0.30, 0.30 and 0.43, respectively, for the CR_h, CR_i, MC_hand MC_i methods (Table 4). This suggested that the MC_i method most accurately reproduced the observed ablation rates, and the least accurate estimates were found using the CR_h method.

Table 4. Statistical measures used to compare the ablation rate estimates from the four different methods with the glaciological data considering only the selected temperature records (see Section 5.1 for details)

The values of RMSE showed that the models performed adequately, with errors in the order of a few cm, and accumulating to 36–58 cm w.e. over a 4-month ablation season (e.g., Fig. 7). The $R^2_{\rm {adj}}$ values were not as compelling; the value of 0.06 for the CR_h method was particularly low, but even the value of 0.43 calculated for the MC_i method was not as high as might be expected. This is due to the penalties applied in the calculation of $R^2_{\rm {adj}}$ for the predictor variables that do not substantially increase the goodness-of-fit of the model. As stated above, the source term, s, did not contribute substantially to ablation rates; therefore, the $R^2_{\rm {adj}}$ of the regression models were penalised for including s as a predictor variable.

Fig. 7. The observed and inferred ablation rates for pit SBP8 during the ablation season of 2017. The ablation estimates from all other pits (Table 1) are in Supplementary Figure S9.

Figure 7 provides an example of the cumulative ablation predicted for an individual pit. Similar plots for the remaining pits can be found in Supplementary Figure S9. It can be seen that the accuracy of sub-seasonal ablation estimates improved when the inhomogeneous methods (CR_i and MC_i) were applied to individual pits. The MC_i method was in near-perfect agreement with the measurements for the first 3 months.

The mean biases (relative biases) between the CR_h, CR_i, MC_h and MC_i estimates of ablation rates and the observed ablation rates were 0.34 (39%), − 0.03 (− 4%), 0.16 (19%) and − 0.10 (− 11%) cm w.e. d⁻¹, respectively (Table 4). The mean biases were higher in the homogeneous methods (CR_h and MC_h) relative to those in the inhomogeneous methods (CR_i, MC_i) (Figs 6, 7). The CR_i method was the least, and the CR_h method was the most, biased.

Figure 6 shows that ablation rates estimated using data where the thermistors were unevenly spaced through the debris layer were less accurate than those made with data where dz ₁ = dz ₂. This was consistent with the result of synthetic experiment 1 (see Section 5.1), which showed that the accuracy of estimates of κ _eff decreased with increasingly uneven sensor spacing.

The uncertainties associated with ablation rates obtained by the glaciological method were 0.1–0.7  cm w.e. d⁻¹ (Supplementary Fig. S10) (Shah and others, 2019). By comparison, the uncertainties in estimated ablation rates that were made using the MC_i method ranged from 0.04 to 0.72 cm w.e. day⁻¹ (Supplementary Fig. S10). The uncertainties in the estimated ablation rates from the other three methods were comparable to that of the MC_i method, but, as discussed above, these methods were less accurate.

An implication of the above findings is that vertical temperature profile measurements can be used to estimate sub-debris ablation (Fig. 6) and its sub-seasonal variability (e.g., Fig. 7) with accuracies comparable to that of the glaciological method in the thickly debris-covered parts of the glacier. The automated temperature sensors provide continuous data, which can be used to obtain the seasonal to sub-seasonal ablation rates with a fraction of the physical labour that is required to obtain the same information using the glaciological method. Therefore, this method may be particularly suitable for relatively inaccessible debris-covered HKKH glaciers, where in situ ablation monitoring is logistically challenging and, currently sparse. The findings suggest that the assumption of purely conductive, vertical heat transport within the debris layer provides a reasonably accurate description of the sub-debris heat fluxes over the debris-covered ablation zone (Section 5.3). The departures from such an idealised model do not lead to errors in the sub-debris ablation estimates that are important given the level of uncertainty typically present in the corresponding glaciological estimates.

5.4.1. Effect of the experimental set-up on the accuracy of sub-debris ablation estimation

Based on the above results, irregular sensor spacing leads to biased estimates of thermal diffusivity as well as sub-debris ablation rates. In the idealised settings of synthetic experiment 1 (Section 4.3), we found 10% bias in dz ₂/dz ₁ leads to <5% biases in both κ _eff as well as κ ₂ in the MC_i or CR_i methods (Fig. 4a). In reality, those biases could be different due to complexities arising from the inhomogeneities in the debris layer, debris thickness, the elevation of the pits, different surface temperatures and so on. That is why we used a comparison of sub-debris ablation with the glaciological observations to check the influence of different parameters related to the experimental setup on the accuracy of sub-debris ablation estimation (Fig. 8).

Fig. 8. (a) In the MC_i method, we plotted the maximum ablation mismatch percentages (Fig. 6d) and the corresponding mean dz ₂/dz ₁. (b) The mean values of the ratios dz ₂/d, Δ₁/d and Δ₂/Δ₁ were plotted for different maximum ablation mismatch percentages (see Section 5.4 for more details). Each point was coloured by the corresponding mean κ _eff.

We considered the set of temperature records where the MC_i method reproduced the observed glaciological ablation rates within a 20% error. The corresponding RMSE between observed and inferred ablation rate was 0.12 cm w.e. d⁻¹. Denoting the separation between the debris surface to the top sensor and that between the bottom sensor to the ice surface by Δ₁ and Δ₂, respectively, the experimental set up for this set was characterised by, dz ₂/dz ₁ = 1 ± 0.03, 0.14 < dz ₂/d < 0.18, 0.47 < Δ₁/d < 0.56 and 0.29 < Δ₂/Δ₁ < 0.39. Beyond these ranges, both the effective thermal diffusivity and the ablation mismatch increased systematically (Fig. 8).

5.5. Recommendations for experimental design

Based on the above discussion, we recommend the following protocol for accurate estimation of sub-debris ablation using vertical debris-temperature profiles.

• • Place three temperature sensors within the debris layer.

• • Maintain an equal spacing between the successive sensors with a tolerance of 3%.

• • Set the sensor spacing to be $\sim {1\over 5}$ th of the debris thickness at the location.

• • The top sensor should be placed approximately at the middle of the debris layer.

• • Discard the debris temperature data for the first 3 days after the installation to allow thermal transients to disappear.

• • Analyse the debris temperature profiles using a finite-difference method that incorporates the vertical inhomogeneity of the debris layer, preferably the MC_i method introduced here.

6. Summary and conclusions

We have presented an analysis of the accuracy of thermal diffusivity and sub-debris ablation estimates made using in situ measurements of temperature within the supraglacial debris on Satopanth Glacier during two ablation seasons. We have compared the four different methods of analysing the debris temperature profiles. The methods are based on one-dimensional heat conduction through a single-layered or a two-layered conductor. The accuracy of the methods was evaluated using idealised synthetic experiments and by direct comparison with glaciological observations of ablation at Satopanth Glacier. We assessed the effects of the vertical spacing of the temperature sensors within the debris layer and vertically inhomogeneous thermal properties. Our main conclusions are:

• • The most accurate estimates of thermal diffusivity and sub-debris ablation are obtained from data collected at equispaced temperature sensors.

• • Taking into account the inhomogeneity of the debris layer by analysing the data based on a multi-layered model improves the accuracy.