A Comparative Study of Methods for Estimating the Thickness of Glacial Debris: A Case Study of the Koxkar Glacier in the Tian Shan Mountains

Abstract

The local or overall mass balance of a glacier is significantly influenced by the spatial heterogeneity of its overlying debris thickness. Accurately estimating the debris thickness of glaciers is essential for understanding their hydrological processes and the impact of climate change. This study focuses on the Koxkar Glacier in the Tian Shan Mountains, using debris thickness data to compare the accuracy of three commonly used approaches for estimating the spatial distribution of debris thickness. The three measurement approaches include two empirical relationships between the land surface temperature (LST) and debris thickness approaches, empirical relationship approach 1 and empirical relationship approach 2, and the energy balance of debris approach. The analysis also explores the potential influence of topographic factors on the debris distribution. By incorporating temperature data from the debris profiles, this study examines the applicability of each approach and identifies areas for possible improvement. The results indicate that (1) all three debris thickness estimation approaches effectively capture the distribution characteristics of glacial debris, although empirical relationship approach 2 outperforms the others in describing the spatial patterns; (2) the accuracy of each approach varies depending on the debris thickness, with the energy balance of debris approach being most accurate for debris less than 50 cm thick, while empirical relationship approach 1 performs better for debris thicker than 50 cm and empirical relationship approach 2 demonstrates the highest overall accuracy; and (3) topographic factors, particularly the elevation, significantly influence the accuracy of debris thickness estimates. Furthermore, the empirical relationships between the LST and debris thickness require field data and focus solely on the surface temperature, neglecting other influencing factors. The energy balance of debris approach is constrained by its linear assumption of the temperature profile, which is only valid within a specific range of debris thickness; beyond this range, it significantly underestimates the values. These findings provide evidence-based support for improving remote-sensing methods for debris thickness estimation.

1. Introduction

Debris refers to the rocky material accumulated on the surface of a glacier, and a glacier where the majority of the ablation area is covered by continuous surface debris is classified as a debris-covered glacier [1]. Debris-covered glaciers are widely distributed across regions such as the Tian Shan, the Himalayas, the Alps, and the Andes, comprising approximately 7.3% of the global glacier-covered area [2,3,4]. The presence of debris alters the energy exchange between the glacier and the atmosphere, significantly impacting glacier ablation. Studies have shown that debris, with its low albedo, significantly absorbs solar radiation. Consequently, a thin layer of debris accelerates the ablation of the underlying ice. However, once the debris exceeds a certain thickness, it acts as an insulating layer, reducing the glacier melt due to its thermal resistance [5,6,7,8]. Therefore, obtaining precise data on the spatial distribution of debris is crucial for understanding the ablation processes and mass balance of debris-covered glaciers [9,10]. Currently, due to the limited data on regional debris thickness, most regional or global glacier models struggle to accurately assess the impact of debris. Some models use a simplified approach by utilizing different ablation factors for bare ice and debris-covered areas, which may significantly reduce the accuracy of glacier ablation estimates [11,12,13]. Therefore, precise measurements of debris thickness are essential for improving predictions of glacier changes and assessing changes in glacier runoff and its impacts [14,15,16].

Field measurement is the most accurate approach for obtaining debris thickness data [17,18,19]. However, this approach is labor-intensive, time-consuming, and particularly challenging in remote or inaccessible glacier regions, limiting its effectiveness in understanding the broader distribution of glacier debris thickness. In response, scholars have increasingly adopted remote-sensing technology, developing various debris thickness estimation methods that leverage its capabilities, such as large-scale monitoring, high spatial and temporal resolution, and cost-effectiveness. These features allow researchers to gather regional or global debris thickness data with relative ease.

Remote-sensing technologies for debris thickness estimation are generally classified into three main categories: the synthetic aperture radar (SAR) surface debris detection approach, the empirical land surface temperature (LST) debris thickness relationship approach, and the remote-sensing and energy balance approach, used in combination with estimation models. The SAR debris detection approach uses long-wavelength radar, which possesses limited penetrative ability, to estimate debris thickness. For example, Huang et al. [20] utilized fully polarized L-band radar, observing that the body-scattering component increases with thicker debris. However, due to the limited penetration depth of SAR, this method is not suitable for estimating thicker debris layers. Recently, the empirical LST and thickness approach and energy balance approach have emerged as the dominant approaches for estimating debris thickness. The strength of the empirical LST and thickness approach lies in its ability to establish linear or nonlinear statistical relationships between the LST and debris thickness [21,22,23]. For example, Mihalcea et al. [24] used ASTER imagery to invert the LST in the study area, correlated the measured debris thickness with the LST, and estimated the debris thickness distribution on the Baltoro Glacier. The energy balance approach, which combines remote-sensing and energy balance analysis, offers a robust physical mechanism. Its enhanced capabilities stem from integrating LST data obtained via remote-sensing with meteorological data to simulate the energy balance process of the debris layer and assess the debris thickness [17,25,26,27]. Stewar et al. [28] and Zhang et al. [29] improved this method by introducing the thermal resistance coefficient, which refines the calculation of the relationship between heat flux and debris thickness, thus improving the estimation efficiency.

Although numerous approaches for estimating debris thickness have been developed, there remains a lack of systematic comparative assessments of their applicability to debris-covered glaciers, primarily due to the scarcity of measured data. In this study, we focus on the Koxkar Glacier on the south slope of Tomur Peak as the study area. Utilizing 556 high-density debris thickness measurements obtained from 2008 to 2009, we evaluate the accuracy and analyze the errors of widely used estimation approaches, exploring potential directions for improving these methods.

2. Study Area and Data

2.1. Study Area

The Tian Shan Tomur Peak region is one of the key areas in Central Asia known for its high concentration of surface debris-covered glaciers. The Koxkar Glacier, a typical debris-covered glacier, is situated on the southern slope of the region [30]. According to the Second Glacier Inventory of China, the average elevation of the Koxkar Glacier is 4348.4 m, with a maximum elevation of 5652.7 m and a minimum elevation of 3000.8 m. The glacier covers an area of 86.88 ± 1.82 km², with a volume ranging from 13.8 to 16.97 km³ [31,32]. The debris-covered area extends across 20.62 ± 1.63 km², primarily below 3800 m, accounting for 23.73% of the overall glacier area and approximately 83% of the total ablation area. The ice tongue is largely covered by substantial, relatively thick debris, with the thickness increasing progressively from the upper part of the ablation zone downwards, reaching over 2 m at the glacier terminus [33]. Due to the debris cover, the melting characteristics of the Koxkar Glacier are primarily manifested as a decrease in thickness, which differs from the melting characteristics of most clean glaciers that exhibit a sharp reduction in area at the termini [34,35]. The most intense ablation occurs in areas with debris less than 10 cm thick and at elevations of approximately 3800 m [34].

2.2. Data Sources

This study utilized 556 debris thickness records from field measurements conducted between 2008 and 2009 through manual excavation, with the measurement points evenly distributed across the debris surface. The data were processed using the ordinary Kriging method to generate a debris thickness distribution map with a 30 m spatial resolution. For the inversion of the LST and related meteorological elements, Landsat 5 TM data from 30 September 2009 (13:00), a cloudless and clear day, were selected.

To simulate the debris surface energy balance, data were collected during a satellite overpass around the Koxkar Glacier area from six automatic weather stations (Figure 1). The automatic weather stations were positioned at elevations of 2900, 3200, 3400, 3700, and 4200 m, recording meteorological elements such as the air temperature, wind speed, anemometer height, relative humidity, and shortwave radiation. Data from two temperature profiling stations were used to observe the temperature variations and characteristics within the debris. Due to damage at the site2 observatory in 2024, field observations from the two temperature profiling stations were conducted over different periods: site1 collected data from 7 to 14 August 2024, while site2 collected data from 3 to 10 August 2023. Measurements were recorded every 10 min. Moreover, at site1, probes were placed at heights of 0, 10, 20, 30, 40, 50, and 60 cm, while at site2, they were placed at heights of 0, 40, 70, and 100 cm.

3. Methods

3.1. Inversion of the LST

The approach used to invert the LST using the radiative transfer equation is widely applicable to thermal infrared data from TM sensors. The LST ($T_S$) is calculated as follows [36,37]:

$$T_S={\displaystyle \frac{T_{sen}}{1+\left(\lambda \frac{T_{sen}}{\rho }\right)\ln \epsilon }}+273.15$$

(1)

$$\rho ={\displaystyle \frac{hc}{k_B}}$$

(2)

where $\lambda$ is the wavelength of the thermal infrared band, and $\epsilon$ is the surface emissivity that is determined with reference to the approach utilized by Chang et al. [38] and applies to the cold region. In addition, $h$ is Planck’s constant, $c$ is the speed of light, and $k_B$ is the Boltzmann constant. The brightness temperature of the sensor (K) is calculated as follows:

$$T_{sen}={\displaystyle \frac{K_2}{\ln \left(1+\frac{K_1}{L_{sen}}\right)}}$$

(3)

$$L_{sen}=GAIN\times DN+Bias$$

(4)

where $K_1$ and $K_2$ are pre-launch calibration constants, and $L_{sen}$ is the sensor radiation intensity. For the Landsat 5 TM band, $K_1$ = 607.76 W·m⁻²·sr⁻¹·μm⁻¹ and $K_2$ = 1260.56 K. Additionally, $DN$ is the calibration coefficient gain and $Bias$ is the offset.

3.2. Approaches for Estimating Debris Thickness

3.2.1. Empirical LST and Debris Thickness Approach

As a poor thermal conductor, debris exhibits the characteristic of ineffective heat transfer. Beneath it, the underlying glacial ice exerts a significant cooling effect. This physical property leads to varying influences on the heat conduction and the cooling effect of glacial ice, depending on the thickness of the debris layer. Specifically, thicker debris is better at insulating and mitigating the cooling effects of the underlying glacial ice, whereas thinner debris layers are more sensitive to these cooling effects. Consequently, debris of different thicknesses displays distinct temperature differentials.

The empirical LST and thickness approach is based on the assumption that the debris thickness is the primary factor influencing the LST [21,24,39]. This allows for the use of the LST to establish a correlation with debris thickness. The approach is straightforward and does not require meteorological data, making it simple to implement and easily applicable to a wide range of debris-covered glaciers. Kraaijenbrink et al. [11] introduced the following empirical equation (empirical relationship approach 1):

$$h_d=e^{{\displaystyle \frac{\left(T_s-T_{min}\right)\ln h_{max}}{T_{p95}-T_{min}}}}$$

(5)

where $h_d$ is the surface debris thickness (cm), $T_{min}$ is the minimum LST, $T_{p95}$ is the LST of the 95th quartile of the debris-covered area, and $h_{max}$ corresponds to the LST of the 95th quartile corresponding to the maximum value of the debris thickness. Based on the LST inversion results and the actual observations from the Koxkar Glacier, $T_{min}$ = 0.12 °C, $T_{p95}$ = 21.71 °C, and $h_{max}$ = 275.32 cm.

Rounce et al. [40] introduced a nonlinear equation (empirical relationship approach 2) based on a debris-covered glacier energy balance model. By deriving the debris LST ($T_s$) under various debris thicknesses ($h_d$), they simulated the relationship between the two as follows:

$$h_d={\left({\displaystyle \frac{T_s{b}^c}{a-T_s}}\right)}^{{\displaystyle \frac{1}{c}}}$$

(6)

where a, b, and c are empirical parameters used to capture the nonlinear relationship between the debris thickness and LST. To optimize the application of this model to the Koxkar Glacier, we fitted 556 measured debris thickness data with Equation (6) to systematically estimate the key parameters in the formula. After optimization, the final precise values for the parameters obtained are a = 33.28, b = 50.80, and c = 0.64.

3.2.2. Energy-Balance of Debris Layer Approach

The energy balance of the debris layer approach has recently become a popular method for estimating debris thickness. This model leverages remote-sensing inversion LST data and meteorological observations to solve the energy balance equation for the debris layer, allowing for debris thickness estimation. According to the energy balance equation, the conduction heat flux within the debris layer ($Q_D$) is calculated as follows [41]:

$$Q_D=R_n+H+LE$$

(7)

where $R_n$ is the net radiation flux, $H$ is the sensible heat flux, and $LE$ is the latent heat flux, as measured in units of W·m⁻².

The heat flux calculation for the debris is based on the following three assumptions: (1) the temperature profile within the debris is linear, (2) the heat conduction within the ice under debris cover is negligible, and (3) the temperature at the interface between the debris and the ice is 0 °C [14,27,42] (F). Under these assumptions, the heat flux at the debris layer ($Q_D$) is expressed as follows:

$$Q_D=K{\displaystyle \frac{T_S}{h_d}}$$

(8)

where $K$ is the effective thermal conductivity of the debris (0.96 W·m⁻¹·K⁻¹) [28].

Combining the above equations, the glacial debris thickness ($h_d$) can be calculated as follows:

$$h_d=K{\displaystyle \frac{T_S}{R_n+H+LE}}$$

(9)

The net radiative flux and longwave radiation are calculated as follows [43,44]:

$$R_n=(1-\alpha )R_S^{\downarrow }+R_L^{\downarrow }+R_L^{\uparrow }$$

(10)

$$R_L^{\downarrow }=\sigma {\epsilon }_a{\left(T_a+273.15\right)}^4$$

(11)

$$R_L^{\uparrow }=\sigma {\epsilon }_s{\left(T_s+273.15\right)}^4$$

(12)

where $R_S^{\downarrow }$ is the downward shortwave radiation, which is obtained from the meteorological station, $R_L^{\downarrow }$ is the downward longwave radiation, and $R_L^{\uparrow }$ is the upward longwave radiation, all measured in W·m⁻². The albedo of the debris layer ($\alpha$) is set at 0.2, and $\sigma$ represents Stefan-Boltzmann’s constant. Additionally, $T_a$ denotes the air temperature, while ${\epsilon }_a$ and ${\epsilon }_s$ are the air specific emissivity and surface-specific emissivity, respectively.

The model assumes that the debris is dry and that the latent heat flux contributes only a small percentage [29,44] to the overall energy balance, as the satellite imagery was acquired under clear sky conditions. Therefore, it can be neglected during the simulation process. The sensible heat flux [45,46] is calculated using the bulk aerodynamic method:

$$H={\rho }_a{\displaystyle \frac{C_p{k}^{2}u\left(T_a-T_s\right)}{\left(ln\frac{z}{z_{0m}}\right)\left(ln\frac{z}{z_{0t}}\right)}}{({\mathsf{\Phi }}_m{\mathsf{\Phi }}_h)}^{-1}$$

(13)

where ${\rho }_a$ is the air density (0.86 kg·m⁻³) [47], $C_p$ is the specific heat capacity of air, $z_{0m}$ is the aerodynamic roughness length (0.098 m) [48], and $k$ is Karman’s constant. Additionally, $z$ represents the observation altitude, and $z_{0t}$ is the thermodynamic roughness length based on the assumption that $z_{0m}$ = $z_{0t}$. The non-dimensional stability functions for the momentum (${\mathsf{\Phi }}_m$), and heat (${\mathsf{\Phi }}_h$) are expressed as functions of ${Ri}_b$ [49]:

$${({\mathsf{\Phi }}_m{\mathsf{\Phi }}_h)}^{-1}={\left(1-5R{i}_b\right)}^2$$

(14)

The bulk Richardson number (${Ri}_b$) is expressed as follows:

$$R{i}_b={\displaystyle \frac{g\left(T_a-T_s\right)\left(z-z_{0m}\right)}{\left(\frac{T_a+T_s}{2}+273.15\right)u^2}}$$

(15)

where $u$ is the wind speed at the observation altitude and $g$ is the gravitational acceleration. The parameters, units and values used in this study are detailed in

Table 1. Abbreviations, units and values of the parameters used in this study.

Parameter	Value (Unit)	Parameter	Value (Unit)
Land surface temperature ($T_S$)	°C	Debris thickness ($h_d$)	cm
Brightness temperature ($T_{sen}$)	K	Effective thermal conductivity ($K$)	0.96 W·m−1·K−1
Wavelength ($\lambda$)	nm	Downward shortwave radiation ($R_S^{\downarrow }$)	W·m−2
Surface emissivity ($\epsilon$)	-	Downward longwave radiation ($R_L^{\downarrow }$)	W·m−2
Planck’s constant ($h$)	6.626 × 10−34 J·s	Upward longwave radiation ($R_L^{\uparrow }$)	W·m−2
Speed of light ($c$)	3 × 108 m·s−1	Albedo ($\alpha$)	0.2
Boltzmann’s constant ($k_B$)	1.38 × 10−23 J·K−1	Stefan-Boltzmann’s constant ($\sigma$)	5.67 × 10−8 W·m−2·K−4
Pre-launch calibration constants ($K_1$)	608 W·m−2·sr−1·μm−1	Air temperature ($T_a$)	°C
Pre-launch calibration constants ($K_2$)	1260.56 K	Air specific emissivity (${\epsilon }_a$)	-
Calibration coefficient gain ($DN$)	0.055158	Surface-specific emissivity (${\epsilon }_s$)	-
Offset ($Bias$)	1.2378	Air density (${\rho }_a$)	0.86 kg·m−3
Minimum LST ($T_{min}$)	°C	Specific heat capacity of air ($C_p$)	1004 J−1·KG−1·K−1
LST of the 95th quartile ($T_{p95}$)	°C	Aerodynamic roughness length ($z_{0m}$)	0.098 m
Maximum of the debris thickness ($h_{max}$)	cm	Karman’s constant ($k$)	0.4
Conduction heat flux within debris ($Q_D$)	W·m−2	Wind speed ($u$)	m·s−1
Net radiation flux ($R_n$)	W·m−2	Observation altitude of wind speed ($z$)	2 m
Sensible heat flux ($H$)	W·m−2	Thermodynamic roughness length ($z_{0t}$)	0.098 m
Latent heat flux ($LE$)	W·m−2	Gravitational acceleration ($g$)	9.8 m·s−2

.

3.3. Evaluation Indicators

In this study, we compare and validate the results of each estimation approach for debris thickness using field-measured data. The accuracy assessment metrics employed include the following. (1) Mean error (ME): This metric reflects the average difference between the estimation approaches, highlighting both high and low estimation scenarios. (2) Mean absolute error (MAE): By disregarding the sign of the error values, the MAE indicates the reliability of the approaches. (3) Root mean square error (RMSE): This metric measures the deviation between the estimated values and the true values, with heightened sensitivity to larger deviations. (4) Median absolute error (MedAE): Less influenced by outliers, the MedAE is suitable for assessing the overall accuracy of the values, particularly in areas with significant spatial heterogeneity. The formulas for these metrics are as follows:

$$ME={\displaystyle \frac{1}{n}}{{\displaystyle \sum }}_{j=1}^n\left(h_i-h_{io}\right)$$

(16)

$$MAE={\displaystyle \frac{1}{n}}{{\displaystyle \sum }}_{j=1}^n\left|h_i-h_{io}\right|$$

(17)

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}{{\displaystyle \sum }}_{j=1}^n{\left(h_i-h_{io}\right)}^2}$$

(18)

$$MedAE=median\left|h_i-h_{io}\right|$$

(19)

where $h_i$ is the value of the height of the debris to be assessed at the $i$th location; and $h_{io}$ is the value of the true thickness at the same location.

To evaluate the overall performance of each estimation approach, we employed the composite rating indicator ($M_R$) [50]. This indicator facilitates a comprehensive ranking of the estimation approaches using the following formula:

$$M_R=1-{\displaystyle \frac{1}{1\times n\times m}}{{\displaystyle \sum }}_{i=1}^n{r}_i$$

(20)

where $n$ is the number of indicators used for the accuracy evaluation, $m$ is the number of estimation methods involved in the assessment, and $r$ is the ranking of each approach for a particular accuracy evaluation indicator. A higher $M_R$ value, closer to 1, indicates the stronger capability of the estimation methods.

4. Results and Analysis

4.1. Results of Debris Thickness Estimation

The spatial distribution of debris on the Koxkar Glacier, as demonstrated in Figure 2, shows an average debris thickness of 0.5 m, with a median value of 0.4 m, and a range of extending from 0.01 to 3.20 m. The spatial distribution of supraglacial debris generally demonstrates a thinning trend as the glacier elevation increases. This inverse relationship between debris thickness and elevation can be attributed to various geomorphological and climatic factors. At the edges of the glacier ablation zone, the complex and unstable terrain serves as an additional source of debris material, facilitating the gradual accumulation of supraglacial debris. As a result, the debris thickness progressively increases from the interior regions toward the outermost boundaries of the glacier.

Figure 3 presents the spatial distribution of debris on the Koxkar Glacier as estimated by the three estimation approaches, showing consistency with the field measurements. The results show the characteristic pattern of debris distribution, where the thickness increases as the distance to the glacier terminus decreases, and the debris thickness gradually increases from the glacier’s interior toward the exterior edges. However, the range of estimated values produced by each approach varies significantly. The thickness measured using empirical relationship approach 1 ranged from 0.01 to 3.54 m, with a 95% quartile thickness of 1.71 m and a mean of 0.55 m. In contrast, empirical relationship approach 2 estimated a maximum thickness of 1.63 m, with a 95% quartile thickness of 0.85 m and a mean thickness of 0.37 m. The energy balance of debris approach results showed a thickness range of 0.01 to 1.18 m, with a 95% quartile thickness of 0.65 m. The debris thickness calculated by the latter two approaches consistently falls below 1 m, while the field measurements indicate the debris thickness on the Koxkar Glacier can exceed 3 m, and at an altitude of 3500 m or lower, the debris thickness on both sides reaches over 1 m. This discrepancy suggests that the accuracy of empirical relationship approach 2 and the energy balance of debris approach is limited, particularly in evaluating areas with significant debris accumulation.

4.2. Spatial Distribution of Errors in Estimation Approaches

Figure 4 illustrates the accuracy of the three approaches in estimating the debris thickness, as well as the spatial distribution of the high and low estimates. The absolute difference can be used as a standard for determining the accuracy. The results show that all three estimation approaches display a pattern of concentrated areas with both overestimated and underestimated debris thickness. For empirical relationship approach 1, the proportion of underestimation is 38.34%, with these areas primarily located along the western edge at elevations above 3500 m and at the terminus of the ice tongue. Regions with an absolute difference of less than 10 cm account for 36.84% of the debris-covered area, while regions with an absolute difference between 10 and 50 cm cover 53.36% of the total debris area. For empirical relationship approach 2, the proportion of overestimation is 40.42%, with the overestimated areas concentrated similarly to those of empirical relationship approach 1, primarily on the eastern side of the ablation area at elevations above 3500 m. However, below this elevation, empirical relationship approach 2 shows higher levels and a greater extent of underestimation compared to empirical relationship approach 1. In terms of the accuracy, 42.37% of the total study area had an absolute difference of less than 10 cm, while 47.92% of the debris-covered area had an absolute difference between 10 and 50 cm. The EDBA exhibited the highest rate of underestimation, accounting for 82.49% of the entire ablation area. Only a small portion of the area above 3500 m showed slight overestimation. For the energy balance of debris approach, areas with an absolute difference of less than 10 cm accounted for 41.32% of all the validation sites, while those with differences between 10 and 50 cm made up 48.04%.

In summary, among the three approaches, empirical relationship approach 2 demonstrates the lowest percentage of absolute differences greater than 50 cm, making it the most accurate method for analyzing the spatial distribution of debris. Regionally, both the empirical LST and thickness approach methods significantly overestimate the debris thickness in the upper regions of the glacier’s surface, while they tend to underestimate the debris thickness at lower elevations. Conversely, the energy balance of debris approach shows almost no significant overestimation but exhibits considerable underestimation in lower-elevation debris areas. For all three approaches, the proportion of areas with absolute differences greater than 50 cm is relatively low, indicating that each approach can effectively assess the debris thickness and distribution on the Koxkar Glacier.

4.3. Effect of Debris Thickness on Estimation Accuracy

To evaluate the accuracy and applicability of the different approaches across different debris thickness ranges, this study classified the ablation zone debris into four thickness layers: thin (0 to 10 cm), medium–thin (10 to 50 cm), medium–thick (50 to 100 cm), and thick (greater than 100 cm) based on the measured thickness. Figure 5 presents the error statistics for each debris thickness layer when different estimation approaches are used. The analysis indicates that empirical relationship approach 1 tends to overestimate the debris thickness across the glacier, whereas empirical relationship approach 2 and the energy balance of debris approach exhibit significant underestimation. All three approaches show considerable underestimation when evaluating thick debris layers, and empirical relationship approach 2 and the energy balance of debris approach underestimate much more than empirical relationship approach 1 (Figure 5a). As the thickness of the debris layers increases, the accuracy of the estimation methods declines significantly (Figure 5b–d), underscoring the clear inverse relationship between the debris thickness and the precision of estimations. This observation emphasizes that achieving accurate estimations for thick debris layers continues to present a substantial challenge, particularly for glaciers characterized by extensive debris accumulation.

The composite ratings indicator (

Table 2. Scores and rankings of the different debris thickness conditions under the three methods based on the composite rating indicator.

Approach	Composite Rating Indicator (${\mathit{M}}_{\mathit{R}}$)/Rank
	Total	≤10 cm	10–50 cm	50–100 cm	>100 cm
Empirical relationship 1	0.42/2	0.17/3	0.00/3	0.50/1	0.67/1
Empirical relationship 2	0.50/1	0.50/2	0.42/2	0.42/2	0.25/2
Energy balance of debris	0.08/3	0.58/1	0.58/1	0.08/3	0.08/3

) revealed significant differences in the accuracy of the three estimation approaches under different debris thickness conditions. Overall, empirical relationship approach 2 had the highest composite rating indicator score (M_R = 0.5), ranking first in terms of the accuracy, while the energy balance of debris approach scored the lowest, indicating the least accuracy for glacier-wide debris thickness estimation. For thin and medium–thin layers, the energy balance of debris approach demonstrated the highest accuracy (M_R = 0.58), followed by empirical relationship approach 2, with empirical relationship approach 1 having the lowest accuracy. When considering medium–thick and thick layers, however, empirical relationship approach 1 provided the most accurate results (M_R = 0.50 and M_R = 0.67), empirical relationship approach 2 came second, and the energy balance of debris approach exhibited significantly lower accuracy compared to the other two approaches.

4.4. Influence of Topographic Factors on Estimation Accuracy

Topographic factors such as the elevation, slope, and aspect interact with the distribution of supraglacial debris thickness, forming a relatively complex relationship [51]. To further explore the influence of each topographic factor of the glacier on each estimation approach, the debris-covered area of the Koxkar Glacier was divided into regions based on the elevation, distance from the glacier flow line, slope, and aspect, respectively. Figure 6 shows the results of the correlation analysis between each factor and the RMSE of the estimation approaches.

The accuracies of all three estimation approaches improve with increasing elevation, exhibiting a strong positive correlation. A distinct threshold marks a turning point in the performance of each method: empirical relationship approach 1 demonstrates superior accuracy at elevations below 3300 m, empirical relationship approach 2 performs best at elevations between 3300 and 3800 m, while the energy balance of debris approach achieves the highest accuracy above 3800 m (Figure 6a).

Figure 6b shows that empirical relationship approach 1 demonstrates significantly higher accuracy compared to the other two approaches, with a strong correlation between its estimation accuracy and the proximity to the mid-stream line. As the distance from the mid-stream line increases, the accuracy of empirical relationship approach 1 decreases noticeably. In contrast, the estimation accuracy of empirical relationship approach 2 and the energy balance of debris approach shows a weaker relationship with the distance from the mid-stream line, displaying no significant variation under the different conditions.

Regarding the slope (Figure 6c), the estimation accuracy of all three approaches declines slightly as the slope increases, likely due to the heavier debris accumulation in areas with steeper gradients. Given the previous discussion on the impact of debris thickness on estimation accuracy, this correlation is to be expected.

As shown in Figure 6d, the aspect influences the estimation accuracy, with variability observed across the different approaches in the same direction. Empirical relationship approach 1 shows higher accuracy in the northeast and north directions, while the energy balance of debris approach performs better in the southeast. However, the maximum RMSE difference between the approaches does not exceed 10 cm, indicating that the differences between the approaches are not substantial.

In summary, glacial features such as the elevation, distance from the glacier centerline, slope, and aspect affect the debris thickness estimation accuracy. Of these factors, the elevation has the most significant impact on the accuracy, making it a key constraint in debris thickness estimation. The accuracy of all three approaches increases significantly with the elevation, suggesting that the elevation and other related topographic factors could serve as valuable indicators for refining debris thickness estimates using remote-sensing methods. This offers a potential pathway for improving the accuracy of debris thickness estimation.

5. Discussion

5.1. Analysis of Empirical Relationship Approach 1 Adaptation

The application of empirical relationship approach 1 to estimate the debris thickness of the Koxkar Glacier confirms the approach’s effectiveness in accurately measuring the debris distribution, particularly due to its simplicity and low data requirements. Empirical relationship approach 1 demonstrates significantly higher accuracy in estimating thick debris (>100 cm) compared to the other approaches. This arises from its use of preset boundary conditions for glacier thickness, enabling more precise regulation in areas with thick debris. However, empirical relationship approach 1 has some limitations. First, it underestimates the debris thickness at the terminus of an ice tongue because, as the layer of debris thickens, the conductive heat flux at the surface weakens. This causes LST saturation at some debris thicknesses [52], preventing accurate differentiation between thickness levels. Second, the model overestimates the size of the upper part of the ablation zone, primarily due to localized warming. This region, located at the confluence of the east and west glacier branches, experiences intense thermal radiation from the surrounding mountains, which is significantly greater than in other areas, especially on clear days [53]. Current models only account for the LST and overlook influential factors, limiting the comprehensiveness and accuracy of debris thickness estimation.

Future studies should consider the incorporation of multivariate variables closely associated with debris thickness [40]. Incorporating topographic, climatic, and other physical variables (e.g., glacier flow rate, debris supply rate, glacier melt) can mitigate the challenges posed by topographic complexity and climatic variability, enhancing the accuracy of estimation outcomes.

5.2. Analysis of Empirical Relationship Approach 2 Adaptation

Empirical relationship approach 2 offers significant advantages in estimating glacial debris thickness, being both straightforward and highly accurate for determining the debris distribution. Its precision is due to its use of the full dataset before fitting the empirical parameters, along with the uniform distribution of data across the study area. In contrast to empirical relationship approach 1, which bases its estimation on only two data points, empirical relationship approach 2 leverages the entire dataset, leading to more accurate results. However, like empirical relationship approach 1, empirical relationship approach 2 is constrained by its reliance on a single variable, the LST. Moreover, the accuracy of empirical relationship approach 2 is highly dependent on the quality of the dataset. Another limitation is that its empirical parameters are typically not generalizable, requiring recalibration to accommodate the specific meteorological conditions of different glacier regions [54]. During the ablation season, solar radiation causes significant diurnal variation in the LST (Figure 7). Empirical relationship approach 2, which uses remotely sensed transient imagery to invert the LST, may therefore have limited applicability beyond specific time windows, reducing the accuracy of the model parameters during other periods.

To address the problems outlined above, the following suggestions are proposed. (1) Normalize the LST for parameter fitting to reduce the errors caused by LST variation. Experimental results show that after normalizing the LST, the estimation accuracy improves significantly (

Table 3. Effect of the LST variation on the accuracy of empirical relationship approach 2.

	Satellite Transit Time		Average Temperature Change (°C)
	20100715 13:11	20130925 13:23	+3	−3	+5	−5
RMSE (cm)	43.82	51.94	38.68	48.04	73.65	51.40
Normalization RMSE (cm)	35.98	39.22	34.53

), demonstrating that this approach enhances the accuracy of empirical relationship approach 2 in debris thickness estimation. (2) Use higher-resolution imagery to model the relationship between the debris thickness and the LST. This would help reduce the errors arising from the mismatch between the raster scale of the imagery and the actual spatial distribution of debris, improving the overall precision of the estimates.

5.3. Analysis of Energy Balance of Debris Approach Adaptation

The energy balance of debris approach is the most suitable approach for mapping debris thickness in the range of 0.0 to 0.5 m. As an energy balance model rooted in physical mechanisms, it comprehensively incorporates climate and heat transfer factors, enabling higher prediction accuracy in regions with thin debris, and is not constrained by spatial or temporal limitations, making it both generalizable and robust. Although the energy balance of debris approach’s implementation requires a substantial number of observational inputs, the use of reanalyzed data can effectively address this challenge. Research has shown that reanalyzed data can supplement unavailable data in the energy balance of debris approach [28,55,56], thus greatly enhancing its applicability. This opens up the possibility of using the energy balance of debris approach to estimate the debris thickness in remote areas without the need for on-site data collection.

Despite its advantages, the energy balance of debris approach leads to significant underestimation of the debris thickness range and magnitude due to the following reasons. (1) It relies on an unreasonable assumption of a linear debris profile temperature. Figure 8 shows the complexity of the debris profile temperatures at the same time, showing how they vary with the debris thickness. The assumption of linearity holds true only for thin debris, but as the thickness increases, it can result in substantial underestimation. Rounce et al. [42] attempted to address this error by introducing a nonlinear approximation factor, but a single factor cannot accurately capture the complex heat transfer dynamics occurring throughout the debris layer across the entire ablation zone. (2) The thermal conductivity of the debris is influenced by several factors, including the density, rock thermal fusion, and porosity [57]. Simplifying the thermal conductivity to a single value reduces the accuracy of the predictions of the model. To improve the performance of the energy balance of debris approach, the following strategies are recommended. First, develop more accurate simulations of the temperature curve within the debris layer to better reflect the internal heat transfer processes. Second, conduct in-depth research on the relationship between the heat transfer coefficients and the physical properties of the debris to enhance the accuracy of these simulations.

This study has uncertainty in evaluating different debris thickness estimation approaches. First, the field measurements of debris thickness often do not align with the raster value space of satellite imagery, as the debris thickness can vary significantly within an area smaller than a single pixel. Second, field measurements of debris thickness are labor-intensive and time-consuming, with only older data available for comparison in many cases. However, due to the slow rate of change in debris, it is assumed that the area and thickness have remained stable over the short term, making the older data still valuable for analysis. Third, this study focused exclusively on the Koxkar Glacier, and the applicability of these results to other glaciers has yet to be validated.

6. Conclusions

In this study, we analyzed the distribution characteristics of the debris thickness on the Koxkar Glacier and compared the accuracy and applicability of three widely used estimation approaches. The key conclusions are as follows.

(1) The overall estimation results suggest that empirical relationship approach 2 has the highest estimation accuracy, while the energy balance of debris approach has the lowest. Although there are significant differences in the range of values, all three approaches effectively captured the distribution characteristics of the debris thickness.

(2) For different thicknesses, the energy balance of debris approach was the most suitable measurement approach, particularly for glaciers with an overall surface debris thickness of less than 50 cm. In contrast, the empirical LST and thickness approach was most effective in areas with thicker surface debris (greater than 100 cm). Among all three approaches, the estimation accuracy generally decreased as the debris thickness increased.

(3) Variations in topographic factors, particularly the elevation, significantly impact the accuracy of debris thickness estimates. The estimation accuracy of all three approaches increases with the elevation.

(4) By analyzing profile data, this study identifies the limitations of each approach and suggests potential directions for improvement, offering valuable guidance for enhancing the measurement accuracy and selecting suitable methods for debris thickness estimation in the future.

In conclusion, this study contributes to understanding the distribution characteristics of the debris thickness on the Koxkar Glacier and provides practical insights for choosing estimation approaches. Future research should focus on developing more accurate and streamlined methods to better support glacier mass balance calculations and modeling. This will also enhance our understanding of the effects of climate change on glaciers.