Modeling Surface Processes on Debris-Covered Glaciers: A Review with Reference to the High Mountain Asia

: Surface processes on debris-covered glaciers are governed by a variety of controlling factors including climate, debris load, water bodies, and topography. Currently, we have not achieved a general consensus on the role of supraglacial processes in regulating climate–glacier sensitivity in High Mountain Asia, which is mainly due to a lack of an integrated understanding of glacier surface dynamics as a function of debris properties, mass movement, and ponding. Therefore, further investigations on supraglacial processes is needed in order to provide more accurate assessments of the hydrological cycle, water resources, and natural hazards in the region. Given the scarcity of long-term in situ data and the difﬁculty of conducting ﬁeldwork on these glaciers, many numerical models have been developed by recent studies. This review summarizes our current knowledge of surface processes on debris-covered glaciers with an emphasis on the related modeling efforts. We present an integrated view on how numerical modeling provide insights into glacier surface ablation, supraglacial debris transport, morphological variation, pond dynamics, and ice-cliff evolution. We also highlight the remote sensing approaches that facilitate modeling, and discuss the limitations of existing models regarding their capabilities to address coupled processes on debris-covered glaciers and suggest research directions.


Introduction
Debris-covered glaciers are common in high-mountain regions, and they play a critical role in governing local water resources, natural hazards, mountain geodynamics, and landscape evolution [1][2][3][4][5][6]. Studies have identified and discussed a variety of controlling factors that govern the dynamics of a debris-covered glacier system including climate, topography, debris load, ice-cliffs, and supraglacial water bodies [7][8][9][10][11]. Supraglacial debris plays a fundamental role in governing the dynamics of a debris-covered glacier because debris load can significantly influence radiation transfer and energy balance at the glacier surface [4,10,[12][13][14][15][16][17]. The thermal properties, porosity, and moisture content variations of debris loads are also coupled with near-surface microclimate which regulates sub-debris melting rate and the glaciers response to precipitation [18,19]. Unfortunately, the properties and transport of supraglacial debris have not been adequately represented in modeling efforts, which poses a challenge to our current understanding of debris-covered glaciers and magnifies the uncertainty in current climate-glacier sensitivity assessments [3,9,11,20,21]. For example, most models characterize debris insulation that reduces ablation, recent observations and simulations, however, show that some debris-covered glaciers have exhibited rapid thinning despite the presence of supraglacial debris [8,22,23]. This "debriscovered glacier anomaly" [24] suggests that the ablation rates for debris-covered glaciers are similar in magnitude to debris-free glaciers, and it appears that debris-covered glaciers could be more sensitive to climate change than previously thought.  [4,38,40,41]. Glacier outlines were acquired from the GLIMS (Global Land Ice Measurements from Space) database.
Supraglacial debris plays a crucial role in glacier surface ablation and morphological evolution [10,13,15]. Supraglacial debris thickness can be spatially heterogeneous due to debris transport by ice-flow, gravitational movement, and hydrological processes [10,29,48]. Scaledependent sediment transport processes also control glacier topography, mass balance gradient, and glacial drainage [49,50]. Although several approaches have been developed to estimate the distribution of debris thickness (see, e.g., in [46]), to understand the effects of debris on surface ablation (see, e.g., in [30,46,51]), and to model debris transport on glacier surfaces (see, e.g., in [9,10]), there are still many unknowns associated with supraglacial debris.

Properties of Supraglacial Debris
Headwall erosion, hillslope sediment transport, basal processes, and mass transport from moraines have been identified as the main sources of supraglacial debris [10,[52][53][54]. Mass movement processes from the adjacent terrain can be significant over time, including landslides [21], avalanches [55], and rockfalls from hillslopes and headwalls [21]. Glacial erosion at the glacier bed also encapsulates debris into the englacial debris load which may eventually emerge at the surface due to advection [10,56]. For many valley glaciers in the Himalaya, headwall and sidewall erosion is considered as the primary source of supraglacial debris [53,57]. Oversteepened hillslopes on headwalls and surrounding sidewalls can deliver large amounts of rocks and sediments onto glacier surfaces and are transported to lower altitudes due to gravity, ice-flow, or supra-fluvial processes. A typical headwall erosion rate ranges from 0.5 to 2 mm year −1 [10]. Hillslope erosion rates are also regulated by the relief and slope angles, and topography controls rock resistance to erosion and the locations where debris can be deposited onto the glacier surface [10]. Moraines and tributary glaciers also contribute to debris input [52,54].
Many studies have documented a negative relationship between debris thickness and ablation rate due to insulation by a thick debris layer (see, e.g., in [12][13][14]). These studies, however, did not account for the melt-enhancing effect of thin debris due to the higher absorption of solar radiation by debris and conduction of heat into the underlying ice [14,58,59]. Therefore, debris-covered ice may exhibit high melt rate than bare-ice if the debris thickness is thin enough. The "critical debris thickness" at which the sub-debris ablation rate equals that of bare-ice usually ranges from 0.02 to 0.1 m [29].
For many temperate mountain glaciers, debris loads over the terminus can limit glacial recession [3,10,52,60] because debris loads are much thicker at the terminus due to increased surface melting at lower altitudes, and the accumulation of the upwelling englacial debris load due to the upward flow of ice in the lower glacier [42,61,62]. Several studies have confirmed that a heavily debris-covered terminus allows the glacier to extend further below the ELA than a clean-ice glacier [3,10].
Several approaches have been developed to estimate the distribution of debris thickness. Mihalcea et al. [46] and Juen et al. [40] derived debris thickness from land surface temperature based on regressions equations using field measurements. Zhang et al. [30] use surface thermal resistance computed from satellite date as a proxy for debris thickness. Rounce et al. [45] developed an integrated method to estimate debris thickness, which involves DEM differencing, ice-flow modeling, and energy balance modeling. It is important to note that these approaches may be site-specific and only yield instantaneous debris thickness distribution that corresponds to the acquisition time of the satellite imagery.
Laboratory experiments on the debris effect were reported by Reznichenko et al. [59], in which the diurnal cycles of radiation forcing and rainfall were found to regulate the melting of debris-covered ice, such that the reduction in ablation rate due to debris cover occurs only when the diurnal cycles of radiation forcing are present. They also found that rainfall on a highly permeable thin debris can accelerate ablation due to effective heat transfer, relatively impermeable debris, however, reduces ablation rate as the rain may freeze within the debris during night before reaching the ice.

Transport of Supraglacial Debris
Understanding the timescales, pathways, and magnitude of debris transport is critical for estimating the variation in supraglacial debris thickness, which significantly governs the mass balance of a debris-covered glacier [10,56]. Modeling debris supply and debris flux is difficult as it requires accurate characterizations of ice-flow dynamics that is responsible for supra-and englacial debris advection [9,10,56], and gravitational and erosion processes that control the supply from headwalls, hillslopes, bedrocks, and moraines [21,54]. These processes must also be allowed to co-evolve and be temporally coincident with the variations in debris properties (such as thickness, thermal resistance, and albedo) that govern ablation. Due to this complexity, the transport of debris is usually neglected or oversimplified in glacier simulations. Nevertheless, several models have been developed to gain insights into these important processes [10,52,54,56].
Ice flow transports debris from active erosion zones to glacier terminus and depositional basins via supraglacial fluxes, englacial fluxes, and basal fluxes [21,60]. Supraglacial debris transport can be investigated using remote sensing approaches, for example, Gibson et al. [29] estimated a debris flux rate range of 0 to 12,000 m 3 a −1 for the Baltoro Glacier based on the analysis of multi-temporal satellite images combined with debris thickness estimates from satellite surface temperature product.
A recent ice-flow-debris transport model was developed by Anderson and Anderson [10]. They modeled englacial and supraglacial debris advection under a steady debris input to understand the mechanisms in the debris-glacier-climate system. Simulations indicated that debris has significant control on glacier length and gradients of ice discharge, ice thickness, and surface velocities. Their model demonstrated that high debris flux slows down the glacier and contributes to extending its length. The basis of this model describes the transport of supraglacial debris in the ice-velocity field as an advection process [10,63]: where h d is the thickness or concentration of surface debris, t is time, and u s is surface iceflow velocity. The use of advection is valid for glacial debris transport because advection is defined as the transport of materials due to the bulk motion of a fluid, and glacier ice is a form of a viscoelastic fluid that is the basis for all modern ice-flow models, although a more accurate model could be achieved by adding a diffusion term to the advection equation [56]. It is worth to mention that englacial debris transport also contribute to the increase in supraglacial debris load as the englacial debris is usually advected towards the surface in the ablation zone and can emerge on the surface within a century as demonstrated in the simulations by Anderson and Anderson [10] and Wirbel et al. [56]. The model presented by Wirbel et al. [56] characterizes 3-D englacial debris advection with the 3-D ice-flow solved using a full-Stokes approach, which is a major improvement to the 2-D model by Anderson and Anderson [10]. The performance of the model was tested using an actual glacier and a hypothetical glacier, given an estimate of debris sources, this model is capable of simulating the timescale, location and the amount of englacial debris that emerge on the surface, which acts as a bridge between englacial and supraglacial debris loads.
Gravitational debris movement (such as sliding and slumping) occurs on hillslopes and on glacier surfaces [21]. Field observations have identified sediment sliding or slumping off steepening ice cliffs due to ice-cliff retreat and supraglacial pond expansion [27,36]. Therefore, gravity-driven debris flux can play a significant role in local debris thickness redistribution especially during the ablation season when the surface topography is constantly changing under rapid melting. Most existing ablation models, however, assume static debris thickness during the simulation period, and existing approaches for deriving realistic debris thickness distribution suffer from high uncertainties [29,40,46], which limits the accuracy of modeled ablation rates.
Another group of models characterizes the mass transport from moraines on debriscovered glaciers. For example, Anderson [52] modeled the evolution of medial moraines on debris-covered glacier tongues and suggested that the thick moraines may finally evolve into a wide thin debris layer in the lower ablation zone. van Woerkom et al. [54] studied the sediment transport from the lateral moraine using high resolution DEM generated from 5 years of UAV data collected over the debris-covered tongue of the Lirung Glacier in Nepal. Their analysis showed that lateral moraines contribute to debris thickening along the margin of the glacier surface; however, the moraines supply only addresses debris thickness increase at glacier margins, the full debris balance of the entire glacier tongue is still unsolved. Better characterizations of debris transport and properties including thickness, grain size, mineralogical mixtures, and water content are required for future studies [15,64].

Surface Ablation
Surface ablation dominates the mass loss of most valley glaciers [8,37]. Although thick debris loads generally suppress melting, ablation rates for a debris-covered glacier can still be high [8,22]. Several studies have reported or predicted high-magnitude downwasting of debris-covered glaciers in the Himalayas [11,23,29]. Researchers also found that the average ablation rate for a glacier to be similar between debris-covered ice and clean ice in the Hindu Kush-Karakoram-Himalaya. Immerzeel et al. [23] presented basin-scale simulations over the Baltoro Glacier, which predicted a significant increase in total runoff, downwasting, and retreat throughout the twenty-first century despite thick debris cover. Simulations by Fujita et al. [8] suggested that the debris cover on a Himalayan glacier can promote ablation such that the debris-covered area contributes more per unit area to the total runoff than the clean-ice portion. Recent studies also found that ice cliffs and supraglacial lakes significantly accelerate differential melting on debris-covered glaciers [3,9,25,47,65].
Studies have also indicated that the spatial distribution of ablation rates on a debriscovered glacier is highly heterogeneous [9,25,33,66]. For example, ablation rate measurements at the same elevation on the Baltoro Glacier can range from about 0.02 md −1 to 0.06 md −1 [66]. The spatial pattern of ablation rates corresponds to the highly variable debris depth distribution over the glacier surface [29,46], which highlights the importance of accounting for debris effects in a surface ablation model.

Surface Energy Balance
The energy balance at the glacier surface is the basis for a physics-based surface ablation model in which solar radiation is the main driving force for ice melt. For a bare-ice glacier, the energy balance is focused on the air-ice interface, for a debris-covered glacier; however, both the air-debris interface and the debris-ice interface need to be taken into account [15,67]. An accurate modeling of the surface energy balance depends on the characterization of multiple surface features, and usually requires inputs and calibrations from in situ measurements or remote sensing analysis (e.g., surface temperature, albedo, debris depth, and surface DEM). The most widely used energy balance model at the debris-air interface can be written as [15,30,41,67] where Q s is the net shortwave radiation flux, Q l is the net longwave radiation flux, Q h is the net sensible-heat flux, Q e is the net latent-heat flux, and Q c is the conductive heat flux into the debris which governs the ablation rate. Note that in a distributed model, heat flux terms may need to be computed at the beginning of each iteration due to the changing surface conditions such as albedo and surface temperature [16,51]. Net radiation fluxes dominate glacier surface energy balance and significantly control surface ablation rates [15,37]. The net shortwave radiation flux and the net longwave radiation flux can be computed as where E b is the direct beam irradiance from the sun, E d is the diffuse skylight irradiance, E t is the adjacent terrain irradiance, and α is the surface albedo. For the longwave radiation fluxes, ε a , ε s , and ε t are the emissivity for the air, glacier surface, and adjacent terrain, respectively; σ is the Stefan-Boltzmann constant; T 0 is the air temperature which is a function of altitude; T s is the glacier surface temperature; and T t is the surface temperature of the surrounding terrain, which represents the longwave adjacent terrain irradiance that is not usually accounted for due to the high complexity and uncertainty. However, studies have suggested adjacent terrain irradiance can make a non-negligible contribution to the highmagnitude melt on some ice cliffs [36,44], the longwave component can also be significant as the temperature of debris surface can be higher than 30 • C during the day [44,68]. Glacier surface albedo regulates net shortwave irradiance and generally decreases over time during the ablation season [37]. It is challenging to model the spatiotemporal variability in surface albedo on a debris-covered glacier as that would require detailed knowledge about the composition, moisture content, thermal, and reflectance properties of supraglacial debris.
The dominant surface irradiance component is the direct beam irradiance (E b ). Many surface ablation models only account for the direct beam irradiance under assumptions of simplified solar geometry, surface albedo, and topographic effects (see, e.g., in [16,30,51]). These assumptions often yield unrealistically homogeneous patterns of ablation rates over glaciers. Complex topography in the High Mountain Asia causes significant variations in surface irradiance over short distances due to extreme relief and cast shadows [8]. The widely accepted parameterization scheme for E b can be written as where E 0 is the exo-atmospheric irradiance [69], which is modified by the Earth-Sun distance and orbital parameters; λ is the wavelength of light; T ↓ is the downward atmospheric transmittance; cosi is the cosine of the incidence angle of illumination (i); and S is a binary coefficient that accounts for cast shadows. E b is a function of wavelength range, and 0.3-3 µm is considered as the typical range for computing broadband shortwave irradiance [70]. Atmosphere attenuation (represented as downward transmittance T ↓ ) accounts for the absorption and scattering of the solar beam [71,72]. Atmospheric transmittance is highly wavelength dependent and controlled by spatio-temporal variations in various atmospheric constituents including aerosols, and water vapor [72,73]. Cloud conditions may also have an impact on the magnitude of direct beam irradiance [74]. The diffuse skylight irradiance E d is governed by both meso-scale and local topographic properties, and is found to have a non-negligible effect on ice-cliff melting [44]. Dozier and Frew [75] developed a simplified model for computing E d : where E dh is the isotropic diffuse skylight irradiance for a horizontal surface defined by Bird and Riordan [69], and V d is the skyview factor as defined by Dozier and Frew [75]. While a more computational expansive solution of the diffuse skylight irradiance was presented Proy et al. [76]. Local topography also effects the magnitude of irradiance, which is usually characterized by the incidence angle of illumination (i) between the sun and normal to the ground surface [77]: where θ s is the apparent solar-zenith angle, θ t represents the terrain slope angle, φ s is the solar-azimuth angle, and φ t is the terrain slope-azimuth angle. Therefore, the calculation of cosi requires digital elevation model (DEM) and solar geometry model. Mesoscale topographic condition of the glacier valley, which is usually characterized by sky-view factors, also has significant impact on surface energy balance by affecting the net radiation fluxes and surface temperature [68,78]. Figure 2 is an example showing the comparison between modeled and measured shortwave irradiance variation over a day at the same location on the Baltoro Glacier in the Karakoram. The measured magnitudes are greater mostly due to the omission of adjacent-terrain irradiance in the modeled result. Turbulent fluxes are difficult to estimate as factors like wind speed can be difficult to obtain over glaciers [79]; nevertheless, several approaches have been developed to model the turbulent fluxes (see, e.g., in [15,33,35]). There are many other factors that influence surface energy balance such as clouds, cast shadows, precipitation, and surface moisture content may be difficult to model; therefore, new approaches that combines field data and remote sensing analysis are needed to better address these issues.

Ablation Dynamics
Several approaches have been developed for estimating ablation rates, such as empirical regression and degree-day modeling from weather station data, and surface energy balance modeling. The regression and degree-day approaches heavily rely on the availability of local data, and therefore are highly site-specific [15,80]. Consequently, more recent studies use the surface energy balance approach for modeling ablation rate. For example, the models developed by Nicholson and Benn [15], Reid and Brock [16] are capable of predicting the internal temperature of debris loads at any location and provide satisfactory ablation rate estimates. Based on the energy balance at the ice/air interface, the melt rate at a bare-ice surface under temperate conditions can be calculated as [15,30] where M i denotes the ablation rate of ice, L f is the latent heat of fusion for ice, and ρ i is the density of ice. The energy balance at the debris/ice interface can be written as [67] where Q m is the heat flux used for sub-debris ice ablation, Q ↓ c is the conductive heat flux from the debris, and Q c is the conductive heat flux towards the ice that is not used for ablation. As most models assume the ice is at its melting point, the second conductive heat term is negligible [37].
Reid and Brock [16] discussed the nonlinear nature of temperature gradient within debris, and developed a model to calculate the internal debris temperature profiles: where T d is the debris internal temperature, which is a function of depth (z) and time (t). ρ d , c d , and k d are the density, specific heat capacity, and thermal conductivity of the debris, respectively. This equation can be solved using Crank-Nicholson scheme, in which the debris cover is broken down to multiple layers to support the computation [16]. Most other models, however, assume constant heat storage and a linear debris-temperature gradient, then Q m during the ablation season can be computed using the one-dimensional heat flux equation described by Fourier's law [15,67]: where T i is the ice temperature, h d is debris thickness, and R is the thermal resistance of the debris layer, which can be estimated from thermal imagery or field data. Surface temperature T s can be measured in the field or through remote sensing, or computed for debris-covered areas at each iteration by solving the surface energy balance equation, and set to melting temperature for bare-ice surfaces. Then, the sub-debris melt rate (M s ) can be computed as [15,67] Nicholson and Benn [15] compared the modeled ablation rates as a function of debris thickness using this approach and field measurements on several glaciers, including the Barpu Glacier [81] and the Rakhiot Glacier [58] in the Himalaya, and the results showed a reasonable agreement. Given the importance of including debris properties in the model, studies have been incorporating remote sensing analysis into modeling to address the distribution of debris thickness [46], and debris thermal resistance [30].
Several recent studies have modeled surface ablation rate on debris-covered glacier based on the same surface energy balance, although some energy flux terms may be calculated differently. Fyffe et al. [68] used a distributed model to compute melt rate on the debris-covered Miage Glacier in the Alps. The energy balance model is not novel, but the distributed modeling highlighted the importance of differentiating clean ice, snow and discontinuous thin debris covered "dirty ice" in surface melt simulations. The results are evaluated using runoff data and ablation stake measurements and the modeled values are in a reasonable range. A detailed parameter sensitivity analysis was conducted, and they found that melt rates are most sensitive to air temperature where debris is thin, as debris gets thicker, the sensitivity decreases as the daily cycle of energy flux is attenuated. The sensitivity of melt rate to albedo was found to be moderate likely due to the high spatiotemporal variability of surface albedo. They also discussed the importance of accounting for realistic wind speed, topographic shielding (sky-view factor), the distribution of dirty ice, and debris thermal properties variations in an energy balance-based surface ablation model. Rounce et al. [41] discussed the significant influences of debris thermal conductivity, albedo, and surface roughness on the surface energy balance of Imja-Lhotse Shar Glacier in Nepal. Sensitivity analysis for these three factors showed that surface ablation is most sensitive to changes in thermal conductivity, moderately sensitive to albedo variation and least sensitive to surface roughness. They also suggest that the water content in the debris plays an important role in regulating energy balance and ablation rates. Minora et al. [82] estimated the ablation rates for debris-covered glaciers over the Central Karakoram using a distributed model, the modeled ablation showed a strong agreement with field measurements on the Baltoro Glacier. The corresponding sensitivity tests suggest that an increase in summer air temperature will significantly increase the ablation, the temporal variations in debris thickness and surface albedo also regulate surface melt.
Vertical thinning dominates the mass loss of debris-covered glaciers [60], and the heterogeneous distribution of debris thickness contributes to the high spatial variability of surface ablation [9,25,33]. Figure 3 is an example of surface ablation modeling that depicts the simulated surface ice-mass loss for a hypothetical bare-ice scenario (Equation (8)) of the Baltoro Glacier versus the realistic debris-covered condition (Equations (11)-(12)), following Mihalcea et al. [46]. The amount of surface ice-mass loss due to ablation is represented in meter water equivalent (m w.e.). This comparison highlights the differential ice loss due to the spatially heterogeneous debris thickness distribution. Note that the overall ablation for bare-ice scenario is higher, but the spatial heterogeneity in surface ablation is much higher on the debris-covered glacier. Other processes on debris-covered glaciers, such as supraglacial ponding may further elevate this differential surface ablation pattern.
Given high melt rates in the areas of clean ice or thin debris cover, more sediments are exposed. Therefore, surface debris concentration is likely to increase over time [10,61,83]. Precipitation including snow and rain also contributes to the total surface ablation [8,68]. Models developed by Reid and Brock [16] and Fujita et al. [8] provide a first-order solution to estimate the melt caused by precipitation. Table 1 provides a comparison between these important models and highlights the contributions made by each of them.
The time step (or temporal resolution) should also be considered in modeling, as the assumption of linear temperature gradient within the debris may only be valid for time steps greater than a day [15]. At a hourly or half-hour time step, the daily cycle in energy flux and the nonlinear temperature gradient in the debris layer regulate the surface energy balance [41,68]. Table 1. Summary and comparison of the reviewed numerical models for debris-covered glaciers.

Model Reviewed Model Type Model Highlights
Nakawo and Young [67] Surface Ablation Classic model accounts for debris temperature, thickness, and thermal conductivity. Nicholson and Benn [15] Surface Ablation Improved debris thermal conductivity calculation affected by pore-space and moisture. Reid and Brock [16] Surface Ablation Supports the calculation of thermal conductivity and internal temperature profiles within debris. Fyffe et al. [68] Surface Ablation Distributed model highlights differential ablation under heterogeneous surface cover conditions. Rounce et al. [41] Surface Ablation Sensitivity tests on debris thermal conductivity, albedo, and surface roughness. Zhang et al. [30] Surface Ablation Estimation of debris surface thermal resistance using remote sensing approaches. Fujita et al. [8] Surface Ablation Accounts for ablation due to precipitation, more accurate estimates of wind speed and albedo. Minora et al. [82] Surface Ablation Sensitivity tests on ablation-air temperature relationship. Collier et al. [84] Glacier-Atmosphere Accounts for glacier-atmosphere interactive coupling.
Anderson and Anderson [10] Debris Transport 2D supraglacial and englacial debris advection governed by ice dynamics. Wirbel et al. [56] Debris Transport 3D englacial debris advection and diffusion governed by ice dynamics. Sakai et al. [25] Supraglacial Pond One of the first energy and mass balance model for supraglacial ponds. Miles et al. [38] Supraglacial Pond Improved version of Sakai et al. [25]'s model, more accurate computation of energy and mass fluxes. Steiner et al. [35] Ice-cliff Energy balance on ice-cliffs affected by complex terrain. Buri et al. [36] Ice-cliff 3D ice-cliff evolution model accounting for water-caused ablation at cliff base.

Surface Ice-Flow
The ice-flow velocity field is the basis for modeling glacier ice movement and glacial debris transport [3,9,10,56]. Ice-flow dynamics also governs the geometry changes (length, width, and thickness) of a glacier, therefore playing a fundamental role in prognostic glacier simulation [21]. The flow velocity of a glacier is predominately governed by ice thickness, and many modelers use the shallow ice approximation to estimate the flow velocity (see, e.g., in [10,[85][86][87]). Modeled ice-flow velocities for valley glaciers, however, usually suffer from high uncertainty which could be caused by the choice of parameters, the unknown basal sliding, and topographic conditions [88]. Given the scope of this review, we only focus on the ice velocity at the glacier surface, which governs the advection of supraglacial debris (Equation (1)) and long-term variations in glacier surface topography. Fortunately, surface ice-flow velocity can be derived from satellite images with relatively high accuracy.
Studies over High Mountain Asia have suggested that the ice-flow can still be very active on a heavily debris-covered glacier even though its terminus is stable (not retreating or advancing) (see, e.g., in [49,89,90]). Most studies use feature tracking based remote sensing analysis to estimate surface ice-flow velocity. Dehecq et al. [91] investigated mountain glacier velocities from a complete satellite archive based on a robust processing strategy. The annual velocity fields were produced for over 76,000 kms 2 of glacierized area over the Pamir-Karakoram-Himalaya. Dehecq et al. [92] delivered glacier surface velocity fields by applying feature tracking to Landsat-7 image pairs using JPL auto-RIFT program [93]. Wendleder et al. [94] presented glacier surface velocities at various time series using intensity offset tracking applied to the Synthetic Aperture Radar (SAR) data from 1992 to 2017.
Rankl et al. [95] derived surface velocity mosaic of Karakoram from TerraSAR-X, ALOS PALSAR, and Envisat ASAR image data acquired between 2007 and 2011. Figure 4 shows the estimated surface ice-flow velocity field of the debris-covered Hispar Glacier using feature tracking technique on Landsat-7 image pair.

Supraglacial Water Bodies
Supraglacial water bodies are commonly found on debris-covered glaciers in the High Mountain Asia [5,26,97]. Supraglacial lakes, ponds, and adjacent ice-cliffs are considered zones of rapid melting on debris-covered glaciers that counter-balances the debris insulation effect [5,11,25,33,35,36,38,44,47,65,98]. Supraglacial water bodies also play an important role in the storage and transport of meltwater in a glacier system [21], which are key processes that alter surface morphology and glacier thinning [27,55,97]. Surface meltwater is usually transported via stream channels during the ablation season and can also be temporarily stored on the glacier surface; therefore, many supraglacial ponds and lakes have a very seasonal behavior [27,31,37]. Supraglacial and englacial fluvial systems can also regulate glacier mass balance by expanding water channels, altering debris thickness distribution, and transporting the heat stored in the debris.
Most existing debris-covered glacier models, however, neglect the surface processes that are related to supraglacial water bodies due to numerical complexity and the lack of data. This is a major limitation that lead to an underestimation of ablation on debriscovered glaciers [5], which could be one of the reasons why most models cannot explain the increasing number of supraglacial water bodies and accelerated downwasting observed on large glaciers in the Himalaya [8,22].
This section reviews recent studies on supraglacial ponds/lakes, ice-cliffs, and surface meltwater drainage with an emphasis on the modeling efforts.

Supraglacial Ponds, Lakes, and Ice Cliffs
Supraglacial ponds and lakes usually form along englacial conduits and in the lowermid ablation zone where glacier-profile slopes are relatively low and topographic depressions are present [33,44,99,100]. Large lakes on Himalayan glaciers can be a few kilometers across [26], while ponds can be as small as a few meters cross [38]. A typical lifespan for a supraglacial lake on Himalayan glaciers is a few years, during which they can form, grow, merge, and be completely drained when they intersect with englacial conduits, and many smaller ponds only exist for one ablation season [27,28,55,101]. Many supraglacial lakes are hydrologically connected [27,28,102], and large lakes usually have higher connectivity with other lakes, and most of them exhibit periodic filling and drainage processes [26]. Most lakes also contribute to a significant amount of englacial ablation by warm water outflow [3,25,38,65]. Large supraglacial lakes can also be destructive, as glacier-lake outburst floods originate from DCGs in the Himalayas, and have caused injuries, deaths and property damage to downstream villages [99,103]. Figure 5 shows a satellite imagery over the lower ablation zone of the Khumbu glacier, on which a large number of supraglacial water bodies can be identified.
Supraglacial ponds and lakes exhibit a lower albedo than surrounding areas [104], which allows them to absorb more solar energy, rapidly melting the ice and expanding the lake area. Previous studies have estimated that the ablation rates around these water bodies can be much higher than that of most debris-covered areas; therefore, supraglacial ponds and lakes are significant contributors to the total ablation [5,25,98]. Furthermore, a dense spatial distribution of lakes often increases glacier thinning because they form a large number of englacial channels, and the collapse of the englacial channel roofs creates new lakes that further accelerate ablation and thinning.
A numerical model for estimating the melting caused by supraglacial pond was developed by Sakai et al. [25] based on their field study on the Lirung Glacier in the Langtang Valley, Nepal. The core of this model is the energy balance of a supraglacial pond that accounts for heat flux input and output due to metlwater, heat storage of the lake, and the bare-ice versus debris covered areas beneath the lake surface: where Q is the net heat flux (input) as discussed in the surface energy balance section; Q in is the input heat flux from meltwater inflow into the lake, which in most cases, can be neglected [25]; Q out is the output heat flux from the water outflow from the lake, which dominates the englacial ablation such as the expansion of conduits [25]; ∆Q t is the change in heat storage of the lake; Q i is the heat flux for ice melting at the underwater ice-cliff; and Q d is the heat flux into the debris layer under the water surface. Miles et al. [38] proposed an improved model of mass and energy balance for supraglac ial ponds following the concept of Sakai et al. [25]. This model uses more accurate approaches for modeling energy flux components, and field measurements of a small pond on Lirung Glacier, Nepal were used to validate the model results. They also compared the results using several different existing approaches by Sakai et al. [25], Lüthje et al. [104], and Steiner et al. [35], and showed a good agreement. Combined with field measurements, this physically based model provides a reasonable estimate of the amount of ice loss due to supraglacial ponds.
The first catchment-scale modeling of the melt enhancement of ponds was presented by Miles et al. [5], in which four glaciers in the Langtang Valley, Nepal in 2014 were modeled using field meteorological data and satellite remote sensing analysis. Results indicate that ponds cause surface melt at a rate of 0.20 ± 0.03 m/year, which is significantly higher than the debris-covered area. These ponds only cover 1.0% of debris-covered area or 0.3% of total glacier area, but are responsible for about 1/8 of total ice loss in the catchment.
The development of supraglacial ponds and lakes on a debris-covered glacier is not only governed by meltwater drainage and filling, but is also strongly controlled by surface topography [99,104], and surface cover conditions (bare-ice or debris) [105]. Studies have suggested that the morphological changes on a glacier surface contributes to the increasing number of supraglacial ponds and lakes [3,106], for example, the undulating surface and the gentle slope of the ablation zone encourages the formation of lakes in depression areas [25][26][27][28]105], and the expansion of lakes further lowers the slope and accumulates more meltwater.
Debris transport around ponds and lakes also plays an important role in their expansion. Studies have identified significantly reduced debris thickness due to lakes development [45], and relatively rapid sediment flow down ice-cliffs during their backwasting [44,53], which explains why the water-facing slopes usually have a thinner debris cover that enhances melting, creating steeper slopes and rapid backwasting of ice-cliff [33,44,65]. Unfortunately, the effects of topography and debris transport have not been adequately characterized in existing models; therefore, we do not have a good understanding of multi-scale surface morphological effects on supraglacial pond/lake development.
Supraglacial water bodies will most likely show an accelerated growth on Himalayan glaciers given current atmospheric temperature trends [5,105]. For example, Gibson et al. [29] estimated that the number of water bodies on the Baltoro Glacier has increased from 234 in 2001 to 570 in 2012, and the area has expanded from 0.66 km 2 to 2.04 km 2 accordingly, which suggests that supraglacial water bodies can expand quickly on debris-covered glaciers and will have a more significant impact on the glacier in the future.
Ice-cliffs usually form along with the expansion of supraglacial ponds and lakes [33,36,44,65]. Sakai et al. [44] classified the ice cliffs on the Lirung Glacier into four types: decayed, temporary, developed, and stable. They also found that the orientation and inclination of ice-cliffs poses a significant control on the ice-cliff ablation because ice-cliffs with certain orientation and inclination angles receive more energy given the solar geometry and cast shadows, and the inclination angle also controls the thickness distribution of debris cover on the ice-cliffs which is another factor that affects the ablation rate. Our current understanding of ice-cliff formation and evolution on debris-covered glaciers is far from sufficient, as they are not only controlled by surface processes, but are also controlled by englacial processes, for example, studies have attributed the formation of some ice-cliffs to the collapses of englacial conduit or the steepening of ice slope due to water undercut [25,44].
Numerical modeling can provide valuable insights into ice-cliff dynamics on debriscovered glaciers. For example, Steiner et al. [35] and Buri et al. [36] used energy balance models for investigating ice-cliff retreat on Lirung Glacier based on a high-resolution DEM. Their models account for the spatial patterns of solar energy fluxes in the complex terrain around the cliff; they found that the topography has a non-negligible effect on ice-cliff dynamics due to the local shading and the adjacent terrain irradiance. Their modeling results confirmed that ice-cliffs can make significant contribution to the total ice mass loss on a debris-covered glacier. The first 3-D ice-cliff evolution model was developed by Buri et al. [36], in which additional processes including water-caused ablation at the cliff base and the reburial by surrounding debris are accounted for. The simulations results show a reasonable agreement with observations and demonstrated the capability of predicting icecliff evolution over an ablation season using numerical modeling. These models, however, have not addressed the pond-cliff coupling [36] and the processes on the ice-cliff surface, such as the redistribution of debris, microclimate, surface roughness, and the potential lowering of albedo due to debris thickness variation and meltwater content. Improved models are needed to address these issues.
The filling and draining cycles also control the evolution of supraglacial ponds and lakes [27,28,105], but they are extremely difficult to investigate as that would require knowledge of subsurface (englacial) structure and hydrological conditions, which cannot be observed by common remote sensing techniques. Modeling water temperatures can also be difficult, which is a challenge for estimating the amount of ice loss that is related to supraglacial water bodies. Collectively, accurate modeling of supraglacial ponds, lakes, and ice-cliff evolution is challenging, and it requires a large amount of field data (e.g., high-resolution surface topography, water depth and temperature, energy fluxes, surface albedo and debris thickness) to facilitate model development and validation.

Supraglacial Drainage and Ponding
Supraglacial ponding which accounts for meltwater filling and drainage can be numerically modeled given the estimates of distributed melt water production and highresolution DEMs of the glacier surface. Basic approaches for flow-path determination have been successfully applied on the Lirung Glacier in the Himalaya to study the connectivity of supraglacial ponds [27]. Lüthje et al. [104] developed a model to investigate supraglacial pond filling and draining on the Greenland ice-sheet margin based on the simulation of meltwater production coupled with topography-controlled surface drainage dynamics: where h w is water depth above the ice surface, He(h) is a Heaviside function to prevent negative water depth, ρ i is the density of ice, ρ w is the density of water, D is a parameter that addresses the dynamic viscosity of water and horizontal permeability of the substance (ice or sediment) in which meltwater travels through [107], and z i is the ice-surface elevation. Field studies suggest that meltwater can travel long distances in active marginal and sub-marginal channels [28]. These channels may be bounded between ice and valley walls, as the surface elevation of glacier margins is typically lower than at the center of the glacier, and meltwater can accumulate in these channels while flowing towards the terminus. The amount of runoff, however, does not represent the total meltwater production due to supraglacial and englacial water storage [37].
The travel time for meltwater to drain out of the glacier ranges from days to hours, depending on the location where the meltwater is generated [37,108]. The time lag is usually smaller in the ablation season when the discharge reaches its maximum in late July or early August [37]. This time delay is due to multiple factors, such as the amount of snow, firn, and debris along the path, as well as the time for subglacial and englacial conduits to develop. Otherwise, the travel speed on ice is similar to the speed in open channels [37].

Feedback and System Couplings
Surface processes on a debris-covered glacier are complicated by a variety of feedback mechanisms and subsystem interactions [11,20,39]. A better understanding of these couplings and feedbacks is important for more accurate assessments of a debris-covered glacier's response to climate change. Many processes operate at very different spatiotemporal scales, making them difficult to investigate in the field. Therefore, numerical modeling plays a significant role in exploring feedback and system couplings.
Studies have identified several major components and their couplings that govern the surface processes on a debris-covered glacier. Here, we summarize them into a conceptual diagram (Figure 6), in which the cause-and-effect relationships and feedback mechanisms are highlighted. In general, climate controls the overall "health" of a glacier, which is typically represented as the mass balance. The mass loss and meltwater production govern the development of supraglacial water bodies and the flow velocity, which alters surface topography and debris distribution. The conditions of debris load, supraglacial water bodies, and surface topography also influence ice mass loss by regulating the surface energy balance as discussed in previous sections. This section reviews the findings about system interactions on debris-covered glaciers with a focus on the feedback mechanisms that are important to assessing climate-glacier sensitivity. Figure 6. Major system couplings and and feedbacks (dashed lines) that govern surface processes on a debris-covered glacier. Positive feedbacks, such as the one between surface mass balance, supraglacial water bodies, and topographic condition, play a critical role in the system as they may be responsible for the high-magnitude ice loss observed on many debris-covered glaciers.
Surface morphological conditions on a debris-covered glacier is dynamic, which influences multiple supraglacial processes. One important controlling factor is topography, studies have indicated that topography has non-negligible effects on glacier surface ablation because slope, slope azimuth, and basin topographic shielding directly control the total amount of shortwave radiation received by the glacier surface which drives ice loss [109,110]. Topography also affects the distribution and transport of supraglacial debris. For example, Rounce et al. [45] found a strong relationship between the decrease in debris thickness and the development of supraglacial lakes and ice cliffs. Collectively, topographic conditions have both a direct and indirect impact on glacier surface ablation. Therefore, a dynamics surface topography should be accounted for in surface ablation modeling. Surface topography also plays an important role in surface ponding [44,99] where supraglacial water bodies are more abundant in zones that exhibit relatively low gradients [33,44,99,100]. Topographic depressions are ubiquitous on debris-covered glaciers due to high spatial variability in ice topography, debris thickness, and ablation rate [30,44,46]. Topographic depressions caused by differential ablation and surface water flow may also contribute to ponding, as they provide topographic sinks for meltwater to accumulate. Therefore, a positive feedback exists on the glacier surface: supraglacial ponds cause spatially heterogeneous surface melt on debris-covered glaciers [24,38,97,105], and this differential downwasting pattern, in turn, facilitates new pond formation, which further enhances glacier ice loss [3,5].
Surface albedo is another factor involved in feedbacks. Field studies have shown that surface albedo exhibits significant spatio-temporal variability on a debris-covered glacier [111,112]. Positive feedback exist between the lowing of surface albedo and enhanced melting: the presence of meltwater lowers surface albedo, which encourages energy absorption and elevates melt rate that generates more meltwater. Pritchard et al. [113] showed that this mechanism has significant influence on millennial-scale mass loss of ice sheets. For certain locations on the glacier surface, such as on ice-cliffs, the decrease in albedo is more significant due to the thinner debris and high moisture content, which is also coupled with the topographic factors (such as the orientation of the cliff wall and the high adjacent terrain irradiance on the cliffs) to form positive feedback that accelerates melting [25,33,53].
Supraglacial ponds play an important role in system coupling and feedbacks on debris-covered glaciers. Multiple studies [3,5,24,33,53] have described a potential positive feedback between surface melting and supraglacial pond formation: the heterogeneous spatial distribution of supraglacial water bodies and associated ice-cliffs enhances heterogeneous surface downwasting, which, in turn, facilities the expansion of existing ponds and the formation of new ponds and ice-cliffs in a positive feedback fashion. Specifically, this positive feedback is controlled by several factors and processes, such as topography, albedo, debris thickness, englacial melting, and calving. As a pond grows, its surface area and surrounding ice-cliff area become larger, and therefore they absorb more solar energy due to the lower albedo, which enhances ablation and leads to further expansion of the pond. Debris thickness distribution is also regulating the development supraglacial ponds in terms of their sizes and locations [41]. Rounce et al. [45] found that debris usually gets thinner as water bodies and ice-cliffs expand, as a result, the melt enhancing effect of thin debris will become significant, which accelerates pond expansion and further expands the thin debris-covered area [33,34]. In addition, the temperature of water bodies and ice-cliffs is affected by the elevated adjacent terrain irradiance due to debris cover [26], and the ablation caused by outflow of high-temperature water from ponds and lakes can be significant [25]. These positive feedbacks may partially explain the observed increase in supraglacial water bodies on some debris-covered glaciers over High Mountain Asia.
The interaction between meltwater and debris is also responsible for chemical properties variations in glacial lakes and ponds [114][115][116]. Based on data collected from 20 glacier-fed lakes located in central southern Himalaya, Salerno et al. [115] found that glacier melting was the main cause of the solute (major ion concentrations and sulfate) increase in these lakes partially due to the enhanced debris/water interaction given the changing climate.
In addition, Collier et al. [19] identified the feedback between atmosphere and glacier climatic mass balance (CMB), and by including the related couplings in their model, they improved the simulation results of surface temperature and snow albedo. Supraglacial conditions influence atmospheric conditions as debris load locally modulates the near surface microclimate due to the sensible and latent heat transfer between the glacier surface and the atmosphere, such that surface heating can be accentuated, possibly initiating atmospheric convective motion [18,19]. Collectively, these interactions and especially the positive feedbacks may result in debris-covered glaciers being more sensitive to climate change than previously thought.

Model Inputs and Validation
Most glacier models are developed based upon realistic initial conditions and are then validated with appropriate data. Therefore, data collection is an essential component of most modeling works. The most valuable data come from field work, because accurate information about key driving parameters such as air temperature, radiation fluxes, and debris thickness can only be measured in the field. In most cases, however, field data are insufficient and limited in time and space. With the advancements in sensor technologies and remote sensing approaches, many researchers now use satellite remote sensing data to support modeling. Unlike the point-based in situ measurements, remote sensing analysis enables rapid mapping over an entire glacier surface and provides valuable information about glacier boundaries, surface cover (e.g., debris, water bodies, ice, and snow), surface topography, surface albedo, and surface temperature, which are all important input parameters in glacier simulations. Another type of data used in glacier modeling derives from numerical simulations that were calibrated with observations, such as the atmospheric reanalysis products. In this section, we review the main data sources for glacier modeling and discuss the remote sensing approaches for extracting information about debris-covered glaciers.

Meteorological Data and Ablation Rate
Meteorological variables make a major component of in situ measurements. For example, Mihalcea et al. [66] reported the use of automatic weather stations (AWSs) for collecting meteorological data on the Baltoro Glacier. The AWSs provided continuous measurements for monitoring supraglacial energy transfer and near-surface climate; the parameters they recorded include air temperature, wind speed and direction, air pressure, humidity, precipitation, and long-and shortwave radiation, which have been proven to be critical for investigating the energy exchange at the glacier surface [46,66]. Meteorological variables are usually recorded at different elevations and at different times to provide information about the spatiotemporal variabilities. The distance between the sensors and the glacier surface need to be measured, as the height of meteorological measurements is required for estimating sensible and the latent heat fluxes.
Ablation rate is usually measured by stakes drilled into the ice at multiple locations during the ablation season. For debris-covered glaciers, ablation rates are often measured at locations with different debris thickness (see, e.g., in [14,41,58]). A major limitation of the ablation stake approach is the noncontinuous reading, some researchers have explored the use of automatic ablation sensor to provide continuous ablation rate measurements [66,117].
Reanalysis climate datasets are becoming more widely available [118], which provide another option for modelers. Several studies have incorporated reanalysis data in their models, for example, Collier et al. [84] obtained the incoming longwave radiation estimates from the ERA-Interim reanalysis and used them in their models jointly with ASW data to understand the coupling between glacier surface and the climate. Fujita et al. [8] utilized the NCEP/NCAR gridded reanalysis data to simulate meltwater production in the Tsho Rolpa Glacial Lake-Trambau Glacier basin in Nepal, and the results show a reasonable agreement with hydrologic records. Researchers also noted that the accuracy of reanalysis products does not meet the requirement of some models if they are not calibrated with in situ measurements [8].

Remote Sensing of Glacier Surface
Remote sensing of glacial surfaces has advanced significantly over the past two decades [119][120][121][122][123]. Classification is a common technique for mapping glaciers, where a predefined semantic label is assigned to each pixel. Multiple classification algorithms have been applied on satellite data for glacier and ice sheet mapping including supervised and unsupervised models, fuzzy classification techniques, and band math and threshold approaches [124,125]. However, accurate mapping of debris-covered glaciers is still a challenging topic because debris loads exhibit similar spectral characteristics as rocks and sediment of the surrounding terrain [119,[126][127][128]. To address the issue, more sophisticated approaches have been developed. For example, Bishop et al. [119] proposed a geomorphometrics-based framework to analyze the topography of debris-covered glaciers where DEM created from two SPOT panchromatic stereopairs was utilized as input. Paul et al. [126] utilized Landsat TM4/TM5 ratio image and slope gradient simultaneously to map debris-covered glaciers. Neighborhood analysis and change detection was applied as well for further improvement of the resulting glacier map. Multiple supervised classification algorithms were compared by Brenning [127], who also investigated the combination of terrain attributes derived from Shuttle Radar Topography Mission (SRTM) DEM and Landsat ETM+ data for mapping rock glaciers in the San Juan Mountains, Colorado.
Supraglacial debris mapping is of great importance to studying debris-covered glaciers. Several approaches have explored to map the spatial distribution of debris loads. For example, Scherler et al. [129] used the TM4/TM5 ratio images of Landsat TM data to discriminate between debris-covered and debris-free ice, and identified the coupling between hillslope and glacier processes in steep terrain. Gibson et al. [29] manually mapped the distribution of debris and supraglacial water bodies on the Baltoro glacier by analyzing multiple Landsat and ASTER satellite imageries, and the geomorphological maps provided insights into the nature of debris loads and the response of debris-covered glaciers to climatic change.
Debris thickness and thermal properties can also be estimated from remote sensing products. For example, Mihalcea et al. [46] found an regression equation between measured debris thickness and glacier surface temperature acquired from the ASTER surface kinetic temperature data over the Baltoro Glacier; Figure 7B shows the estimated debris thickness based on this approach. Juen et al. [40] used three different regression approaches to estimate debris thickness for Koxkar Glacier from ASTER land surface temperature data. Similar to the regression solution applied to the Baltoro Glacier [46], they also found that the exponential regression is more accurate for estimating debris thickness from remotely sensed surface temperature by comparing with field radar measurements. However, both approaches suffer from high uncertainty for areas covered by thin debris. Similar approaches allow a fast mapping of debris-depth conditions over a large area, but the limitations are obvious given the empirical and site-specific nature of the determined relationship: the surface temperatures from satellite data typically have a resolution of over 30 m; therefore, it is almost impossible to register that value to a specific location on the glacier surface where the debris thickness is measured. These empirical approaches also do not account for the variations in albedo, particle size, topographic shadings, and moisture content in the debris layer that could compromise the accuracy. Thermal properties of the debris (e.g., the thermal resistance, which is defined as debris thickness divided by its thermal conductivity) are often needed in modeling works, and they can also be derived from remote sensing analysis. For example, Zhang et al. [30] developed a method to estimate debris thermal resistance based on the albedo computed from visible/near-infrared of ASTER data, the brightness temperature computed from the thermal band, and the reanalysis climate data. Although the above-mentioned issues still exist, such remote sensing approaches provide a first-order approximation of the debris properties that can be used as model inputs.
Remote sensing also plays a vital role in mapping supraglacial water bodies [130][131][132]. Fujita et al. [130] computed normalized difference water index (NDWI) from satellite data, and they successfully delineated lake boundaries on Himalayan glaciers by thresholding the NDWI images. Gibson et al. [132] also used NDWI to assist manual mapping of supraglacial water bodies on the Baltoro Glacier. Miles et al. [27] identified hydrological sinks using DEM to understand the spatial pattern of supraglacial ponding and the development of conduits. Figure 7C depicts the extracted water bodies from an ASTER multi-spectral imagery over the tongue of Baltoro Glacier ( Figure 7A) using the NDWI approach.  [46]. (C) Supraglacial water bodies extracted from the satellite imagery using the normalized difference water index approach by Fujita et al. [130].
With the advances in artificial intelligence, machine learning schemes are now being introduced in glacier mapping. Inspired by the semantic segmentation, Xie et al. [133] proposed a convolutional neural network (CNN)-based semantic segmentation framework to extract debris-covered areas. Landsat-8 images, geomorphometric parameters, and DEMs were selected as the input datasets. Even though limitations exist such as the cost of manual labeling, machine learning shows great application potential in mapping supraglacial features.
Surface topography is one of the most important model inputs for simulating supraglac ial processes. Digital elevation models are widely used to represent the topography in a grid or raster data structure. All major global DEM datasets are acquired by satellite remote sensing via Radar or photogrammetry. Popular global DEMs including the SRTM (Shuttle Radar Topography Mission) and the ASTER GDEM (Advanced Spaceborne Thermal Emission and Reflection Radiometer Global DEM) have a spatial resolution of 30 m, which meet the requirements of most large-scale glacier models. Some modeling efforts (such as supraglacial pond or ice-cliff models), however, require higher resolution DEM that must be generated using data collected in the field. For example, Buri et al. [36] used 0.2 m resolution DEMs generated from UAV photographs to model the backwasting of ice-cliffs on the debris-covered Lirung Glacier. DEM differencing is also very useful, as surface elevation changes can be detected by comparing multi-temporal DEM datasets that are spatially well-registered. Several studies [27,100,122,134] have used DEM differencing to study ice mass loss, the elevation of glaciers' equilibrium line altitude (ELA), and the evolution of supraglacial water bodies; all of these information are very valuable for model validation.
Ice-flow velocity field at the glacier surface is another vital model input that is typically derived from remote sensing analysis (Figure 4). Feature tracking using spectral imagery provides a solution for ice-velocity derivation by measuring the displacements based on the location of identical surface features from sequential images (image pairs) [135][136][137][138][139][140]. Surface features are represented by the patterns of a group of individual pixels. By shifting small search windows across each single band image pair, the displacement of the dominant feature within the window is computed through the normalized covariance correlation method [141]. Although image matching is not limited by topography, a prerequisite of enough cumulative displacements, appropriate time intervals between images should be considered in data selection in order to obtain high-accuracy ice velocity estimates. Satellite microwave remote sensing, particularly Synthetic Aperture Radar (SAR), provides advantages of high accuracy flow-rate estimation and all-weather monitoring [142]. Two SAR images of the same surface area acquired at slightly different orbit configurations are usually combined to exploit the interferometric phase difference occurring between the acquisition time intervals.

Issues and Research Directions
Numerical modeling has greatly facilitated the study of surface processes on debriscovered glaciers. However, given the difficult to accurately characterize multiple factors (and their interactions) that operate at different scales, more sophisticated models are always warranted to better address supraglacial features, processes and feedback mechanisms that are oversimplified or neglected in existing models. From an implementation perspective, new numerical schemes and computational optimization strategies are also needed to support more complex numerical simulations. We identified the following issues based on the review of existing works, and they point to the directions for future studies.
Supraglacial debris: (1) Debris thickness distribution is one of the most important controlling factors in modeling debris-covered glaciers. Most existing approaches for estimating debris thickness (see, e.g., in [30,40,46]), however, suffer from high uncertainty and are more or less site-specific. (2) The input and transport rates of debris are largely unknown, models are need to characterize and differentiate debris load contribution from different sources, such as headwall erosion, hillslope process, and englacial transport. Hypothetical scenarios have been investigated [10], but more realistic modelings are needed.
(3) The pathways and rates of glacial debris discharge from the ablation area are not well constrained. Current models cannot explain the full sediment balance at the glacier terminus [54]. (4) The topographic and hydrological controls on debris transport is largely unknown [143]. (5) The mixing of debris with different lithology, grain size, and moisture content has not been accounted for existing models. (6) The melt-enhancing effect of very thin debris layer or "dirty ice" has not been adequately explored. (7) Improved models also need to account for the non-turbulent debris flux and fluvial debris flux [21].
Surface ablation dynamics: (1) The ablation dynamics of ice-cliffs and supraglacial ponds/lakes are oversimplified or neglected in many existing simulations [5,38,44]. (2) The influence of clouds on the longwave radiation has not been fully addressed [74]. (3) The nonlinear profile of thermal properties within the debris layer has not been accounted for. (4) Surface albedo distribution requires more accurate and continuous estimates rather than being restricted by instantaneous remote sensing data.
Supraglacial water bodies and ice-cliffs: (1) The filling and drainage processes that control the evolution of many supraglacial ponds and drainage efficiency [27,28,55,105] have not been adequately modeled for debris-covered glaciers. (2) Existing models do not account for the formation of supraglacial ponds due to the collapse of water channel roofs [25]. (3) Ice-flow control on the long-term evolution of supraglacial ponds and icecliffs if often neglected. (4) The adjacent terrain irradiance on the energy balance of ice-cliffs and the melt caused by the outflow of high-temperature water need to be accounted for in improved models [25,26,38,44]. (5) The relationship between the decrease in debris thickness and the development of supraglacial ponds and ice-cliffs has been identified [45], but the controlling processes are still largely unknown.
A variety of feedback mechanisms and system couplings complicate a debris-covered glacier's response to climate change [4,9,11]. Therefore, an integration of existing models that focuses on different components of the glacier system (e.g., surface energy balance and ablation, debris transport, supraglacial ponding, and ice-cliff backwasting) is required in future studies in order to better understand the system couplings and associated time scales, and to provide a more accurate assessment of glacier sensitivity to climate change [4]. The realistic distributions of debris thickness, surface albedo, lithological components, supraglacial ponds, ice-cliffs, and ice-flow velocity are necessary for improved models. Many glacier surface features and properties can be estimated using remote sensing analysis; therefore, glacier surface characterization using high-resolution remote sensing data can also be very useful for validation of simulation results, which represents another research direction.

Conclusions
Debris-covered glaciers in High Mountain Asia make up a major component of the regional water cycle that governs water resources for millions of people. The mass balance conditions of these glaciers are significantly governed by their surface processes that are related to supraglacial debris, water bodies, ice-cliffs, and topographic conditions. These features and processes complicate glacier response to climate change and contribute to the high degree of uncertainty in our current assessments of glacier-climate sensitivity in the region.
Recent studies have demonstrated the capability of numerical modeling in understanding complex supraglacial processes and showed that simulations can help explain conflicting views about these glaciers, such as the "debris-covered glacier anomaly", and the rapid heterogeneous downwasting observed on many glaciers. Therefore, this review focused on the numerical modeling efforts in characterizing supraglacial features and processes that govern ice mass loss and morphological changes on debris-covered glaciers. We also discussed the associated feedbacks and system couplings on the glacier surface, and the remote sensing approaches that facilitate modeling. In addition, we identified some limitations in existing models and suggest research directions for future studies. Specifically, this review highlighted the importance of the following components in modeling debris-covered glaciers.

•
The properties and transport of supraglacial debris. • Surface energy balance and ablation dynamics. • Supraglacial water bodies and ice-cliffs. • Surface ice-flow. • Systems couplings and feedbacks between debris load, ablation, water bodies, and topographic conditions.