1. Introduction
Glaciers in High Mountain Asia, including the >10% with debris-covered surfaces (Herreid and Pellicciotti, Reference Herreid and Pellicciotti2020), have been losing mass during recent decades (e.g. Brun and others, Reference Brun, Berthier, Wagnon, Kääb and Treichler2017; Shean and others, Reference Shean, Bhushan, Montesano, Rounce, Arendt and Osmanoglu2020; Miles and others, Reference Miles, McCarthy, Dehecq, Kneib, Fugger and Pellicciotti2021). A thick debris mantle insulates the underlying ice, whereas a thin debris layer tends to enhance ice ablation (e.g. Østrem, Reference Østrem1959; Fyffe and others, 2014; Collier and others, Reference Collier, Maussion, Nicholson, Mölg, Immerzeel and Bush2015). It has been suggested that the insulating effect of debris reduces glacier mass loss; however, several studies have reported that debris-covered glaciers show comparable or greater thinning rates relative to clean glaciers, even when controlling for elevation biases (e.g. Nuimura and others, Reference Nuimura, Fujita, Yamaguchi and Sharma2012; Lamsal and others, Reference Lamsal, Fujita and Sakai2017; Brun, Reference Brun2019). This phenomenon is known as the debris-covered anomaly (Pellicciotti and others, Reference Pellicciotti, Stephan, Miles, Herreid, Immerzeel and Bolch2015; Salerno and others, 2017) and is generally attributed to a combination of melt enhancement and ice dynamics (e.g. Rounce and others, 2021). Previous studies have suggested that ice cliffs can cause local enhancement of debris-cover glacier surface melt rates (e.g. Sakai and others, Reference Sakai, Nakawo and Fujita1998; Buri and others, Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021; Miles and others, Reference Miles, Steiner, Buri, Immerzeel and Pellicciotti2022), contributing to the debris-cover anomaly, but the glacier-wide contribution of cliffs to mass balance is known for only very few sites (e.g. Buri and others, Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021).
Previous studies have attempted to quantify ice cliff backwasting (horizontal retreat) and melt rate (perpendicular ice loss) using in situ measurements (e.g. Sakai and others, Reference Sakai, Nakawo and Fujita1998; Steiner and others, Reference Steiner, Pellicciotti, Buri, Miles, Immerzeel and Reid2015; Anderson and others, Reference Anderson, Armstrong, Anderson and Buri2021a), terrestrial photogrammetry (e.g. Brun and others, Reference Brun2016; Watson and others, Reference Watson, Quincey, Smith, Carrivick, Rowan and James2017a; Kneib and others, 2022), airborne remote sensing (e.g. Brun and others, 2018; Mishra and others, 2022), and energy balance modelling (e.g. Han and others, Reference Han, Wang, Wei and Liu2010; Steiner and others, Reference Steiner, Pellicciotti, Buri, Miles, Immerzeel and Reid2015; Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016a, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeelb). In particular, the widespread use of Unoccupied/Uncrewed Aerial Vehicles (UAVs) and Structure from Motion (SfM) technology has dramatically advanced ice cliff studies (e.g. Brun and others, 2018; Immerzeel and others, 2014; Mishra and others, 2022;Zhao and others, 2023). Although the digital elevation model (DEM)-differencing approach can potentially derive precise mass loss of an ice cliff, this approach requires repeated photogrammetric surveys and consideration of ice dynamics, either directly or through mass conservation approaches that require measuring/modelling physical parameters (ice thickness and horizontal/vertical velocity distributions; Brun and others, 2018; Mishra and others, 2022; Zhao and others, 2023). This approach also lacks the temporal resolution to explain how ice cliffs fluctuate between the initial and subsequent timing of aerial photogrammetric surveys and provides only the total amount of melt. Other studies have developed energy balance models to quantify ice cliff mass loss (e.g. Sakai and others, Reference Sakai, Nakawo and Fujita1998; Reid and Brock, Reference Reid and Brock2014). Recent studies have combined UAV-based high-resolution terrain data with energy balance modelling (e.g. Steiner and others, Reference Steiner, Pellicciotti, Buri, Miles, Immerzeel and Reid2015; Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016a), expanding the potential for estimating melt patterns on the highly heterogeneous surfaces of ice cliffs. The energy balance approach can extend our knowledge of the temporal/sub-seasonal variability of ice cliff ablation and its mechanisms (i.e. the relationship with meteorological factors and ablation characteristics), rather than simply deriving mass loss. Such models have been shown to adequately quantify ice cliff melt at fine temporal scales and to be transferable across different climatic settings (Kneib and others, 2022). Buri et al. (Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021) first applied an energy balance model for ice cliffs that takes into account the cliff morphological changes to the catchment scale. Although these authors succeeded in quantifying the contribution of ice cliff mass loss to four glaciers in the Langtang Valley of the Nepal Himalaya, glacier-scale applications of ice cliff energy balance models are rare. Hence, the transferability of such a model to other regions with different climatic conditions needs to be confirmed. The reason why such studies are limited is because the application of an ice cliff energy balance model at the glacier scale requires high-resolution topographic data and an ice cliff inventory covering the entire debris-covered area to reconstruct the heterogeneous glacier surface. In addition, meteorological data (especially incoming shortwave and longwave radiation) over debris-covered glaciers at high altitude and the temperature characteristics of the supraglacial debris are also essential to force such process-based models. Furthermore, because previous studies focused on only short-term (Sakai and others, Reference Sakai, Nakawo and Fujita1998; Buri and others, Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021), the monthly fluctuations in ice cliff geometry and energy balance, and their relationship to morphological/energy-balancing characteristics on glacier-wide and annual scales, remain largely unknown.
In this study, we combine high-resolution aerial photogrammetric data and a process-based energy balance model to (1) characterise the monthly melt rate of ice cliffs, (2) estimate the annual ice cliff mass loss and its contribution to the glacier mass balance, and (3) quantify the contribution of different energy fluxes at the cliff surface on debris-covered Trakarding Glacier in the eastern Nepal Himalaya.
2. Study site, data and methods
2.1. Study site
Debris-covered Trakarding Glacier (27.9°N, 86.5°E) is located in Rolwaling Valley in the eastern Nepal Himalaya (Fig. 1a). This glacier spans elevations of 4530–6670 m above sea level (a.s.l.). The terminus of Trakarding Glacier is Tsho Rolpa, one of the largest glacial lakes in Nepal and a site of glacial lake hazard management. Trakarding Glacier covers an area of 8.21 km2, and 41% of its surface is covered by debris (Herreid and Pellicciotti, Reference Herreid and Pellicciotti2020). The glacier accumulates mass mainly from avalanches from the southeastern rock wall. Measured debris thickness on the glacier ranges from a few centimetres to ∼70 cm, with the maximum thickness being observed near the terminus. Trambau Glacier is situated above Trakarding Glacier, and the two glaciers have been disconnected since 1970 (Fig. 1a). Therefore, previous studies have tended to consider these glaciers together as the ‘Trambau–Trakarding glacier system’ (Sunako and others, Reference Sunako, Fujita, Sakai and Kayastha2019). We have been conducting in situ measurements on Trakarding Glacier since May 2016 using meteorological instruments and ablation stakes. We also conducted aerial photogrammetric surveys in October 2018 and October 2019 (Sato and others, 2021; Section 2.2), which defined the period of interest for the model simulation. This study investigates ice cliffs over a 2.93 km2 area of Trakarding Glacier (Fig. 1a, b, and d), covering most of the debris-covered part. We modified glacier outlines obtained from the GAMDAM glacier inventory (Nuimura and others, 2015; Sakai, Reference Sakai2019) to match the glacier width and terminus positions in October 2018.
Overview of the Trambau–Trakarding glacier system. (a) Location of the rolwaling region (inset) and outline of Trambau glacier, Trakarding glacier, and the study area; (b) an ice cliff in the debris-covered area; (c) the automatic weather station (AWS) beside Trakarding glacier; and (d) locations of air (ta) and debris surface (ts) temperature sensors, AWS, ablation stakes, time-lapse (TL) camera, and ice cliffs. glacier outlines in (a) are from the GAMDAM glacier inventory (Nuimura and others, 2015). Blue circles in (d) are air temperature sensors that were used to calculate the temperature lapse rate (LR; Fig. 2 and Table 1), and orange dots are temperature sensors for air and debris surface (Figs. 3, S1 and Table 1).
2.2. Terrain data and ice cliff inventory
We used high-resolution photogrammetry-based terrain data and an ice cliff inventory to compute cliff melt using an energy balance model (Section 2.3). Details of the datasets have been described by Sato and others (2021). The photogrammetric datasets were obtained on 18 October 2018 and 18–19 October 2019. We used a helicopter in 2018 and a fixed-wing UAV in 2019. These photogrammetric datasets covered most of the debris-covered area of the glacier. We generated high-resolution orthoimages and DEMs with 0.2 m resolution (hereafter SfM-orthoimages and SfM-DEMs) using the SfM software package Agisoft Metashape Professional Edition version 1.5.1. We conducted kinematic GPS surveys on-/off-glacier, and vertical uncertainties of these SfM-DEMs were estimated at ±1.82 m in 2018 (26,142 validation points) and ±2.35 m in 2019 (based on 8,790 points), respectively. From these SfM-orthoimages and SfM-DEMs, we manually delineated all ice cliffs and supraglacial ponds on Trakarding Glacier using ArcGIS. We found 481 ice cliffs in the 2018 orthoimages, covering an area of 1.38 × 105 m2 (Fig. 1d). The delineation uncertainty in individual cliff map-view areas (planimetric areas) was assessed by five operators and did not exceed 10% area of the cliff inventories. In addition, with respect to the sensitivity of the inclined area (actual area of ice cliff slope) to the quality of the DEM, the slope angle uncertainty of ±1° does not cause more than ±2% of inclined area changes (Sato and others, 2021).
2.3. Meteorological observations and debris surface temperature estimation
We obtained meteorological data from an automatic weather station (AWS) located beside Trakarding Glacier (4806 m a.s.l.; Fig. 1c and Table 1). Details of the meteorological observations have been described by Sunako et al. (Reference Sunako, Fujita, Sakai and Kayastha2019) and Fujita et al. (Reference Fujita, Sunako and Sakai2021). Air temperature, ground surface temperature, wind speed, relative humidity, and upward/downward shortwave radiation were recorded from May 2016 to November 2021. Air temperatures over the debris-covered area were also measured to calculate the temperature lapse rate. Unfortunately, the air temperature stations on the glacier were not operational during the period of interest (October 2018 to October 2019); therefore, we estimated the seasonal temperature lapse rate from air temperatures measured on the debris-covered area from March 2020 to October 2021 (Figs. 1d and 2, and Table 1). We divided this period into four seasons and calculated the seasonal/hourly temperature lapse rate following Heynen et al. (Reference Heynen, Miles, Ragettli, Buri, Immerzeel and Pellicciotti2016).
Details of meteorological instruments on Trakarding Glacier, the locations of which are shown in Fig. 1d. Details of meteorological observations at the AWS have been described by Sunako et al. (Reference Sunako, Fujita, Sakai and Kayastha2019) and Fujita et al. (Reference Fujita, Sunako and Sakai2021)
Temperature lapse rates over the debris-covered area of Trakarding glacier in the pre-monsoon (PRM; 1 march to 14 June), monsoon (M; 15 June to 30 September), post-monsoon (POM; 1 October to 30 November), and winter (W; 1 December to 28 February) seasons. The locations of the air temperature sensors are shown in Fig. 1d.
To calculate longwave radiation flux from heated debris, it is necessary to estimate the spatiotemporal surface temperature distribution, which we obtained on the basis of correlation with air temperature on the glacier (Foster and others, Reference Foster, Brock, Cutler and Diotri2012; Steiner and Pellicciotti, Reference Steiner and Pellicciotti2016). We measured air temperature (1.5 m above ground) and debris surface temperature at the same four locations (Fig. 1d and Table 1). These sensors were installed in May 2016, and they recorded air/debris temperature until November 2017 (Fig. S1). We then applied piece-wise linear regression between air temperature and debris surface temperature (Fig. 3). The air/debris temperature data were classified into four seasons, and the most appropriate piece-wise linear regressions were estimated from all temperature observation sites. Significant correlations were obtained for all seasons (Fig. 3). The resultant empirical equations were combined with the air temperature lapse rate to estimate the spatial distribution of debris surface temperature.
Relationship between air temperature and debris surface temperature in the pre-monsoon (a), monsoon (b), post-monsoon (c), and winter (d) seasons. Different colours indicate different locations (Figs. 1d and S1, and Table 1). Dashed vertical lines represent the threshold temperature (tc) in the piece-wise regression. Letter ‘r’ is the mean correlation coefficient between debris surface temperature and air temperature. The time series of air and debris surface temperatures are plotted in Figure S1.
2.4. Dynamic 3D ice cliff backwasting model
We employed a process-based dynamic 3D ice cliff backwasting model (hereafter; dynamic cliff backwasting model) developed by Buri et al. (Reference Buri, Miles, Steiner, Immerzeel, Wagnon and Pellicciotti2016b) to estimate year-round ice cliff mass loss. A full description of the model has been provided by Buri et al. (Reference Buri, Miles, Steiner, Immerzeel, Wagnon and Pellicciotti2016b), Buri and Pellicciotti (Reference Buri and Pellicciotti2018), and Buri et al. (Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021); therefore, only a brief summary is provided here. This grid-based model estimates the mass loss of ice cliffs by calculating the hourly energy balance on the cliff surface. The energy balance is calculated as follows:
\begin{equation}{Q_m} = {I_n} + {L_n} + H + LE{\text{ }}\end{equation}where Qm is the energy flux for cliff surface melt, In and Ln are net shortwave and longwave radiation fluxes, H is sensible heat flux, and LE is latent heat flux. Units of all fluxes are in W m–2. Previous studies neglected heat from precipitation and conductive heat flux into the ice cliff surfaces (e.g. Reid and Brock, Reference Reid and Brock2014). We also assume the conductive flux to be negligible when Qm is positive or zero. In some rare cases the energy balance can be negative, indicating either a non-negligible conductive heat flux or refreezing. In these cases, we treat Qm as zero in the melt calculation, and we discuss this effect later (Section 4.5).
The dynamic cliff backwasting model uses a high-resolution DEM to calculate radiation fluxes that account for complex cliff slope geometry and debris-covered surface topography. The ice cliff and glacier surface topography were reconstructed from SfM-DEM-2018 (Section 2.2), and the valley/mountain topography surrounding the calculation area was obtained from ASTER-GDEM3. The SfM-DEM was resampled to 1.0 m resolution owing to the computational cost. Then, two ice cliffs were excluded from the simulation because they were too small and contained only a few grid cells, resulting in 479 ice cliffs being used as model input. These topographic datasets allowed the calculation of sky- and debris view angles for each grid cell on the cliff surface to estimate net radiation. Incoming shortwave radiation consists of direct shortwave radiation (Is), diffuse shortwave radiation from the sky (Ds), and incoming shortwave radiation reflected from terrain (Dt). The incoming shortwave radiation observed at the AWS was split into Is and Ds following Reindl et al. (Reference Reindl, Beckman and Duffie1990). This approach uses a clearness index, which is the ratio of incoming shortwave radiation at the AWS to theoretical extraterrestrial solar radiation, which has been used in previous cliff modelling studies (Han and others, Reference Han, Wang, Wei and Liu2010; Reid and Brock, Reference Reid and Brock2014). Incoming longwave radiation comprises longwave radiation from the sky (Ls) and from surrounding debris (Ld). These longwave fluxes were modelled using the Stefan–Boltzmann relation forced by air and debris surface temperature (Section 2.3). Finally, the cliff volume loss (Vcl; m3 w.e.) was calculated as follows:
\begin{equation}{V_{cl}} = \tfrac{{{Q_m} \cdot \Delta t \cdot S}}{{{\rho _i}{L_f}}}\end{equation} where 
$\Delta t$ is the simulation time step (3600 s), S is a map-view area of ice cliff (m2), ρi is ice density in the debris-covered area (900 kg m−3), and Lf is the latent heat of fusion of ice (334 kJ kg−1).
Furthermore, the dynamic cliff backwasting model incorporates changes in ice cliff geometry associated with its ablation. This allows a realistic long-term energy balance to be calculated for ice cliffs, which are known to morphological change considerably over time (Watson and others, Reference Watson, Quincey, Carrivick and Smith2017b; Kneib and others, 2021; Sato and others, 2021; Kneib and others, 2022). In our study, we set nine geometry updates over one year. The first geometry update was on 28 February 2018, the end of the winter season, followed by geometry updates on the last day of each month until the end of the simulation period (18 October 2019). The geometry updates consider expansion due to cliff backwasting and shrinkage/burial due to debris slumping. Some ice cliffs were completely buried and disappeared during the simulation period. In this study, the angle at which the ice cliff surface was buried by debris (slope threshold) was set at 35° after some trials comparing model output and observed ice cliff shape at the end of the simulation period in October 2019 (Sato and others, 2021). As previous studies have reported that the angle of repose for debris mantle varies between 30° and 45° (Sakai and others, Reference Sakai, Nakawo and Fujita2002; Kraaijenbrink and others, Reference Kraaijenbrink, Shea, Pellicciotti, de and Immerzeel2016; Moore, Reference Moore2018; Westoby and others, 2020; Sato and others, 2021), this threshold is considered reasonable. Further detailed model parameters/settings are given in Table S1. This study is the first attempt to apply this dynamic model to an entire glacier on a year-round time scale. When we calculated annual cliff melt rates, the initial/final ice cliff areas largely changed in this dynamic cliff backwasting model. Hence, the daily melt rate of each ice cliff surface (Md; m w.e. d−1) was defined as follows:
\begin{equation}{M_d} = \frac{{2\Delta V}}{{\left( {{S_1} + {S_2}} \right) \cdot \Delta t}}\end{equation} where 
$\Delta $V is the cumulative cliff volume loss in one year (m3 w.e.), S1 and S2 are the initial and final ice cliff areas during the simulation period (m2), and 
$\Delta t$ is the number of days (365 days). In this study, ice cliff melt rate is defined as the rate of vertical ablation at each grid cell.
In the case where ice cliffs are covered in snow, the energy flux on the cliff surface should be used to melt the snow layer without melting the ice cliffs. To determine the effect of snow on ice cliff melt, we examined daily time-lapse photographs taken near the glacier terminus to assess snow cover, obtained using a Brinno TLC200 time-lapse camera that has a 112° field of view (Fig. 1d and S2a–d). Photographs were taken each day at 13:46 (Nepal time, UTC + 5:45) and were used to estimate the number of snow-covered days in each month. We assumed that all ice cliffs had zero melt for an entire day if snow cover was visible at the glacier terminus. We applied the number of snow cover days in the target month (dsnow; days) to reduce the monthly ice cliff volume loss (Vcl, month; m3 w.e.) accounting for snow cover as follows:
\begin{equation}{V_{cl,{\text{ }}month}} = \frac{{{V_{cl,{\text{ }}est}} \cdot \left( {{d_{month}} - {d_{snow}}} \right)}}{{{d_{month}}}},\end{equation}where Vcl, est is the monthly ice cliff volume loss without considering snow cover (m3 w.e.) and dmonth is the number of days in the target month (∼30 days). Then we defined overestimated ice cliff volume loss (Vcl, overest; m3 w.e.), that would not have occurred due to the presence of snow as follows:
\begin{equation}{V_{cl,{\text{ }}overest}} = {V_{cl,{\text{ }}est}} \cdot \frac{{{d_{snow}}}}{{{d_{month}}}}\end{equation}2.5. Model validation and sensitivity analysis
The dynamic cliff backwasting model was developed based on studies at Lirung Glacier in the Langtang Valley, Nepal Himalaya, and has shown good performance against stake- and photogrammetry-based cliff melt measurements, though these have been limited to four specific cliffs (Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016a, Reference Buri, Miles, Steiner, Immerzeel, Wagnon and Pellicciotti2016b). The model has also shown an excellent ability to reproduce weekly melt rates at other monsoonal debris-covered glaciers (Kneib and others, 2022). We did not install stakes on the cliff surface at Trakarding Glacier; therefore, we validated the model performance by comparing the modelled ice cliff morphology with the observed ice cliff morphology from aerial photogrammetry. First, we detected more than 200 ice cliffs (hereafter ‘surviving ice cliffs’) that remained between the 2018 and 2019 manually digitised inventories (Sato and others, 2021). Then, we randomly selected 100 of these surviving ice cliffs, avoiding merged/split ice cliffs. We compared the cliff map-view area, cliff inclined area (actual slope surface), slope, and aspect of the modelled cliffs at the end of the modelling period with those of the observed cliffs. The cliffs were categorised as stream- or pond-influenced when located closer than 40 m to either of these features (Kneib and others, 2023), and this validation exercise was conducted for each type of ice cliff.
We also used the same 100 selected cliffs to estimate the sensitivity of the model to its parameters. We ran the cliff backwasting model 22 times, individually increasing/decreasing the five fixed physical and six meteorological parameters (Table S2). We change surface roughness, debris/ice emissivity, and debris/ice albedo as physical parameters. We chose air temperature, incoming short-/long-wave radiation, debris surface temperature, relative humidity, and wind speed as meteorological parameters. These meteorological parameters were obtained at AWS (Figs. 1d and Table 1) and estimated/interpolated on each cliff using temperature lapse rates and correlation between debris and air temperature (Figs. 2 and 3). Owing to our higher spatial resolution (1.0 m instead of 3.0 m) and longer modelling period (year-round rather than five months), we compared the results of sensitivity tests only on upper/lower parameter bounds taken from Buri et al. (Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021) to reduce the computational cost. The parameters and ranges used in the one-at-a-time sensitivity test were taken from Buri et al. (Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021) using the same dynamic cliff backwasting model. In the previous study, 100 Monte Carlo simulations were run on 40 ice cliffs to estimate the model uncertainty. We also tested the model using a coarser resampled DEM (3.0 m) to assess the effect of resolution.
2.6. Glacier mass balance and mass loss contribution of ice cliffs
To estimate the mass loss contribution of ice cliffs to the glacier-scale mass balance, we employed the surface elevation change provided by Hugonnet and others (2021), as our photogrammetric datasets do not cover the entire glacier. The glacier-wide mass balance on Trakarding Glacier was calculated using a modified GAMDAM glacier inventory (Section 2.1) and the surface elevation change rate during 2015–2019. The uncertainty of this glacier thinning rate has been reported as ±0.17 m a−1 in the target region (Hugonnet and others, 2021). The glacier-averaged surface elevation change rate was converted to glacier mass balance (m w.e. a−1) using the bulk density of ice in debris-covered areas (0.9) relative to the density of water (Miles and others, Reference Miles, Willis, Buri, Steiner, Arnold and Pellicciotti2018). We also calculated the mass balance of the debris-free Trambau Glacier to evaluate the contribution of ice cliff mass loss to the mass loss of the whole Trambau–Trakarding glacier system. Then, we used 0.85 as the density of ice (accounting for firn zones) to calculate the glacier-wide mass balance (Huss, Reference Huss2013). These glacier-wide mass balances indicate net ablation, including accumulation. Hence, we calculated mass loss only in the ablation zone of the debris-covered tongue (Fig. S3) of Trakarding Glacier to evaluate the contribution of ice cliff melt to only debris-covered part. To estimate mass loss on the debris-covered tongue, we calculated the mean emergence velocity (Ve, m a−1) from the upper and lower ice flux as follows (e.g. Nuimura and others, Reference Nuimura, Fujita, Fukui, Asahi, Aryal and Ageta2011; Miles and others, Reference Miles, Willis, Buri, Steiner, Arnold and Pellicciotti2018; Fig. S3):
\begin{equation}{V_{\text{e}}} = \frac{{\left( {{q_{{in}}} - {q_{{out}}}} \right)}}{A},\end{equation}where qin and qou t are the ice fluxes at the upper and lower boundaries (m3 a−1; Fig. S3), respectively, and A is the area of the target zone (m2). The ice fluxes q (qin and qout, m3 a−1) were calculated by
\begin{equation}
q = W \cdot h \cdot v \cdot\ \cdots\end{equation}where W is the width of the ice flux gate (m), h is ice thickness (m), and v is depth-averaged glacier velocity (m a−1). We employed the published ice thickness (h) from Farinotti and others (2019) (Fig. S3). The depth-averaged velocity (Vc) was assumed to be 90% of the surface velocity estimated from photogrammetry-based orthoimages (Sato and others, 2021), following previous studies (e.g. Miles and others, Reference Miles, Willis, Buri, Steiner, Arnold and Pellicciotti2018; Sato and others, Reference Sato, Fujita, Inoue and Sakai2022). We then set simple upper/lower flux gates to cover most of the debris-covered area (Fig. S3). Finally, the mass balance of the debris-covered tongue (MB, m w.e.) was calculated as follows:
\begin{equation}MB = \left( {\frac{{dh}}{{dt}} - {V_e}} \right)\cdot\frac{{{\rho _i}}}{{{\rho _w}}},\end{equation}where dh/dt is the mean surface elevation change rate on the debris-covered tongue (m a−1), and ρw is the density of water (1000 kg m−3).
Trakarding Glacier is a lake-terminating glacier; therefore, its terminus has retreated and been losing mass through calving and frontal ablation. Almost all previous studies dealing with ice cliff melt effects have focused only on land-terminating debris-covered glaciers (e.g. Brun and others, 2018; Buri and others, Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021; Zhao and others, 2023). To compare our results with these previous studies, we estimated terminus volume loss and eliminated it from the glacier mass balance. Terminus volume loss Vt (m3) can be theoretically determined from the change in ice front position and ice flux at the terminus (e.g. Wei, Reference Wei2021), as follows:
\begin{equation}{V_t} = {V_f} + {V_r},{\text{ }}\end{equation}where Vf and Vr are volume losses (m3) derived from terminus ice flux and terminus retreat, respectively. Although the lake depth at the terminus is required to calculate Vf and Vr, the elevation change dataset used to calculate the glacier-wide mass balance in this study (Hugonnet and others, 2021) does not account for glacier mass loss below the lake surface. Therefore, only the mass loss above the lake is needed to exclude the calving effect from the glacier mass balance. We calculated ice flux above the lake surface (Vf) at the glacier terminus as follows:
\begin{equation}{V_f} = {W_t}\cdot{h_l}\cdot{v_t},{\text{ }}\end{equation}where Wt is the glacier width at the terminus (m), hl is the ice thickness height above the lake surface (m), and vt is the depth-averaged velocity at the terminus (m a−1). We used photogrammetry-based DEM and surface velocity data for the study period to compute hl and vt (Sato and others, 2021). We also used elevation change above the lake surface of the retreating terminus portion area obtained from photogrammetry-based DEMs to calculate the terminus volume loss (Vr).
Finally, we calculated the enhancement factor (EF) (e.g. Brun and others, 2018; Buri and others, Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021), which measures the relative cliff melt contribution to the relative cliff area for the debris-covered tongue, Trakarding Glacier, and Trambau–Trakarding glacier-system-scale mass loss. The EF is defined as follows:
\begin{equation}EF = \frac{{\frac{{{V_{cl}}}}{{{A_{cl}}}}}}{{\frac{{{V_{gl}} - {V_{cl}}}}{{{A_{gl}} - {A_{cl}}}}}},{\text{ }}\end{equation}where Vgl and Vcl are the volume losses, and Agl and Acl are the planimetric areas of the glacier surface and ice cliffs, respectively.
3. Results
3.1. Model output and sensitivity
We compared modelled and observed ice cliff morphology for the selected 100 surviving ice cliffs used for validation (Fig. 4). The mean residuals in observed ice cliff shape and model output were 165 m2 for the map-view area and 154 m2 for the inclined area (Fig. 4a and b), corresponding to 32% of the mean cliff map-view area and 22% of the inclined area, respectively. Considering the mean perimeter length of the ice cliffs (143 m), the mean residual of the map-view area represents a 1.1-pixel overestimate of the outward expansion of the geometry of the cliffs in the dynamic cliff backwasting model. Excessive expansion of the cliff area can lead to overestimation of cliff mass loss; conversely, it can also lead to underestimation of the yearly ice cliff melt rate [Eq. (3)] The model tends to make the cliffs downwasting (surface lowering) more than they would in reality, resulting in shallower slopes after one year (mean residual of −4.5°; Fig. 4c). However, there is no notable difference in the ice cliff aspect (mean residual of 3.0°; Fig. 4d). The errors of these validations were aggregated in cliff aspect and stream/pond-influenced (or not) type ice cliffs (Fig. 4a–d). However, there were no significant differences between categories.
Comparison of observed (obs.) and modelled (mod.) ice cliff shapes for 100 surviving ice cliffs in October 2019, showing map-view area (a), inclined area (b), slope difference (c), and aspect difference (d) of ice cliffs. The red outlines in (a) and (b) indicate stream- or pond-influenced ice cliffs (section 2.5), and plot shapes/colours show original ice cliff aspects in October 2018. The R-squared values in (a) and (b) are from linear regressions between observed and modelled ice cliff map-view/inclined areas. Boxes in (c) and (d) show the 25th–75th percentiles, and black dots depict individual errors on the modelled ice cliff slope and aspect. In (c) and (d), ‘stream/pond’ (‘normal’) indicates whether ice cliffs are influenced (or not) by a stream/pond.
Examples of modelled ice cliffs extracted from the stagnant part (surface velocity less than 3 m a−1; Figs. 5 and S4) show how cliff geometry changes with each computational step (Fig. 5a and b). The model outputs seem to reconstruct the morphological changes on the ice cliff that have become less steep (Fig. 5b). Furthermore, as the simulation was performed at a high spatial resolution (1.0 m), the spatial heterogeneity of incoming radiation due to the complex morphology and the melt gradient within individual ice cliffs are well represented (Figs. 5c–e and S4).
Examples of cliff backwasting model outputs. (a) Outlines of changes in ice cliff geometry and (b) changes in elevation profiles of ice cliff slope, compared with observations. (c–e) Monthly averaged values of (c) melt rate, (d) incoming direct shortwave radiation, and (e) longwave radiation from surrounding debris in May 2019. The black arrow across the cliff in (a) indicates the elevation profile in (b). The colours of the ice cliff outlines and elevation profiles in (a) and (b) depict the respective geometry updates. Pink and grey polygons in (a) are observed ice cliff shapes. Hashed lines in (b) are observed ice cliff elevation profiles from Sato and others (2021). Arrows in (c) indicate aspects of the two example ice cliffs.
The model sensitivities were tested on 12 parameters (23 patterns considering upper/lower bounds) for the 100 surviving ice cliffs. We calculated the relative residual in each cliff volume loss for the original parameters and the upper/lower bounds of the tested parameters (Fig. S5a–c). For almost all tests in the chosen range, the median absolute deviation in cliff volume loss did not exceed 15%, except for the incoming longwave radiation and DEM resolution tests (Fig. S5a–c). When the input DEM was resampled to 3.0 m resolution, the mean residual of relative volume difference was +82%, which is a substantial increase in volume loss (Fig. S5c). This remarkable increase in volume loss is not due to an increase in the melt rate of ice cliffs but rather to unrealistic expansions in cliff size during the simulation period (Fig. S6).
3.2. Volume loss of ice cliffs
We calculated the volume loss of 479 ice cliffs spread over the entirety of Trakarding Glacier with 1 hour time steps. The total volume loss of the ice cliffs during one year (18 October 2018 to 18 October 2019) was 1.34 × 106 m3 w.e. The monthly ice cliff volume losses of the entire glacier are summarised in Fig. 6. The largest volume loss was observed in June (2.32 × 105 m3 w.e., 17% of the annual volume loss). The monthly volume loss of ice cliffs increased after March, at the end of the winter season. Monthly cliff volume losses from May to August (summer monsoon season) each exceeded 12% of the annual volume loss. We were able to clearly confirm snow cover in time-lapse images between January and April (Fig. 6). Most of the days in February and March were snow-covered, which reduced the annual volume loss of ice cliffs by 10% compared with the case where snow cover was ignored (i.e. potential volume loss). Our simulation of the whole-year range showed that 18% of the annual mass loss occurred outside of the ablation season (November to April; Fig. 6).
Monthly ice cliff volume loss and its ratio to total cliff volume loss. The unfilled dashed bars from January to April represent overestimated ice cliff volume loss without considering snow cover. The October 2018 period is from the 18th to the 30th (12 days), and the October 2019 period is from the 1st to the 18th (18 days). The POM, W, PRM, and m indicate post-monsoon, winter, pre-monsoon, and monsoon season (grey shaded areas depict winter and monsoon season).
We divided the ice cliff volume loss into eight aspect bins (Figs. 7a and S7a). The largest volume loss was observed in the north-facing cliff group (337.5°–22.5°), with 2.88 × 104 m3 w.e. (21% of the total volume loss). The north-, northeast-, and northwest-facing cliff groups account for >50% of the total volume loss of all ice cliffs on Trakarding Glacier. However, it should be noted that the area of ice cliffs was not uniform across aspect bins, and north-facing ice cliffs predominated (Fig. 7b).
Aspect distribution of (a) ice cliff volume loss from October 2018 to October 2019 and (b) ice cliff map-view area.
3.3. Ice cliff melt rate
Cliff melt rate was calculated using monthly/annual time steps for each ice cliff. The mean and median melt rate of all ice cliffs on Trakarding Glacier were 2.7 × 10−2 and 2.6 × 10−2 m w.e. d−1, respectively. The ice cliff melt rate varied, ranging from 1.1 × 10−2 to 7.3 × 10−2 m w.e. d−1 for the 479 ice cliffs. We estimated the relationship between ice cliff melt rate and cliff aspect/elevation (Fig. 8a and b). Then, we applied harmonic and linear regression to each parameter and found that both parameters had a significant relationship with ice cliff melt rate (p < 0.001). However, ice cliff aspect better explained the melt rate (R 2 = 0.48; Fig. 8a) than did cliff elevation (R 2 = 0.03; Fig. 8b). In this harmonic regression between the aspect and melt rate of the ice cliffs, the maximum melt rate occurred for south-facing cliffs (3.9 × 10−2 m w.e. d−1 at 181°), 1.8 times higher than the minimum value for north-facing cliffs (2.1 × 10−2 m w.e. d−1 at 1°; Figs. 8a and S7b).
Individual ice cliff melt rates over the simulation period as a function of (a) cliff aspect and (b) cliff elevation. dashed red lines in (a) and (b) are harmonic/linear fittings of ice cliff melt rate. The colour scale indicates the volume loss at individual ice cliffs from October 2018 to October 2019.
We also calculated monthly ice cliff melt rates and aggregated them with respect to these aspects (Fig. 9 and Table S3). The results showed clear seasonality in ice cliff melt rate: the highest melt rate was observed in June (median of 6.5 × 10−2 m w.e. d−1), and the lowest in December (3.4 × 10−3 m w.e. d−1), with an estimated difference of almost 20-fold between them. The ice cliff melt rate continued to decrease from the start of the simulation until December, recording its minimum value. From January approaching the ablation season, the melt rate increased, reaching its peak in June (Fig. 9).
Monthly ice cliff melt rate plotted over time with respect to ice cliff aspect. Coloured boxes show the 25th–75th percentiles, and black dots represent outliers. Red asterisks indicate a significant difference in melt rate between north- and south-facing ice cliffs in each month (p < 0.05; student’s t-test).
Notably, we detected a relationship between aspect and the seasonal dependence of ice cliff melt rate. During the cold season (October 2018 to February 2019), melt rates show a strong aspect dependence, with the melt rate decreasing in order from south- to west- to east- to north-facing cliffs (Fig. 9). However, the difference in ice cliff melt rate between cliffs with different aspects decreased approaching the melt season. From June to August 2019, there was no significant difference in melt rate between north- and south-facing cliffs (via Student’s t-test; p > 0.05). In June, when the maximum melt rate was observed, the melt rate was almost equal for all ice cliff aspects (Fig. 9). Focusing on the contrasting south- and north-facing ice cliffs, the strongest aspect dependence of melt rate (relative difference) was in December, when the median value of the melt rate of south-facing ice cliffs (1.9 × 10−2 m w.e. d−1) was more than 20 times higher than that of north-facing ice cliffs (8.5 × 10−4 m w.e. d−1). In absolute terms, the difference between the melt rates of the south- and north-facing ice cliffs reached its maximum at the end of October 2018, when the melt rate of south-facing ice cliffs (3.6 × 10−2 m w.e. d−1) was 2.7 × 10−2 m w.e. d−1 higher than that of north-facing ice cliffs (9.0 × 10−3 m w.e. d−1). It should be noted that this period (October 2018) does not cover the whole month (Table S3).
3.4. Energy balance
We aggregated the component of the energy balance that leads to ice cliff ablation each month (Fig. 10a and Table S4). Fig. 10a shows the daily mean energy balance for each month. It should be noted that almost allmelt of ice cliffs occurred in daytime (between 08:00 and 17:00). In addition, during the timesteps when the energy balance became negative, Qm was treated as zero (no melt; Steiner and others, Reference Steiner, Pellicciotti, Buri, Miles, Immerzeel and Reid2015). In terms of the contribution to ice cliff melt, incoming shortwave radiation (Iin) was the most significant component in the energy balance, which consists of direct shortwave radiation (Is) and diffuse shortwave radiation from the sky (Ds) and terrain (Dt) (Fig. 10a). The annual mean of In is 163 W m−2, which explains 63% of the melt component (positive energy flux for melt). Iin shows seasonal variation between ∼ 100 and 250 W m−2, and the highest Iin was observed in May (249 W m−2). There was also seasonal variability in the ratio of Iin components (Fig. 10d). The annual mean energy fluxes of Is and Ds were 63 and 11 W m−2 (55% and 39% of Iin), respectively (Table S4). Is was generally more predominant than Ds, although Ds became larger in July and September, accounting for more than 50% of the total incoming shortwave radiation (Fig. 10d). Dt was the minor component and remained below 10% relative to Iin throughout the target period (Fig. 10d).
Mean energy balance of all (a), north-facing (b), and south-facing (c) ice cliffs. Components of incoming shortwave radiation (is, ds, and dt) of all (d), north-facing (e), and south-facing (f) ice cliffs. (g) Observed air temperature (ta, AWS) and incoming shortwave radiation (swin, AWS) at the AWS. (h) Clearness index at the AWS (section 2.4) and normalised solar altitude for Trakarding Glacier. Note that almost of the melt of ice cliffs were occurred in daytime. During the time steps when the energy balance became negative, the heat for ice cliff melt (qm) was treated as zero.
We compared the monthly Qm with monthly air temperature and shortwave radiation observed at the AWS (Fig. 10g and Table S5). A significant correlation was found between Qm and daytime shortwave radiation (r = 0.97, p < 0.001) but not between Qm and air temperature at the AWS. For incoming longwave radiation (Lin), the mean of longwave radiation fluxes from the sky (Ls) and debris (LD) were 135 and 102 W m−2 (57% and 43% of Lin) during the target period, respectively (Table S4). In this study, debris surface temperature was forced by air temperature (Section 2.3; Fig. 3); even considering the topography (debris view angle) for individual ice cliffs, there was a significant correlation between monthly AWS air temperature and Ld (r = 0.77, p < 0.05). Mean net longwave radiation (Ln) was constantly negative throughout the year, which indicates a consistent negative influence on cliff ablation at Trakarding Glacier (Fig. 10a). Sensible and latent heat fluxes either inhibited or promoted cliff melt depending on the season but did not exceed 10% of the positive/negative energy fluxes in most months (Fig. 10a).
The energy balance shows contrasting patterns depending on aspect, especially between north- and south-facing ice cliff groups. The energy balance components for the north-/south-facing cliff aspects are summarised in Fig. 10b and c. The annual averages for north- and south-facing Qm were 76 and 138 W m−2, respectively (Table S4). The average Qm for south-facing ice cliffs is 1.8 times higher than that for north-facing ice cliffs. This result supports the harmonic regression of the annual melt rate estimated from all ice cliffs (Fig. 8a). Is values of north- and south-facing ice cliffs are notably different, with the annual mean of Is of south-facing ice cliffs (133 W m−2) being ∼ 2 times higher than that of north-facing ice cliffs (Fig. 10b–c and Table S4). This difference in Is contributes to the contrasting Qm between north- and south-facing ice cliffs. The mean annual ratios of Is and D s for north-facing cliffs were 49% and 44%, and those for south-facing cliffs were 65% and 30%, respectively (Table S4). The aspect dependence of the Iin component amplified/decreased according to season, being more pronounced in the cold season (Fig. 10e and f). The average Is from December to February, defined as the winter season in this study, was 26 W m−2 for north-facing ice cliffs and 115 W m−2 for south-facing ice cliffs. During the winter season, the Is of north-facing cliffs accounted for less than 35% of In, whereas south-facing cliffs accounted for ∼ 70% (Fig. 10e and f). However, the difference in shortwave radiation flux components decreased towards the monsoon season, with Is accounting for ∼ 50% of Iin for both north- and south-facing ice cliffs (Fig. 10e and f).
3.5. Mass loss contribution and melt enhancement of ice cliffs
We estimated the glacier mass balance and compared it with the mass loss of ice cliffs (Table 2). It is noted that these results consider the effect of snow cover, which reduced ice cliff mass loss in the cold season (Section 3.2). The glacier-wide mass balance of Trakarding Glacier in 2015–2019 was −0.86 ± 0.15 m w.e. a−1. Comparing the glacier-scale and ice cliff mass losses, ice cliffs accounted for 23.2% of the total glacier mass loss, despite covering only 1.7 % of the glacier surface. When the terminus mass loss (Vt) was excluded from the total glacier mass loss, the total ice cliff mass loss accounted for 25.8% of the glacier-scale mass loss (Table 2). The EF [Eq. (11)] of the ice cliffs relative to Trakarding Glacier was 15.6 (including Vt) and 20.3 (excluding Vt). With respect to the mass loss of the debris-covered tongue (ablation zone), ice cliffs covered 5.1% of the surface area and contributed 33.1% of the glacier-scale mass loss (excluding Vt; Table 2). The EF of ice cliff melt was 9.1 on the debris-covered tongue. We also compared our modelled ice cliff ablation with mass loss measured at seven stakes in the debris-covered area (Figs. 1d and 11). The mean annual mass balance at stakes was −1.2 m w.e. (ranging from −2.2 to −0.6 m w.e.), and the mean debris thickness at the stake points was 15 cm (12–65 cm) (Fig. 11). The mean annual melt rate of the ice cliffs (−9.9 m w.e.) was 8.3 times higher than the average measured for debris-covered ice. From these results, the ice cliff EF on the debris-covered tongue that was estimated from remote-sensing data (9.1) is considered reasonable compared with the observed EF (8.3).
Annual ablation of ice cliffs (coloured crosses) and mass balance at stakes (coloured circles) from October 2018 to October 2019. The box plots are ice cliff ablation rates summarised for 100 m elevation bins. Cross colours indicate ice cliff aspect, and circle colours denote debris thickness (DT; cm) measured at stakes. Grey boxes show the 25–75th percentiles, and unfilled dots represent outliers.
Summary of area, mass loss, and enhancement factor (EF) of ice cliffs on the debris-covered tongue, Trakarding glacier, and the Trambau–Trakarding glacier system. The terms ‘incl.’ (‘excl.’) vt indicate that the terminus volume loss is contained (not contained) in glacier mass loss. Glacier mass losses was calculated from the surface elevation change dataset of Hugonnet and others (2021). The target regions of the debris-covered tongue, Trakarding glacier, and the Trambau–Trakarding glacier system are shown in figures 1a and S3. Note that the mass loss of ice cliffs on snow-covered days was excluded from all results
We also estimated the contribution of ice cliff mass loss to that of the Trambau–Trakarding glacier system (Table 2). The ice cliffs had an area relative to the glacier system of 0.4% and accounted for 7.3% (including Vt) and 7.8% (excluding Vt) of the total glacier system mass loss. In addition, the cliff EF in the glacier system was 18.0 (including Vt) and 19.3 (excluding Vt) (Table 2). The cliff EF for the Trambau–Trakarding glacier system was smaller than that of only Trakarding Glacier when calving volume loss was excluded (Table 2).
4. Discussion
Previous studies have estimated ice cliff melt using energy balance approaches (e.g. Sakai and others, Reference Sakai, Nakawo and Fujita2002; Han and others, Reference Han, Wang, Wei and Liu2010; Buri and others, Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021) for the main ablation season (May to October). However, our simulation is the first to estimate full-year and glacier-scale ice cliff melt. Our simulation also considers cliff geometry updates, which is the only effective method for estimating long-term ice cliff melt by reconstructing realistic cliff fluctuations (Buri and others, Reference Buri, Miles, Steiner, Immerzeel, Wagnon and Pellicciotti2016b). Our study is the first attempt to apply such a dynamic model to another site that has different conditions from the Langtang region. Furthermore, this study is the first to use high-resolution (1.0 m) terrain data to reconstruct the morphology of ice cliffs as input with validated its geomorphological change with as many as a hundred ice cliffs. Here, we discuss the implementation of our modelling approach, the role of ice cliffs as melt hotspots on Trakarding Glacier, and the seasonal melt characteristics of ice cliffs.
4.1. Performance of the dynamic cliff backwasting model
We validated the cliff shape of the modelled output and showed that of all metrics (size, slope, and aspect; Fig. 4a–d). We showed that the cliff area was slightly overestimated by the dynamic model (1.1-pixel overestimate of the outward expansion; Sect. 3.1). As the morphology of the ice cliffs was updated in a non-linear manner every month, it is difficult to quantify the volume loss uncertainty associated with the overestimation of the ice cliff area expansion at the final output. Here, we assume that the overestimation of ice cliff area expansion occurred linearly during the simulation period, which means that the overestimation of ice cliff area changes could cause up to a 16% overestimation of cliff volume loss. Furthermore, cliff slope was underestimated by the model (Fig. 4c). The lack of representation in the model of physical collapse and thermal undercutting caused by channels and ponds might have contributed to this change in ice cliff morphology beyond that expected from the energy balance. Such changes in ice cliff morphology have been observed on other debris-covered glaciers (e.g. Kraaijenbrink and others, Reference Kraaijenbrink, Shea, Pellicciotti, de and Immerzeel2016; Kneib and others, 2022; Petersen and others, Reference Petersen, Hock and Loso2024). The median residuals of slope angle between stream/pond-influenced and normal (not stream/pond-influenced) ice cliffs were −5.7° and −3.1°, respectively (Fig. 4c), indicating no significant difference. However, this lack of difference may be due to the limited sample size (only 21 normal ice cliffs) and to the difficulty of assessing the strength of the link between hydrology and cliff evolution from remote-sensing data only (Kneib and others, 2023). Petersen et al. (Reference Petersen, Hock and Loso2024) reported that stream-influenced ice cliffs have incisions formed by supraglacial channels in their underpart and steeper slopes than those of ice cliffs that are not stream influenced at Kenicott Glacier, Alaska. This stream effect also leads to the generation of wider angles between cliff slopes and backwasting ramps (the opposite side mound of the cliff surface). Although the dynamic cliff backwasting model used in this study does not take into account such an effect of supraglacial channels, this effect should be incorporated into future dynamic models to estimate more realistic ice cliff fluctuations.
We tested the sensitivity of the model to 12 physical and meteorological parameters (22 patterns) for the 100 ice cliffs used for validation (Fig. S5a and b). The uncertainty in our simulation does not appear to differ from that of previous work: the sensitivity analysis shows a similar response to parameters as reported by Buri et al. (Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021), using the same dynamic cliff backwasting model. Buri et al. (Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021) reported an uncertainty of 3355 m3 w.e. per cliff in the cliff backwasting model for the Langtang region based on 100 Monte Carlo simulations. In the work of Buri et al. (Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021) in the Langtang region, the source ice cliff inventory was delineated from SPOT-6 satellite imagery (1.5 m resolution; Steiner and others, Reference Steiner, Buri, Miles, Ragettli and Pellicciotti2019), and those authors reported that the median extracted ice cliff size was 845 m2, seven times larger than the median ice cliff size of Trakarding Glacier (112 m2). Application of the uncertainty in volume loss of ice cliffs estimated by Buri et al. (Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021) to the median cliff size of Trakarding Glacier yields an uncertainty of volume loss per cliff of 445 m3 w.e in this study. Multiplication of the uncertainty per ice cliff by 479 ice cliffs yields a value of 16% of the total volume loss of all ice cliffs on Trakarding Glacier (2.21 × 10−4 km3 w.e.).
We also tested the effect of the resolution of the input DEM, and this simulation showed a remarkable increase in mass loss compared with the baseline setting due to cliff expansion (Fig. S5c). The resolution of the DEM is one parameter that significantly influences the estimated expansion of ice cliffs (Fig. S6). If coarse satellite-based DEMs (2–3 m) are used in this dynamic model, the parameters suppressing the expansion/shrinkage of ice cliffs should be optimised for each case. In particular, two parameters should be considered: the slope threshold, which controls cliff reburial, and the negative buffer size, which controls the shrinkage and expansion of cliffs (Buri and Pellicciotti, Reference Buri and Pellicciotti2018).
4.2. Ice cliffs as melt hot spots of Trakarding Glacier and the Trambau–Trakarding glacier system
Our study indicates that ice cliffs were responsible for a substantial portion of mass loss at Trakarding Glacier (Table 2). Although ice cliffs covered only 1.7% of the surface of Trakarding Glacier, they contributed 25.8% of the glacier-scale mass loss when the terminus mass loss was considered. Buri et al. (Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021) used the same energy balance model and estimated the contribution of ice cliffs to glacier-scale mass loss for four glaciers in Langtang Valley. Those authors reported that the contribution of cliffs to glacier-scale mass loss ranged from 7% (Shalbachum Glacier) to 17% (Langtang Glacier). Our simulation used a higher-resolution original orthoimage/DEM (0.2 m/0.2 m spatial resolutions; Sato and others, 2021) compared with the Langtang study (1.5 m/3.0 m; Steiner and others, Reference Steiner, Buri, Miles, Ragettli and Pellicciotti2019) to generate the ice cliff inventory. Hence, more small cliffs were extracted at Trakarding on account of the difference in resolution of images for extracting ice cliffs (Kneib and others, 2020; Section 4.1). In addition, we targeted a longer period (1 year) than that of Buri et al. (Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021), who computed the melt only for the ablation season. The greater contribution of ice cliff mass loss to glacier-scale mass loss at Trakarding Glacier relative to the Langtang Valley glaciers can be attributed to the following reasons. Trakarding Glacier has a higher density of ice cliffs compared with the glaciers in Langtang Valley (Steiner and others, Reference Steiner, Buri, Miles, Ragettli and Pellicciotti2019; Sato and others, 2021; Kneib and others, 2023). In addition, the mass loss in the debris-covered area might be effectively compensated by avalanche accumulation from the headwall, suppressing the mass loss on the glacier surface (Fig. S8). The upper part of Trakarding Glacier, which has progressively separated from the clean Trambau Glacier during the past few decades, has no widespread transition zone from the debris to clean areas (Fig. S8), where melting should be enhanced (e.g. Fyffe and others, 2014; Fyffe and others, Reference Fyffe, Woodget, Kirkbride, Deline, Westoby and Brock2020). These glacier characteristics appear to have resulted in a higher melt contribution of ice cliffs to glacier-scale mass loss at Trakarding Glacier than at other debris-covered glaciers. Quantifying avalanche accumulation and tracking terminus mass loss will be essential aspects of future work to understand the glacier regime.
Considering the entire Trambau–Trakarding glacier system, ice cliffs accounted for 0.4% of the total area and contributed 7.3% of the glacier-scale mass loss (Table 2, including Vt). Although Trambau Glacier has debris-covered medial moraines and contains some ice cliffs (Kneib and others, 2023), these were not covered by our photogrammetric survey. If these ice cliffs were included, then the ice cliff mass loss and their contribution to the mass loss of the glacier system should be increased. The EF of ice cliffs for the glacier system was 18.0, which is similar to that for Trakarding Glacier (20.3). In general, the EF of ice cliffs is greatly increased when upper clean ice areas are considered, as accumulation by solid precipitation suppresses the glacier mass loss (Buri and others, Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021). In contrast, decreased EF values occur in the case of extended ablation zones. In our target region, the EF of Trakarding Glacier changed little when including the accumulation zone of Trambau Glacier (Table 2). Sunako et al. (Reference Sunako, Fujita, Sakai and Kayastha2019) reported that the equilibrium line altitude of the Trambau–Trakarding glacier system is ∼5800 m a.s.l. and that most of the clean ice part is in the ablation zone. In addition, the mass balance of the lower part of Trambau Glacier (the clean part) is more negative (approximately −5 to −4 m w.e. a−1) than that of Trakarding Glacier (Sunako and others, Reference Sunako, Fujita, Sakai and Kayastha2019). As Trambau Glacier (the clean part) has a large ablation zone, the EF of the entire glacier system might therefore show a similar value to the EF of Trakarding Glacier despite including an accumulation zone (Table 2). Nevertheless, the mass loss of ice cliffs relative to that of the entire glacier system (>7%) can be considered significant in terms of total mass balance or discharge.
4.3. Aspect characteristics of ice cliff melting
Previous studies have addressed the aspect characteristics of cliff melt and reported higher ablation rates on south- to southeast-facing ice cliffs using energy balance approaches (Sakai and others, Reference Sakai, Nakawo and Fujita1998; Buri and Pellicciotti, Reference Buri and Pellicciotti2018). Compared with these experimental studies, our study estimated ice cliff melting rates in a more realistic setting, in which we confirmed a clear dependence of ice cliff melt rate on aspect (Figs. 8a, 11 and S7b). We estimated a harmonic regression between ice cliff aspect and annual melting rate, with the maximum value for south-facing ice cliffs, which was almost twice as high as the lowest value (for north-facing ice cliffs; Fig. 8a). This aspect dependence is ascribed to the difference in incoming direct shortwave radiation between south- and north-facing ice cliffs (Fig. 10b, c, e, and f and Table S4; Section 3.4).
Despite south-facing ice cliffs having higher melt rate, the south-facing ice cliffs accounted for only about 20% of the total volume loss, more than 50% of melt occurred in the population of north-facing ice cliffs (Figs. 7a, S7a, and S7b). This volume loss pattern is attributed to the persistence of north-facing ice cliffs, resulting in their larger proportion relative to other aspects ice cliffs on Trakarding Glacier (Fig. 7b). Previous studies have shown that cliff lifespan is strongly related to cliff aspect. South-facing ice cliffs receive intense shortwave radiation on their cliff top, leading to a gentler slope, so they tend to be buried by debris (Sakai and others, Reference Sakai, Nakawo and Fujita1998, Reference Sakai, Nakawo and Fujita2002; Buri and Pellicciotti, Reference Buri and Pellicciotti2018). Such a characteristic of incoming shortwave radiation and melt pattern was also confirmed in our simulations (Fig. 5c and d). Conversely, north-facing ice cliffs provide shade on their own surface, thereby extending their lifespan (Sakai and others, Reference Sakai, Nakawo and Fujita2002). The predominance of north-facing ice cliffs has been confirmed on Trakarding Glacier (Fig. 7b; Sato and others, 2021) and has also been observed in other Himalayan regions (e.g. Watson and others, Reference Watson, Quincey, Carrivick and Smith2017b; Kneib and others, 2023). Therefore, south-facing ice cliffs, despite having a much higher melt rate than that of north-facing cliffs, are in the minority and do not exceed the total mass loss of north-facing ice cliffs.
4.4. Seasonal fluctuations in ice cliff melt
Our full-year simulation of 479 cliffs provided seasonal fluctuations and aspect characteristics of ice cliff melt. The monthly ice cliff melt rate was not correlated with air temperature but with shortwave radiation observed at the AWS, and the highest melt rate and total volume loss occurred during the pre-monsoon season (Figs. 6, 9, and 10).
In addition, the aspect dependence of ice cliff melt rates was magnified during the winter season but dissipated towards the monsoon season (Fig. 9). These characteristics are ascribed to regime changes in direct shortwave radiation, which accounts for the majority of the energy flux of cliff melt (Fig. 10a; Section 3.4). In the Himalayan region, the solar elevation angle lowers in winter (i.e. the solar zenith angle increases; Cooper, Reference Cooper1969; Fig. 10h). This leads to a strong contrast in the energy balance of ice cliffs, with substantially more direct shortwave radiation being received on south-facing cliffs than on north-facing ice cliffs (Fig. 10b, c, e, f, and h; Section 3.4). This contrast of incoming shortwave radiation leads to a pronounced difference in winter melt rate between north- and south-facing ice cliffs (Fig. 9). Conversely, the solar elevation angle becomes higher from the end of the pre-monsoon to the monsoon season, which provides more direct shortwave radiation to a broader range of ice cliff aspects (Fig. 10e, f, and h). In addition, during the monsoon season, diffuse shortwave radiation is predominant due to cloud cover shelter in the Himalaya (Fig. 10a; Section 3.4; e.g. Sakai and others, Reference Sakai, Nakawo and Fujita2002). Thus, the aspect dependence of ice cliff melt rate is expected to disappear during the monsoon season, and the melt rate decreased despite higher temperatures than during the pre-monsoon (Figs. 9 and 10g). This study, therefore, revealed the aspect dependence of ice cliff melt rate (i.e. dependence on direct shortwave radiation) and confirmed that south-facing ice cliffs maintain a high melt rate despite the cold season (Figs. 8a, 9, and 10b–c). Such cliff melts in the cold season should thus be considered in estimating debris-covered glacier melt for glaciers with abundant cliffs at their surface.
4.5. Limitations and future work
This study did not consider ice cliffs that were newly formed during the simulation period. At Trakarding Glacier, almost half of all ice cliffs (235 cliffs; 2.4 × 10−2 km2) were replaced between 2018 and 2019 (Sato and others, 2021). The mean individual area and total area of these newly formed ice cliffs were relatively small compared with surviving ice cliffs. However, the newly formed ice cliffs exhibited a wide variety of cliff aspects, and the proportion of south-facing newly formed cliffs was higher than that of surviving ice cliffs. We estimated the melt of newly formed cliffs by applying the empirical relationship derived between ice cliff melt rate and aspect to the inventory of newly formed ice cliffs detected by Sato and others (2021) (Fig. 8a). The newly formed ice cliffs had a slightly higher mean melt rate (2.8 × 10−2 m w.e. d−1) than that ice cliffs used in simulation (2.7 × 10−2 m w.e. d−1). When considering mass loss of newly formed ice cliffs, the monthly mass loss of the total ice cliffs increased by 16%. It is impossible to determine when the ice cliffs formed during the one-year simulation period. If we hypothesised that cliffs were formed in the middle of the simulation period (i.e. the annual mass loss of newly formed ice cliffs was reduced by half), then newly formed ice cliffs accounted for only 2% of glacier-scale mass loss. Although it is worth considering the presence of newly formed ice cliffs to estimate the actual contribution of ice cliffs to debris-covered glaciers, this may not be particularly significant on Trakarding Glacier.
Although trials were conducted involving changing model parameters, it was challenging to reconstruct disappeared ice cliffs (i.e. those buried by debris during the simulation period) completely. We compared simulated and observed ice cliffs (Sato and others, 2021), focusing on cliff disappearance (Table S6). Although the model usually represents surviving cliffs well, only 25% of cliffs that should have disappeared at October 2019, disappeared during the simulation period. To estimate ice cliff persistence, not only energy balance should be considered but also changes in surface albedo (e.g. Kneib and others, 2022), cliff calving (e.g. Miles and others, 2017; Watson and others, Reference Watson, Quincey, Smith, Carrivick, Rowan and James2017a), debris mobility (e.g. Westoby and others, 2020), and maintenance effects from adjacent supraglacial ponds/channels (e.g. Anderson and others, Reference Anderson, Armstrong, Anderson, Scherler and Petersen2021b; Kneib and others, 2023; Petersen and others, Reference Petersen, Hock and Loso2024). Although these processes are difficult to represent, future models of cliff evolution should attempt to incorporate them. The temperature distribution of the debris surface, which controls incoming longwave radiation, may also influence cliff persistence (Sakai and others, Reference Sakai, Nakawo and Fujita2002). Previous studies have reported that debris surface temperature is controlled by incoming shortwave radiation and surface aspect (e.g. Kraaijenbrink and others, 2018). In this study, we estimated debris surface temperature from the empirical relationship with air temperature. However, further geometric factors with respect to debris-covered surfaces should be considered when estimating surface temperature for the dynamic model. In terms of glacier scale cliff dynamics, one previous study has developed an ice cliff tracking method to detect their dynamics automatically (Kneib and others, 2021). Such work has the potential to be further automated, possibly with machine learning approaches, which may allow model development to also predict ice cliff dynamics.
This study assumed that all ice cliff melting ceased for entire days when snow cover was observed from time-lapse images (Fig. S2). These time-lapse images covered only the terminus portion of Trakarding Glacier, and the distance of time-lapse camera from the debris-covered surface limited the information that could be obtained from the images. Therefore, the effect of snow cover on individual ice cliffs could not be considered. Girona-Mata et al. (Reference Girona-Mata, Miles, Ragettli and Pellicciotti2019) estimated snow line altitude fluctuation at Langtang catchment and reported that snow cover on south-facing slopes melts/retreats more efficiently than that of north-facing slopes owing to incoming shortwave radiation flux, especially in winter. If such mechanisms were applicable to the microtopography of debris-covered glaciers (ice cliff surfaces), then snow cover on north-facing ice cliffs may persist longer than that on south-facing ice cliffs, inhibiting ice cliff melt. This would lead to the aspect dependence of cliff melt rate in cold seasons becoming even more pronounced.
While our simulation did not consider the refreezing of cliff melt water, the incoming energy flux sometimes became negative in the cold season and/or night time. Some previous studies suggested the importance of taking into account refreezing at cliff surfaces (Steiner and others, Reference Steiner, Pellicciotti, Buri, Miles, Immerzeel and Reid2015). Following the Steiner et al. (Reference Steiner, Pellicciotti, Buri, Miles, Immerzeel and Reid2015) approach, we estimated the energy available for refreezing at the cliff surfaces. Using the aggregated hourly energy balance for each month (Fig. S9), we first detected the time that energy balance transitioned from positive to negative values. Then, we calculated the net energy for ice melt one hour before/after the energy balance became negative, assuming that all the negative energy could be used for refreezing, but only in the first hour with energy balance, when there would still be water at the cliff’s surface. These results were used to estimate the negative energy flux available for daily refreezing (Qrf; W m−1) and are summarised in Table S7. The relative contribution of refreezing was significant during the winter season. When refreezing was taken into account, the amount of ice cliff melt in December was suppressed by 18.5% (Table S7). However, the contribution of refreezing is small in terms of ice cliff volume loss because of the low ice cliff volume losses observed during the winter season (Fig. 6). From June to September, the mean hourly energy flux does not become negative, and no refreezing is expected. When cumulative Qrf is compared to the cumulative Qm, its annual melt suppression effect (ice replenishment) was counted as less than 1% of total ice cliff melt (Table S7). Although the year-round scale contribution of refreezing was estimated to be close to negligible, the processes of meltwater discharge and water retention on ice cliff slopes are so complex that more observation and quantification of this component might be required for future development of the dynamic cliff backwasting model.
Although we identified a clear aspect dependence of ice cliff melt rate at Trakarding Glacier, there are cases in other regions where no aspect dependence has been observed. Anderson et al. (Reference Anderson, Armstrong, Anderson and Buri2021a) conducted in-situ measurements of ice cliff backwasting rates and found no aspect dependence at debris-covered Kennicott Glacier in Alaska. Hence, those authors concluded that a simple degree-day factor approach could be applied to estimate cliff mass loss across the entire glacier. It would be valuable to apply our approach to other glaciers/regions in order to understand such regional differences in the cliff melt process.
5. Conclusion
This study presented an application of a process-based dynamic ice cliff backwasting model for Trakarding Glacier in the eastern Nepal Himalaya. The study represents the first attempt to estimate full-year-scale cliff mass loss using high-resolution photogrammetry-based terrain data, cliff inventory, and cliff geometry updates. The main conclusions can be summarised as follows:
1. Variations in ice cliff melt rate are more strongly influenced by aspect than by elevation, and the melt rate of south-facing ice cliffs is ∼2 times higher than the rate of north-facing ice cliffs. However, the predominance of north-facing cliffs results in that these cliffs still account for the majority of ice cliff melt.
2. Between July and September, diffuse shortwave radiation accounted for more than 50% of incoming shortwave radiation to ice cliff surfaces. The amount of shortwave radiation flux is strongly dependent on ice cliff aspect with seasonal variability. During the winter season, the direct shortwave radiation of south-facing ice cliffs reached 4.5 times that of north-facing ice cliffs.
3. The aspect dependence of ice cliff melt rate increased during the cold season, but there was no significant aspect-related dependence during the monsoon season. Such seasonal changes are ascribed to changes in solar altitude and monsoonal cloud cover, which are strongly related to direct shortwave radiation to the cliff surface.
4. Although ice cliffs cover less than 2% of the glacier surface, they account for ∼26% of the glacier surface mass loss on the entirety of Trakarding Glacier. Even at the scale of the Trambau–Trakarding glacier system, ice cliffs are a non-negligible melt component, accounting for >7% of the total mass loss of the system, despite occupying <1% of the surface area.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.17.
Acknowledgements
This study was supported by JSPS-KAKENHI (grant numbers 17H01621, 18KK0098, 22H00033, 23KK0064), JSPS-SNSF Bilateral Programmes project (HOPE, High-elevation precipitation in High Mountain Asia; grant numbers 183633 for Swiss, 20191503 for Japan), JSPS-Overseas Challenge Program for Young Researchers (grant number 202180207), the Sasakawa Scientific Research Grant from The Japan Science Society. We are indebted to Guide For All Seasons Trek for logistical support during the fieldwork. We thank H. Inoue, R. B. Kayastha, R. Kayastha, A. Tsushima, E. A. Podolskiy, S. Takenaka, and R. Hazuki for their supports during the field observations. We also thank C. Fyffe, T. Shaw, M. McCarthy, C. Zhao, S. Fugger, and A. Jouberton for their help and guidance. We acknowledge C. Tijm-Reijmer and the two anonymous reviewers for their comments that helped improve the manuscript.
Author contributions
YS, PB, MK, ESM, FP, and KF designed this study. YS, SS, TRG, AS, and KF conducted the field observations. YS analyzed the remote-sensing data. PB and MK provided the ice cliff energy balance model and helped apply it. YS, SS, and KF analyzed meteorological and debris temperature data. YS, PB, MK, ESM, and FK wrote the manuscript. All of the authors contributed to the discussion.
Competing interests
The contact author has declared that neither they nor their co-authors have any competing interests.
Open access
