1. Introduction
Glaciers in the Hindu-Kush Karakoram Himalaya (HKKH) constitute an important freshwater reserve for downstream populations (Azam and others, Reference Azam2021). 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, Reference Herreid and Pellicciotti2020). 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, Reference Østrem1959). 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, Reference Banerjee and Shankar2013; Anderson and Anderson, Reference Anderson and Anderson2016; Banerjee, Reference Banerjee2017; Ferguson and Vieli, Reference Ferguson and Vieli2020).
The influence of supraglacial debris on ablation depends on its thickness and thermal properties (Mihalcea and others, Reference Mihalcea2006; Reid and others, Reference Reid, Carenzo, Pellicciotti and Brock2012; Fyffe and others, Reference Fyffe2020). While field observations of debris properties in HKKH are scarce (Conway and Rasmussen, Reference Conway and Rasmussen2000; Nicholson and Benn, Reference Nicholson and Benn2013; Rounce and others, Reference Rounce, Quincey and McKinney2015; Chand and Kayastha, Reference Chand and Kayastha2018; Rowan and others, Reference Rowan2021), they indicate that debris properties vary within and between glaciers at a range of scales, and through time (Mihalcea and others, Reference Mihalcea2008; Nicholson and others, Reference Nicholson, McCarthy, Pritchard and Willis2018; Shah and others, Reference Shah, Banerjee, Nainwal and Shankar2019). 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, Reference Nicholson, McCarthy, Pritchard and Willis2018). This results in considerable uncertainty in many glaciological applications, such as satellite-based debris-thickness estimation (e.g., Mihalcea and others, Reference Mihalcea2008; Foster and others, Reference Foster, Brock, Cutler and Diotri2012; Schauwecker and others, Reference Schauwecker2015) and glacio-hydrological modelling (e.g., Fujita and Sakai, Reference Fujita and Sakai2014; Hagg and others, Reference Hagg2018; Zhang and others, Reference Zhang2019; Steiner and others, Reference Steiner, Kraaijenbrink and Immerzeel2021).
It is generally accepted that the global population of glaciers is undersampled (Mernild and others, Reference Mernild, Lipscomb, Bahr, Radić and Zemp2013) 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, Reference Wester, Mishra, Mukherji and Shrestha2019), yet, direct measurements of surface ablation have been reported for fewer than 20 debris-covered glaciers in this region (Winter-Billington and others, Reference Winter-Billington, Moore and Dadic2020). 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, Reference Wang, Liu, Shangguan, Radić and Zhang2019) 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, Reference Cogley2010). This is a labour-intensive method (Kaser and others, Reference Kaser, Fountain and Jansson2003). From our direct experience (Shah and others, Reference Shah, Banerjee, Nainwal and Shankar2019), 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, Reference Nicholson and Benn2013).
Conway and Rasmussen (Reference Conway and Rasmussen2000) 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, Reference Conway and Rasmussen2000), has been adapted in several studies in HKKH (e.g., Haidong and others, Reference Haidong, Yongjing and Shiyin2006; Nicholson and Benn, Reference Nicholson and Benn2013; Chand and Kayastha, Reference Chand and Kayastha2018; Rowan and others, Reference Rowan2021) and elsewhere (e.g., Nicholson and Benn, Reference Nicholson and Benn2006; Anderson and others, Reference Anderson, Armstrong, Anderson and Buri2021).
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, Reference Shah, Banerjee, Nainwal and Shankar2019) 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 km2, and it spans an elevation range of 3900–6200 m a.s.l (Shah and others, Reference Shah, Banerjee, Nainwal and Shankar2019). 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, Reference Banerjee and Bilal2018). Avalanches contribute a major fraction of the total accumulation of the glacier (Laha and others, Reference Laha2017). 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, Reference Shah, Banerjee, Nainwal and Shankar2019).
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, Reference Valdiya1999). 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, Reference Sakai and Fujita2010).
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, Reference Shah, Banerjee, Nainwal and Shankar2019). 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 (Reference Shah, Banerjee, Nainwal and Shankar2019). 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 (Reference Shah, Banerjee, Nainwal and Shankar2019) 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, Reference Evatt2015) 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, Reference Conway and Rasmussen2000) 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.
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 1, z = 0 and z = dz 2) are indicated with solid black arrows. Figures (a) and (b) show the thickness l 1 (l 2) of the top (bottom) layer as used in the CRi, and MCi 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−3, C = 750 J kg−1 K−1 and ϕ=0.3 are constant in z (Conway and Rasmussen, Reference Conway and Rasmussen2000) and (b) K(z) and s(z) are constant in each layer. Then, for each layer, Eqn (1) can be written as,
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:
Here, the temperature gradient at the debris-ice interface is ${dT\over dz}$
, the latent heat of fusion L f = 3.34 × 105 J kg−1 and the density of water ρ w = 1000 kg m−3.
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, Reference Conway and Rasmussen2000) to a two-layered conductor. We denote the original one-layered model as CRh, and the generalisation to the two-layered model as CRi. 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 MCh and MCi, 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 1 (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, Reference Collier2014; Giese and others, Reference Giese, Boone, Wagnon and Hawley2020), 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 1, z = 0 and z = dz 2 as shown in Fig. 3) are available, the second-derivative term in the right-hand side of Eqn (2) can be approximated as follows,
The lowest order error term is linear in (dz 1 − dz 2) when dz 1 ≠ dz 2, and O(dz 3) otherwise. This suggests that in this finite-difference scheme, the errors are larger in the case of non-uniform sensor spacing, e.g., dz 1 ≠ dz 2.
4.1.1.1. Homogeneous CR method (CRh)
Given the discrete time series of temperature measured at the three different sensors (as shown in Fig. 3), this method (Zhang and Osterkamp, Reference Zhang and Osterkamp1995; Conway and Rasmussen, Reference Conway and Rasmussen2000) employs the following finite-difference approximations to evaluate the derivative terms in Eqn (2) assuming higher-order correction terms (Eqn (4)) to be negligible:
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, Reference Conway and Rasmussen2000). The standard errors of the fitted parameters were taken as the corresponding uncertainties.
4.1.1.2. Inhomogeneous CR method (CRi)
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 CRh method (Eqn (6)) to a two-layered model:
where κ 1 and κ 2 are the thermal diffusivities of the top and bottom layers with thicknesses dz 1 and dz 2, 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 κ 1, κ 2 and s. The uncertainties in these parameters were computed in the same manner as described in the CRh 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
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, Reference Gelman2014). We describe our implementation of this method for the one-layer case (MCh) and the two-layer case (MCi) below.
4.1.2.1. Homogeneous Bayesian method (MCh)
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, Reference Slingerland and Lee2011). The upper boundary condition was fixed using the observed temperature at z = −dz 1 (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 2 were used to compute the sum of squared errors (δ2) relative to the observed temperatures (T obs) as follows.
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( − δ2), 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 mm2 s−1, and −6 × 10−4 to 6 × 10−4 Ks−1, respectively. We searched the two-dimensional space of the parameters (κ, s) using a Monte Carlo procedure in order to minimise δ2(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 (MCi)
This method used two layers, with the boundary being equidistant from the sensors at z = 0 and z = dz 2 (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, κ 1, κ 2, s 1 and s 2, 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 MCh were used for the MCi.
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 (σ 0) were combined with the fit uncertainties (σ 1) 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 2 (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 1 and dz 2. The value of κ was held constant at 1.0 mm2 s−1, and the value of dz 2/dz 1 was varied within the range 1.0–5.7, as reported in the relevant literature (Conway and Rasmussen, Reference Conway and Rasmussen2000; Nicholson and Benn, Reference Nicholson and Benn2006; Brock and others, Reference Brock2010; Reid and others, Reference Reid, Carenzo, Pellicciotti and Brock2012; Nicholson and Benn, Reference Nicholson and Benn2013; Rounce and others, Reference Rounce, Quincey and McKinney2015; Chand and Kayastha, Reference Chand and Kayastha2018; Rowan and others, Reference Rowan2021). In our field measurements, the values of dz 2/dz 1 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 (CRi and MCi), we computed the effective thermal diffusivity (κ eff) using
where l 1 and l 2 are the thickness, and κ 1 and κ 2 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 1 and dz 2, and assessed the trend of κ eff with dz 2/dz 1 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 (Reference Rowan2021). Rowan and others (Reference Rowan2021) 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 (Reference Rowan2021) 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 CRh method to estimate κ for different values of dz 2/dz 1 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 2/dz 1 = 1, s = 0, and used two horizontally homogeneous layers with different κ (κ 1 for the top layer and κ 2 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 κ 1 and κ 2 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, Reference Macfarlane, Hodges and Lux1992; Nicholson and Benn, Reference Nicholson and Benn2013). In both CRh and CRi methods, following Conway and Rasmussen (Reference Conway and Rasmussen2000), the value of ${dT\over dz}$
was obtained by a linear fit to the mean temperature at all three sensors. In the MCh and MCi methods, ${dT\over dz}$
was obtained by a linear fit to the mean simulated temperature at each of the grid points between z = dz 2 and the debris-ice interface (Fig. 3). In the CRi and MCi methods, the thermal diffusivity estimated for the bottom-most layer, κ 2, 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 (Reference Shah, Banerjee, Nainwal and Shankar2019), 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 (Reference Shah, Banerjee, Nainwal and Shankar2019).
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 2/dz 1 = 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 2/dz 1 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 2/dz 1 > 1 (Figs 5a–5d). The trend of increasing κ eff with increasing dz 2/dz 1 > 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 2/dz 1 > 1 was also observed in the debris temperature data from Khumbu Glacier, which were analysed using the CRh 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 2/dz 1 (see the text for details). The grey horizontal line denotes the true value of κ eff (1 mm2 s−1) used in the synthetic experiments. The vertical grey band is the range of dz 2/dz 1 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 κ 2/κ 1 used in the forward model. Here we kept dz 2/dz 1 = 1.
Fig. 5. (a)–(d) Estimated values of κ eff from four methods are plotted as a function of dz 2/dz 1 for all the records. The colours of the symbols denote the month of the temperature measurement. (e) The CRh method estimates of κ eff plotted for different dz 2/dz 1 in a pit (KH1) from Khumbu Glacier during ablation season of 2014 (Rowan and others, Reference Rowan2021).
While the above dependence of the biases in inferred κ eff on dz 2/dz 1 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, Reference Anderson and Anderson2010), 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 κ 2/κ 1, 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 CRh, CRi, MCh and MCi methods were 0.62, 0.08, 0.07 and 0.03 mm2 s−1, respectively. As these values show, the accuracy of the CRh method was substantially lower than that of the other methods. The performance of the CRi, MCh and MCi methods was comparable, but ultimately, the MCi 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 1 = dz 2. 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 CRh and CRi 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 (MCh and MCi), and 30 records with respect to the finite-difference methods (CRh and CRi).
In the MCi method, the mean, standard deviation and range of κ 1 (κ 2) were 1.2 (2.0), 0.6 (1.0), 0.3–0.32 (0.7–4.3) mm2 s−1, 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, Reference Conway and Rasmussen2000; Nicholson and Benn, Reference Nicholson and Benn2013; Rounce and others, Reference Rounce, Quincey and McKinney2015; Chand and Kayastha, Reference Chand and Kayastha2018). 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 (CRi and MCi) 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 2 = dz 1, 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 CRi and MCi methods, and that from the literature were 1.3 (0.8), 1.5 (0.9) and 0.9 (0.5) mm2 s−1, 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, Reference Mihalcea2008; Foster and others, Reference Foster, Brock, Cutler and Diotri2012; Schauwecker and others, Reference Schauwecker2015). 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, Reference Nicholson and Benn2006; Reid and Brock, Reference Reid and Brock2010; Lejeune and others, Reference Lejeune, Bertrand, Wagnon and Morin2013) or glacio-hydrological models (Fujita and Sakai, Reference Fujita and Sakai2014; Hagg and others, Reference Hagg2018; Zhang and others, Reference Zhang2019), 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, Reference Rowan2021) (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 1 and s 2 obtained using the MCi 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 δ2 of the fits were not sensitive to s 1 and s 2, while it had a sharp minima in the κ 1 and κ 2 plane. We also found that the MCi method is not sensitive enough to estimate vertical inhomogeneities in s. This was indicated by an anti-correlation between the observed values of s 1 and s 2, so that s 1 + s 2 ≈ 0 at almost all the points. This implied that a similar δ 2 can be obtained by simply exchanging the sign of s 1 and s 2 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 CRi method, we found that the net estimated heating of the layers due to the sources, 6 ± 9 W m−2, was not substantial relative to the uncertainty, by comparison to the mean conductive flux of 36 ± 6 W m−2; similar results were found using both the CRh and MCh methods. The same was evident using the data from Khumbu Glacier (Rowan and others, Reference Rowan2021), where the net estimated heating of the layers due to the sources varied from 4 ± 12 to 16 ± 17 W m−2.
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 2 = dz 1 that were used for the analysis, the glaciological ablation rates varied between 0.2 and 1.3 cm w.e. d−1 (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 (Reference Shah, Banerjee, Nainwal and Shankar2019) and other studies (e.g., Winter-Billington and others, Reference Winter-Billington, Moore and Dadic2020). The CRh, CRi, MCh and MCi 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−1, respectively.
Fig. 6. Comparison of ablation rate estimates from the CRh (a), CRi (b), MCh (c), MCi (d) methods, with that obtained from the observed glaciological method using ablation stakes. Each point is coloured by dz 2/dz 1 of the corresponding pit. The asserted RMSE and $R^2_{\rm {adj}}$
were estimated using the selected temperature records with dz 2/dz 1 = 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−1 using the CRh, CRi, MCh and MCi 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 2 ($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 CRh, CRi, MChand MCi methods (Table 4). This suggested that the MCi method most accurately reproduced the observed ablation rates, and the least accurate estimates were found using the CRh 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 CRh method was particularly low, but even the value of 0.43 calculated for the MCi 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 (CRi and MCi) were applied to individual pits. The MCi method was in near-perfect agreement with the measurements for the first 3 months.
The mean biases (relative biases) between the CRh, CRi, MCh and MCi 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−1, respectively (Table 4). The mean biases were higher in the homogeneous methods (CRh and MCh) relative to those in the inhomogeneous methods (CRi, MCi) (Figs 6, 7). The CRi method was the least, and the CRh 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 1 = dz 2. 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−1 (Supplementary Fig. S10) (Shah and others, Reference Shah, Banerjee, Nainwal and Shankar2019). By comparison, the uncertainties in estimated ablation rates that were made using the MCi method ranged from 0.04 to 0.72 cm w.e. day−1 (Supplementary Fig. S10). The uncertainties in the estimated ablation rates from the other three methods were comparable to that of the MCi 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 2/dz 1 leads to <5% biases in both κ eff as well as κ 2 in the MCi or CRi 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 MCi method, we plotted the maximum ablation mismatch percentages (Fig. 6d) and the corresponding mean dz 2/dz 1. (b) The mean values of the ratios dz 2/d, Δ1/d and Δ2/Δ1 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 MCi 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−1. Denoting the separation between the debris surface to the top sensor and that between the bottom sensor to the ice surface by Δ1 and Δ2, respectively, the experimental set up for this set was characterised by, dz 2/dz 1 = 1 ± 0.03, 0.14 < dz 2/d < 0.18, 0.47 < Δ1/d < 0.56 and 0.29 < Δ2/Δ1 < 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 MCi 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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2022.35.
Data availability
All the observed debris-temperature data are available using this link : https://osf.io/kysg6/.
Acknowledgements
We acknowledge help from Gajendra Badwal, Surendra Badwal, Sunil Singh Shah, Nepalese porters, and people from Mana and Badrinath villages during the fieldwork. S.L. acknowledges financial support from the DST-INSPIRE fellowship (IF170526). The field campaigns on Satopanth Glacier were supported by The Institute of Mathematical Sciences, Chennai. A.W.-B. and M.K. were supported by a Canada Foundation for Innovation grant, a MITACS Globalink Research Award, and a Canadian National Science and Engineering Council Discovery Grant to M.K. We thank the two anonymous reviewers and the Associate Chief Editor Nicolas Cullen and Scientific Editor Dan Shugar, for providing insightful comments and suggestions.
Author contributions
A.B., S.L., M.K., A.W.-B., R.S. and H.C.N. designed the field experiments. A.B., S.L., R.S. and A.W.-B. did the field experiments. A.B. and S.L. designed the theoretical analysis with inputs from R.S. and A.W.-B. S.L., A.B., R.S. and A.W.-B. wrote the paper with help from H.C.N. and M.K.
Open access
