Understanding Complex Debris-Covered Glaciers: Concepts, Issues, and Research Directions

Understanding the climate-glacier dynamics of debris-covered glaciers is notoriously difficult given a multitude of controlling factors and feedback mechanisms involving climate forcing, debris-load properties, supraglacial water bodies, and multi-scale topographic effects. Recent studies have provided insights into controlling factors, and have reported the presence of anomalies that contradict the general consensus of the protective influence of debris loads on ablation dynamics. Nevertheless, numerous processes that regulate glacier dynamics at various spatial and temporal scales have not been adequately accounted for in statistical and numerical modeling studies. Furthermore, important feedbacks involving ablation, topography, irradiance, gravitational debris flux, and supraglacial ponding are often neglected or oversimplified in existing models, which poses a challenge to our understanding of conflicting field observations such as the accelerated mass loss of many Himalayan glaciers, and glacier-subsystem responses (ice-flow, debris flux, surface morphology, and supraglacial water bodies) to climate forcing. This paper provides insights into the complexity of debris-covered glacier systems by addressing concepts and issues associated with forcing factors and glacial processes, and highlights the importance of understanding system couplings and feedbacks. Specifically, we review recent studies on debris-covered glaciers and utilize simulation results based on the Baltoro Glacier in the central Karakoram to discuss important concepts and issues. Our results demonstrate that climate forcing, the properties and transport of debris, topography and supraglacial water bodies are the key controlling factors in a debris-covered glacier system, and that their coupled effects and positive feedbacks may increase the ice loss of a debris-covered glacier. We also recommend new research directions for future studies.


INTRODUCTION
Debris-covered glaciers (DCGs) present a series of unique and complex challenges to those investigating climate-glacier dynamics and their role in downstream impacts such as freshwater availability, hydroelectric power generation, natural hazards and landscape evolution (Bishop et al., 2002;Seong et al., 2009;Immerzeel et al., 2010;Reid and Brock, 2010;Bush et al., 2020). While debris-free glaciers are relatively well-understood and modeled, DCGs introduce a number of additional complex dynamical and thermodynamical feedbacks between the underlying glacier ice, the debris cover, the adjacent-terrain geomorphological conditions, and the overlying atmosphere. The thermal and radiative characteristics of the debris cover itself depend on geographic location and the surrounding lithology. The thickness of the debris cover, which so critically determines whether the debris acts as an insulator for the ice or accelerates ice ablation, depends on the steepness of the surrounding topography, regional erosion characteristics, and glacier dynamics, which in turn depend on atmosphere-glacier interactions that have not yet been fully quantified Dobreva et al., 2017;Farinotti et al., 2020).
Currently, we have not achieved an agreement on the state and fate of many DCGs in High-Mountain Asia due to uncertainties related to climate, topographic evolution, debris-load forcings, and the development of supraglacial water bodies and ice-cliffs. Most researchers agree that debris insulation makes a glacier less sensitive to climate change (e.g., Reid and Brock, 2014;Pratap et al., 2015;Anderson and Anderson, 2016;Vincent et al., 2016). Recent observations and simulations, however, show that many DCGs exhibit comparable ice-loss to debris-free glaciers despite the presence of supraglacial debris Immerzeel et al., 2013;Fujita et al., 2014;Salerno et al., 2017). Therefore, it appears that DCGs may be more sensitive to climate change than previously thought, as a multitude of processes and feedbacks have not been accurately characterized. Existing DCG models do not provide a satisfactory explanation to these conflicting interpretations of climate-glacier dynamics, which is partially due to the over-simplified assumptions of important processes and feedbacks that involve ablation dynamics, debris flux and property variation, surface ponding, and topographic evolution.
Mountains that exhibit DCGs typically have substantial vertical relief and erosional forces that produce significant mass flux from the sidewalls onto the glacier. In the central Karakoram in Pakistan, extreme relief and complex climate systems generate some of the largest DCGs in the world (Bishop et al., 2001;Hewitt, 2005;Copland et al., 2009). Precipitation rates can be highly variable in space and time caused by the dominance of regional climate systems and topography. A particular system such as the Westerlies or the Monsoon could dominate over an area over a specific time frame, and it is also possible that an area could receive precipitation from both systems for a period. Furthermore, topographic variation partially controls the degree of orographic precipitation and this would also vary depending upon which system dominants over an area (i.e., moisture conditions, wind velocity and direction). This very high degree of coupled climate-topographic interactions presents a challenge to the numerical modeling of climate-glacier dynamics.
The surface-energy balance determines the ablation rates on DCGs, and net long-wave and short-wave fluxes play a major role in that balance, but the atmospheric conditions through which that radiation must pass play an equally important role: if the atmosphere is cloudy then the incoming short-and long wave radiations are fundamentally altered. Additionally, atmospheric temperature above the surface is critical to determining precipitation type (rain, snow, or some combination of frozen/liquid water). For example, debris cover can provide a surface heat source that drives convection and precipitation in the overlying atmosphere (e.g., Collier et al., 2015); depending on the temperature, that precipitation may fall as snow which would accumulate on the glacier ice in the absence of the debris, but given the presence of the warm debris the snow may melt and percolate down to the glacier surface, thereby enhancing ablation.
Supraglacial ponding is another critical process on DCGs. Supraglacial ponds and lakes play an important role in glacier ice-loss as well as glacial hydrology (Fountain and Walder, 1998;Sakai et al., 2000;Wessels et al., 2002;Miles et al., 2016Miles et al., , 2018. Researchers have identified an increasing number of ponds on some large DCGs (Gibson et al., 2017). Studies have also indicated that the increase in spatial density of supraglacial ponds can be facilitated by localized thinning, insufficient drainage, collapse of water channels, and overall lowering that decreases the slope of the glacier altitude profile in the ablation zone (Sakai et al., 2000). Unfortunately, many of these important mechanisms have not been adequately studied and are often neglected in existing glacier models and simulations.
Such complexities are only now beginning to be understood through the use of numerical models (e.g., Collier et al., 2013;Rowan et al., 2015;Anderson and Anderson, 2016;Miles et al., 2018). Field observations of en-debris and sub-debris processes can be difficult to obtain given time, economic, and logistic constraints. Nevertheless, field observations are essential for understanding the magnitude of specific properties and processes, providing constraints for input into numerical models, and for validation and accuracy assessment of numerical modeling efforts. However, many mechanisms that govern glacier dynamics can not be fully understood through field observations and remote sensing data. Consequently, numerical modeling that evaluates new parameterization schemes to account for processes and important feedback is sorely needed to better understand complex climate-glacier dynamics of DCGs.
The objectives of this paper are therefore to discuss a multitude of complex interactions and uncertainties that limit our understanding of how DCGs may be responding to climate change. We review our current understanding of debriscovered glaciers and identify gaps in knowledge about some processes and feedbacks that appear to play important roles with respect to ablation dynamics. We also utilize numerical simulations to demonstrate various concepts and processes. We specifically discuss various issues and use our simulation results to identify future research pathways that can help eliminate uncertainty and promote new empirical and numerical-modeling research. We acknowledge that we cannot account for all processes, nor utilize fully comprehensive parameterization schemes that accurately predict process rates (an objective of numerical modeling research). Rather, we attempt to provide for a first-order characterization and understanding of processes and feedbacks based upon the simulations that we have incorporated into our discussion. Furthermore, it should be noted that this work should not be viewed as a formal modeling paper, complete with parameter sensitivity and uncertainty analyses, as this is beyond the scope of our objectives.
We synthesize our treatment of this complex topic by organizing the paper based upon concepts related to climateglacier system components. We first address climate forcing and specifically address radiative and precipitation forcing, where we identify some processes that should be considered regarding the surface energy balance in high relief mountains. We then discuss debris load properties and debris transport dynamics. We finally address various feedback mechanisms and system couplings on the glacier surface with respect to ablation, supraglacial water bodies, ice cliffs, and surface topography.

CLIMATE FORCING
Radiative and precipitation forcing govern glacier dynamics, although the partitioning of these forcing components in any particular geographic area is not well-known with any degree of certainty, given a multitude of multi-scale topographic effects and our inability to accurately characterize and account for radiative and precipitation processes (Bishop et al., 2019;Bush et al., 2020). For example, it is thought that radiative forcing dominates in the Eastern Himalaya, while precipitation forcing has the biggest influence in the western Himalaya . Nevertheless, we require a much more detailed understanding of climate-glacier dynamics from an internal and external scaledependent forcing perspective, as various aspects of climate, glacial, and geomorphological systems need to be accounted for. Unfortunately, we lack basic data regarding radiation transfer and precipitation components that govern glacier states, and we do not know the responses and sensitivity to climate change for many regions of the world Dobreva et al., 2017;Bishop et al., 2019;Bush et al., 2020;Farinotti et al., 2020). Furthermore, the data sets that we tend to rely upon for such information (e.g., climate simulations and re-analysis data) do not accurately characterize local-scale radiative and precipitation fluxes due to inadequate representation of topography and/or processes (Bush et al., 2020).
Energy balance calculation at the glacier surface (with or without debris) is the first step toward determining the low-level air temperature above the glacier, but it neglects potential temperature advection by atmospheric circulations (e.g., anabatic, katabatic, or frontal winds, which all advect temperature). Furthermore, during the daytime, heating of the debris cover itself can generate atmospheric convection above the debris leading to strong converging surface winds and potential precipitation, which can further alter the surface energy balance. Consequently, numerous radiative-transfer (RT) parameters, circulation dynamics and cloud properties, and convective, cyclonic and orographic precipitation mechanisms must be considered to understand glacier response and sensitivity to forcings.

Irradiant Fluxes
The short-wave and long-wave radiant energy components of the surface-energy balance strongly regulate ablation dynamics. The short-wave surface irradiance (E) components include the direct irradiance from the sun (E b ), the diffuse-skylight irradiance due to atmospheric scattering (E d ) and the adjacentterrain irradiance (E t ). All of these components are space-timetopography-wavelength dependent, and numerical modeling is the only way to account for the inherent scale-dependent operation of orbital, atmospheric and lithospheric processes that govern the variability of the magnitude of these RT parameters over glaciers (Bishop et al., 2019). These radiation fluxes significantly contribute to ablation at a particular location depending upon the surrounding land cover and topographic conditions. Studies have estimated that the shortwave radiation usually accounts for 75% or more of total irradiance on glacier surfaces (Oerlemans and Klok, 2002).

Short-Wave Direct Irradiance
The direct irradiance is the dominant short-wave irradiance component (Proy et al., 1989;Olson and Rupper, 2019). It is strongly controlled by local and meso-scale topographic effects. It can be represented as: where λ l is wavelength, E 0 is the exo-atmospheric irradiance adjusted for Earth-Sun distance, T ↓ is the downward atmospheric transmittance of radiation governed by atmospheric absorption and scattering, and θ s is the solar zenith angle. Transmittance also accounts for the atmospheric optical depth that is partially governed by the relief structure of the landscape (i.e., altitude). The cosine of the incidence angle (cos i) characterizes the local topographic effects on irradiance given solar and terrain geometry relationships. Specifically, the incidence angle of illumination (i) between the sun and normal to the ground surface is defined as: where θ s is the apparent solar zenith angle that accounts for altitude variations and atmospheric refraction given the atmospheric temperature and pressure profiles, θ t represents the terrain slope angle, φ s is the solar azimuth angle, and φ t is the terrain slope-azimuth angle. Calculation of cos i is possible with the use of a digital elevation model, and negative values must be corrected to 0.0. Local topographic properties can have a significant effect on the magnitude of the direct irradiance that reaches the surface of a glacier. Ablation and meltwater transport, and gravitational sediment fluxes all alter the altitude, slope and slope azimuth of the ice and sediment surfaces of a glacier. A feedback between ablation-topography and surface irradiance controls the topographic evolution of a glacier surface . The importance of the local topographic conditions on the direct irradiance cannot be overlooked. The cosine of the incidence angle can range from 0 to 1 based upon slope angle and slope-azimuth angle variations caused by sediment transport and ablation processes. Based upon our simulations in the Karakoram (see Figure 1), direct irradiance values can be greater than 1,000 Wm −2 . Consequently, rapidly evolving topography should be accounted for as this sub-parameter strongly regulates ablation, all other factors such as debris properties being equal. For many valley glaciers, the topographic shading on the direct solar radiation can be significant due to complex topography (Olson and Rupper, 2019). The S parameter in Equation (1) accounts for the meso-scale relief structure of the topography that governs the presence of cast shadows on the landscape. Cast shadows on the landscape are governed by the size of the solar disk and the Earth-Sun distance. More specifically, the solar angular width (α s (t), [degrees]) in relation to topographic relief, will govern the overall length of a cast shadow in the direction of φ s . This angle and other geometry parameters can be used to compute the planimetric length of the umbra and penumbra subregions that collectively represent the total castshadow region. The umbra is the region over the landscape where E b is totally obstructed by the topography. Conversely, the penumbra region receives a fraction of E b because a portion of the solar disk is obstructed by the topography. In the penumbra region, the fraction of E b is governed by the pixel location. Consequently, a location in the umbra region would exhibit a shadow coefficient of 0.0, while a location in the penumbra region would exhibit a coefficient value ranging from 0.0 < S < 1.0. Therefore, cast shadow locations can also include irradiance from E b . Furthermore, because of the eccentricity of the Earth's orbit around the sun, the length of the cast shadow region varies nonlinearly with time. Given the extreme relief in many mountain environments, cast shadows can significantly modulate direct irradiance over the course of a day, let alone over the course of an ablation season. The importance of cast shadows should also not be underestimated, as the S parameter strongly regulates the direct irradiance during the early morning and late afternoon, as governed by basin relief and solar geometry. Large areas of a glacier's surface can be in cast shadow (umbra) for hours during the day, thereby blocking the direct irradiance and decreasing the surface irradiance and ablation potential. It is more difficult to evaluate the significance of the penumbra on ablation, other than to say that the length of the penumbra in the direction of the solar azimuth angle causes a direct irradiance gradient that influences ablation depending upon its length. Numerical modeling and sensitivity analysis are required to have a better understanding of the entire S parameter. Consequently, the S parameter can have a significant to insignificant influence on glacier ablation depending upon environmental conditions and the total length of the umbra and penumbra and time that cast shadows are found over the glacier surface. Its modulating effects on ablation will increase during the ablation season as solar zenith angles increase.
Simulating the magnitude and distribution of E b on a glacier surface is relatively straight-forward, given that topographic evolution resulting from ablation dynamics is accounted for with respect to local topographic properties and relief. We demonstrate relatively accurate simulated estimates of diurnal direct irradiance (1 day; July 24, 2004) using a SRTM 30 m DEM, and compared our simulation to field-based surface irradiance measurements collected on the same day by Mihalcea et al. (2008) for the same location on the Baltoro Glacier (Figure 1). We FIGURE 1 | Simulated short-wave direct-beam irradiance (E b ) component (wavelength range: 0.3 µm to 3 µm) for the Baltoro Glacier, compared to field measurements of short-wave irradiance at the same location on July 24th 2004. The discrepancy in magnitude is mainly caused by the diffuse-skylight irradiance and the adjacent-terrain irradiance.
accounted for orbital parameters based upon Berger (1978) and used the apparent solar zenith angle accounting for parallax and atmospheric refraction, assuming the temperature and pressure profiles of the atmosphere. In this simulation, however, we did not utilize ray-tracing to account for cast shadows. The simulation results are reasonable, the modeled E b is less than measured E values in the field, as we did not account for E d and E t that is represented in the E values. Perhaps it is more important to notice the discrepancy in magnitude in the early morning and late afternoon, caused by cast shadows (i.e., shape at the base of the distribution is different). Rapidly changing magnitudes at these times clearly reveals the influence of cast shadows in E, vs. simulated E b values. We would expect these differences to become more pronounced during the later stages of the ablation season as the solar zenith angles increases.
Challenges in estimating this dominant component more accurately involve better characterization of atmospheric conditions, use of higher-resolution DEMs for establishing initial conditions, and the development and evaluation of different parameterization schemes to account for those controlling parameters that govern the energy-balance at the glacier surface.

Short-Wave Diffuse-Skylight Irradiance
Atmospheric scattering will produce a hemispherical source of irradiance that contains isotropic and anisotropic components (Perez et al., 1986;Proy et al., 1989). The anisotropic nature of E d has been the largest source of error for estimating this component, as it is governed by circumsolar brightening due to forward scattering by aerosols, and horizontal brightening due to Rayleigh scattering (Perez et al., 1986). A Two-stream approximation model assumes that there is a Rayleigh-scattering component (E r ), an aerosol-scattering component (E a ), and a secondary ground/sky-backscattering component (E g ) caused by multiple interactions between the ground surface and the atmosphere (Proy et al., 1989;Gueymard, 1995;Zhang et al., 2015). One can account for single-scattering or a multiplescattering Rayleigh atmosphere. The single-scattering albedo is required to account for aerosol scattering which is a function of both wavelength and humidity. Aerosol type and models such as rural, urban, maritime or tropospheric can also be accounted for Gueymard (1995). The ground backscatter component can be modeled by accounting for the zonal ground reflectance from the overall ground and sky reflectance (Gueymard, 1995). Numerous parameterization schemes are available that enable a reasonable characterization of this irradiance component for a horizontal surface (e.g., Bird and Riordan, 1986;Gueymard, 1995). Multiscale topographic effects, however, need to be considered at each location on the landscape (Proy et al., 1989;Bishop et al., 2019).
Both local and meso-scale topographic properties govern E d and Proy et al. (1989) provides a computation solution such that: where L ↓ is the downward radiance from sky-hemispherical incident directions θ i and φ i , and θ i and φ i are the zenith and azimuth angles from the hemisphere, α s is the solarelevation angle (α s = π/2 − θ s ), and I is the incidence angle of the direction defined by sky-hemisphere and terrain geometry, similar to Equation (2) using incident geometry. Such a computation is computationally expensive given that E d is wavelength dependent, and hemispherical topographic shielding must be accounted for by limiting the integration based upon the maximum relief angle for all hemispherical azimuth directions. Basin relief production is primarily governed by paleo-glacier erosion and uplift history. There is a time disconnect such that the majority of the relief production occurred in the past, however the resulting relief (i.e., topographic shielding) partly influences the modern-day diffuse-skylight irradiance. Its altitudinal spatial variability can therefore govern the spatial variability in ablation. We demonstrate this topographic control using radiative and ablation simulations over the Baltoro Glacier in the central Karakoram (Figure 2). For this simulation, we did not perform the full integration of E d because of the computational intensity of doing so over all wavelengths. Rather we used the algorithms from Bird and Riordan (1986) to estimate E d on a horizontal surface, and then multiplied the results by the so-called skyview-factor coefficient that represents the relative magnitude of topographic shielding. The computation of this topographic parameter is highly scale dependent, and Bishop and Dobreva (2017) demonstrated that it is highly variable across mountain environments given the polygenetic nature of erosion dynamics and topographic evolution (for more parameter details see Dozier et al., 1981;Wilson and Bishop, 2013). Regarding the ablation simulations, we account for debris thickness distribution over the Baltoro Glacier from Mihalcea et al. (2008) and meltwater ponds to generate a temporally-averaged ablation rate over the glacier. We then compared the skyview-factor coefficient to simulated ablation over the glacier (Figure 2). Clearly there is a decrease in the FIGURE 2 | Topographic shading effect on the simulated surface ablation rates for the Baltoro Glacier in 2004. The ablation rates are averaged over 30 m altitude bins. The skyview factor represents the topographic shielding that governs the diffuse-skylight irradiance. The abrupt drop in the skyview is associated with the minimum ablation along the profile, which highlights this meso-scale topographic effect on glacier surface ablation.
ablation rate with altitude that accounts for the collective influence of cooler atmospheric temperatures and an increase in topographic shielding (lower skyview coefficient). At about 4, 400 − 4, 500 m (below Concordia; Concordia is at about 4,691 m) topographic shielding begins to decrease and is associated with the confluence of glacier tributaries. It then abruptly increases with maximum topographic shielding being spatially coincident with the minimum ablation rate for the profile. These simulation results suggest that meso-scale topographic effects also play an important role governing the magnitude of E d and ablation depending upon paleo-glacier erosion dynamics and basin relief production (i.e., erosion-uplift dynamics).
The significance of this secondary surface-irradiance parameter on ablation is difficult to ascertain, as numerical simulation sensitivity analysis focused on irradiance partitioning is sorely needed. We have demonstrated, however, that mesoscale topographic effects on the diffuse irradiance potentially regulate ablation rates based upon glacier erosion histories (i.e., relief production). Simulation results indicate that the magnitude of this parameter can range from relatively low or moderate energy levels in the visible portion of the spectrum (generally 5-200 Wm −2 ) depending upon topographic conditions. Spatial distribution patterns of higher diffuse energy are most likely found at higher altitudes where basin relief significantly decreases. Furthermore, the influence of landcover conditions can also increase the magnitude of this component given high albedo surfaces (i.e., snow) and secondary scattering back to the surface.

Short-Wave Adjacent-Terrain Irradiance
The aforementioned irradiance components and surface Bi-Directional Reflectance Distribution Function (BRDF) interact with the surrounding terrain geometry to produce the adjacentterrain irradiance (E t ). This irradiance component is extremely complicated but important to characterize, because variations in land cover and topographic complexity strongly regulate its magnitude. Highly reflective materials such as felsic minerals, vegetation, and ice and snow coupled with steep slopes can cause significant irradiance that impacts ablation dynamics, and especially the retreat of ice-cliffs (Sakai et al., 2002;Buri et al., 2016).
The characteristics of surface reflectance properties can range between diffuse and specular (Iqbal, 1983;Zhang et al., 2015). We know, however, that the isotropic assumption of reflectance is not valid in mountains, and the surface BRDF (ρ brdf ) must be accounted for at each location on the landscape, as reflectance from the surrounding terrain and other glacier surface locations can contribute to E t . The BRDF is used to describe the anisotropic nature of surface reflectance characteristics. It is a scattering function that describes anisotropic reflectance given all inputoutput angles. It is characterized as: where, L is the surface radiance and E is the surface irradiance are a function of solar geometry and terrain geometry. Consequently, the topography has a significant affect on the BRDF (Zhang et al., 2015), and the intimate mixture of materials and structure of the surface also influence scattering direction. The BRDF, however, is difficult to measure accurately, and is also the basis for the computation of surface albedo, another important RT parameter that is difficult to estimate.
The irradiance component E t can be characterized in the following fashion: where L is the adjacent surface reflected radiance coming from the effective incident direction, φ i is the hemispherical incident azimuth angle, θ i is the incident vertical hemispherical zenith angle to account for terrain radiance above and below a pixel location, T ↓↑ is the atmospheric transmittance given the optical depth of the atmosphere due to relief and propagation zenith angle through the atmosphere (θ v ) between two locations, cos I t represents the terrain incidence angle given the influence of the local terrain geometry in relation to the incident directional geometry, and S t representing terrain blocking of the adjacent surface radiance between any two points of the surrounding terrain. As a general rule-of-thumb, the computation of E t should account for the terrain conditions extending out to some distance (e.g., 5 km) from each pixel location (Proy et al., 1989). However, the appropriate distance in any direction can vary as a function of the anisotropic nature of the topography, given the orientation of the structural fabric of the topography related to differential erosion and tectonics. As such, the adjacent hemispherical distribution of locations that can contribute to E t for any particular location can be extremely scale dependent, given that it is also a function of altitude, such that the contributing zones of adjacent-terrain irradiance may be spatially contiguous or spatially fragmented. To our knowledge, this level of complexity has not been evaluated with respect to determining its impact on ablation dynamics, although it most likely plays a role in ice-cliff retreat when considering the coupled influence of the long-wave adjacent-terrain irradiance term.
The short-wave adjacent-terrain irradiance is thought to be a tertiary surface irradiance component that does not significantly affect ablation. However, given our parameterization scheme and the nature of highly directional and diffuse reflectance components associated with the BRDF, coupled with changing glacier surface topography, solar zenith and solar azimuth angles, the E t component is likely to be highly variable in space and time, and facilitate ablation in general. The significance of this component could be greater than the diffuse irradiance given specific topographic and adjacent land cover conditions, such as at depression areas that are surrounded by highly reflective felsic minerals and snow.

Precipitation
Precipitation in mountainous regions is partly modified by avalanche redistribution of snow, and is determined by uplift in traveling frontal zones or by orographic uplift. The former will move spatially with the frontal system while the latter will remain stationary given the topography and the wind direction. Determining where near-constant orographicallyforced precipitation occurs may lead to an explanation of why a particular glacier exists at one location but not in another nearby location.
Separating orographic precipitation from frontal precipitation is difficult from observations (in which all precipitation is combined). For orographic precipitation one needs to know the wind direction, the atmospheric moisture content, and the topography itself but one also has to track the properties of the air parcels along their trajectories to keep track of the amounts of condensation, evaporation, and precipitation as the air travels up and perhaps over multiple peaks, potentially losing all of its water content along the way. Some attempts at quantifying this process have been made using a combination of a high resolution DEM and downscaled reanalysis climate data (e.g., Dobreva et al., 2017).
Regardless of the driving mechanism behind mountain precipitation, the precipitate-formation process is mostly a coldbased one given typical surface temperatures at high elevations and the adiabatic lapse rate, i.e., the precipitate-formation process occurs at temperatures below freezing and so is governed by the Bergeron-Findeisen process. According to this process, in a cloud that has a mix of ice crystals and supercooled water droplets, water vapor in the cloud preferentially deposits on the ice crystals because of the lower saturation vapor pressure over ice than over water. As the water vapor decreases through deposition on the ice, the supercooled droplets evaporate and lose their water to the growing ice crystals, which may then grow large enough to overcome any updrafts and fall to the surface. This process dominates mid-latitude precipitation (even during the summer) and at high altitudes.
Most precipitation in mountainous regions therefore begins frozen; whether it remains frozen on its way down to the surface is dependent on the atmospheric temperatures through which it falls. At high altitudes, it is more likely to remain frozen than, for example, at sea level in mid-latitudes, but if conditions are warm enough then it will melt completely and fall as rain (or, if the precipitate does not have enough time to melt completely it will fall as graupel, sleet, or freezing rain).
Snowfall above the equilibrium line altitude (ELA) is paramount for a glacier's survival. Climate change can alter the ELA up or down, depending on local conditions, but with a general warming trend ELAs on average are rising. If the ELA rises above the peak of the mountain then no glacier can exist. If atmospheric conditions above the ELA become conducive to rain then glacier ablation will be accelerated.
When debris cover is present additional factors must be taken into account. If the debris temperature is below freezing, snowfall would accumulate as on a debris-free glacier. If the debris cover is above freezing, however, then it will melt the snow leading to percolation of liquid water onto the glacier surface that can enhance ablation. Refreezing of the percolated water at the glacier surface is possible, with the release of latent heat into the interstitial spaces of the debris cover (Collier et al., 2014(Collier et al., , 2015. The debris cover fundamentally alters the fluxes of sensible and latent heat to the overlying atmosphere (Collier et al., 2014). In addition, a relatively warm debris cover (through daytime solar heating) can trigger convection in the atmosphere above and produce precipitation that melts on contact and increases surface ablation (Collier et al., 2015). So DCGs have the potential to increase their own ablation rates through dynamic interactions with the overlying atmosphere as well as the underlying glacier.

Composition and Properties
Investigators are increasingly recognizing the significance of the supraglacial debris load as a major controlling factor in regulating ablation and the morphological evolution of DCGs (Adhikary et al., 2000;Nicholson and Benn, 2006;Reid et al., 2012;Rowan et al., 2015;Anderson and Anderson, 2016;Zhang et al., 2016;Gibson et al., 2017). Field data and numerical-modeling efforts (given various assumptions) demonstrate the general relationship between debris thickness and ablation rate, as thicker debris loads decrease the thermal energy available for ablation at the debris/ice interface, while thin and moisture laden debris can significantly accelerate ablation beyond that of debris-free conditions (Mattson, 1993;Kayastha et al., 2000;Reznichenko et al., 2010). This general relationship, however, does not account for a variety of debris-load properties, processes, and feedbacks regarding ablation dynamics, as supported by findings related to the "debris-covered glacier anomaly" (Salerno et al., 2017) and numerical modeling efforts (Anderson and Anderson, 2016;Zhang et al., 2016;Gibson et al., 2017;Huo et al., 2020). Consequently, the multi-faceted nature of this controlling factor is not adequately understood, as there are many unknowns with regard to properties and processes that govern ablation dynamics, and research clearly reveals the relatively high spatiotemporal variability in debris properties and processes that should be accounted for Scherler et al. (2011b); Reid et al. (2012); Rowan et al. (2015); Yue et al. (2017). The complexity of the problem is highlighted by the relatively large number of debris, near-surface and scale-dependent topographic parameters that partially govern ablation. Table 1 lists various physical and radiative-transfer properties that all contribute toward regulating the magnitude of ablation in space and time. It is also important to realize that location/scale specific parameters related to the topography and lithology beneath and surrounding a glacier also contribute to spatial and property variations in the debris load, as debris sources and mineralogical composition of debris and rocks is governed by the nature of the geological setting, and the geomorphological system within the glacial basin. An important first-order property that requires consideration is the particle size/distribution of the debris load. This property can be highly spatio-temporally variable over a glacier, and is regulated by source and sediment-transport dynamics. It can have a significant impact on ablation dynamics, as "particles" can include dust, fines, gravel, cobbles, rocks, isolated boulders, and boulder fields (Figure 3; the full range of debris particle sizes). Such heterogeneous conditions are common on many glaciers around the world, and "particle" distributions can be significantly different on glaciers within a region, depending upon paleoclimate-glacier dynamics (Bush et al., 2020). Such variations are the result of polygenetic sediment-transport dynamics involving mass movements, aeolean, suprafluvial, englacial, gravitational, and ice-flow transport mechanisms. Spatio-temporal variations of this aspect of the debris load are usually not considered, and homogeneous unconsolidated debris conditions are assumed, which causes uncertainties in ablation estimates.
The composition of the debris load is another important property. It is rarely accounted for, although it can be taken into consideration by using estimates of the debris bulk conductivity. It not only governs the bulk thermal properties of the debris, but also the surface reflectance properties and the surface energy balance of the short-wave irradiance (i.e., albedo). Consequently, remote sensing analysis and absorption-feature modeling can be used for mineral detection and determining the composition of the debris using hyperspectral data. This approach can provide more explicit estimates of the surface composition and distribution of moisture-laden debris. Such data regarding the debris load are too logistically difficult to acquire in the field, but high-resolution hyperspectral remote-sensing capabilities are now available using uncrewed aerial systems (UAS) to better assess "particle" distributions, surface mineralogy and other parameters that are not usually accounted for (e.g., Fyffe et al., 2020).
Debris load thickness has been shown to be highly variable (Mihalcea et al., 2008;Fyffe et al., 2020). Many glaciers exhibit low-frequency debris thickness variation (the generally thicker debris toward the terminus) controlled by surface ice-flow transport dynamics and englacial transport of sediment and rocks into the terminus region (Anderson and Anderson, 2016;Wirbel et al., 2018). Higher-frequency debris thickness variations must also be considered, as the debris load is modulated by suprafluvial transport, topography-driven gravitational transport, and other ablation dynamics including supraglacial lake development and evolution. Furthermore, a high percentage of debris cover lowers the glacier surface albedo that partially counters the debris insulation effect (Fyffe et al., 2020). Such complexities have not been adequately investigated, and it is essential to note that complex feedbacks exist between sediment fluxes, ablation, topography and the surface-energy balance . Although several approaches have been developed to estimate the distribution of debris thickness (e.g., Mihalcea et al., 2008;Zhang et al., 2011;Reid et al., 2012;Rounce et al., 2018), higher-frequency variations in debris thickness and the spatial variability of ablation is not known with any degree of certainty on most DCGs.
Field measurements suggest that debris temperature gradients are approximately stable and linear with depth during the ablation season (Rowan et al., 2021), although the linearity and stability depend on the time scale, and linearity is also dependent upon the debris thickness. Assuming linearity, debris surface temperature can be used to estimate sub-debris ablation, which is also controlled by debris thickness, and thermal conductivity (Nakawo and Young, 1981;Nicholson and Benn, 2006). Given shallow debris thicknesses, relative homogeneity in particle distribution, and lack of moisture-laden sediment, this first-order approximation assumes isotropic conditions with respect to various properties. With greater debris load thicknesses, however, porosity variation has the potential to alter the thermal gradient given the translocation of the fine size-fraction of particles and presence of moisture-laden sediment. Therefore, the thermal conductivity of debris can be difficult to determine given the complex composition and pore-space conditions. Nicholson and Benn (2006) presented an approach to address this issue using measured temperature profiles. Based on this model, the bulk thermal conductivity of the debris layer (k d ) can be written as: where c d , ρ d , c v , ρ v are the specific heat capacity and density of debris, and specific heat capacity and density of the porespace filler, respectively. ϕ is the porosity, and κ is the apparent thermal diffusivity of the debris, which can be determined using the thermal diffusion equation: where T d is the debris temperature, which is a function of depth (z) and time (t). κ can then by determined from the measured temperatures-time and temperatures-depth variations. Another method was presented by Reid and Brock (2010), based on which, the debris was broken down to multiple thin layers to estimate the thermal conductivity and internal temperature profiles: Debris moisture conditions, mineral mixing, debris surface structure, and surface topography can also significantly influence the surface albedo that regulates debris temperature and ablation (Pelto, 2000;Nicholson and Benn, 2006). Accurate estimation of the albedo based upon the surface BRDF is difficult to compute, and rapidly changing glacier-surface dynamics should cause the surface albedo to exhibit relatively high spatio-temporal variability (Takeuchi and Li, 2008;Yue et al., 2017). Therefore, better estimates of debris surface reflective parameters are required using more reliable parameterization schemes. We demonstrate this using numerical simulations and accounting for the spectral reflectance of minerals that make of the debris load of the Baltoro Glacier. Specifically, we utilize a radiation-transfer approach (Coakley, 2003) to compute the broadband surface albedo (α) based on the mixing of sediments and water at each grid cell, such that: where E is the total short-wave irradiance and L is the reflected surface radiance, which, based on an isotropic (Lambertian) assumption, can be written as (Bishop and Colby, 2011): where r mc is the spectral reflectance of the mineralogical composite representing the debris, which is considered as a composite reflectance of different types of minerals. The debris lithological distribution over the Baltoro Glacier was based on Gibson et al. (2017), which indicates that the surface is mostly covered by a composite of gneiss (about 53%), granite (about 30%), and schist metasediment (about 17%). Simulated results in Figure 4 reveal significant spatial variations in the surface albedo over the lower ablation zone of the Baltoro Glacier caused by mineralogical and moisture variations. We performed the narrow-band to broad-band conversion of surface albedo using ASTER imagery based on the work of Liang (2001), and confirmed that our simulated values are very close to the remotesensing approach. In the future, the use of semi-empirical BRDF models should be used with multi-angle reflectance imagery to generate higher quality albedo estimates for glacier surfaces.

Debris Transport
Debris loads are mobile due to ice-flow, gravitational movement, and fluvial processes (Anderson, 2000;Benn and Evans, 2014;Anderson, 2016, 2018;Zhang et al., 2016;van Woerkom et al., 2019;Fyffe et al., 2020). Although the exact glacial debris production mechanisms have not been fully understood, most investigators consider the main sources are from headwalls, sidewalls, basal erosion and moraines (Anderson, 2000;Benn and Owen, 2002;Anderson and Anderson, 2016;van Woerkom et al., 2019). Studies have observed mass transport from landslides (Benn and Evans, 2014), avalanches (Benn and Lehmkuhl, 2000), and rockfalls from valley walls (Benn and Evans, 2014). Englacial debris loads and rock encapsulated by basal erosion also contribute to the total debris mass balance, and several models have been developed to investigate the transport of englacial debris loads (Anderson and Anderson, 2016;Wirbel et al., 2018). One of the first models of supraglacial debris flux was developed by Anderson (2000), which describes the diffusive debris flux originated from medial moraines in valley glaciers. Anderson and Anderson (2016) and Anderson and Anderson (2018) developed a similar model that couples with ice dynamics to investigate the longer-term evolution of debris cover. Moore (2018) presented a comprehensive study on the gravitational debris transport across the ice surface, as the topography evolves due to melt. A theoretical framework was developed for assessing slope stability and gravitational mass transport accounting for meltwater balance on ablating ice. Modeled results provided insights into the geometry of stable slopes as a function of debris thickness and texture. Another recent study by Fyffe et al. (2020) provides valuable field observations on the rates and forms of debris transport. van Woerkom et al. (2019) studied the debris transport from lateral moraines using high resolution DEMs and found that the lateral moraines contribute to debris thickening along the margin of the glacier surface.
Here we discuss two dominant processes that govern debris movement on a glacier at different spatio-temporal scales: the gravitational process that operates at smaller spatio-temporal scales given complex glacier topography, and the advectivediffusive processes that control debris flux over larger spatiotemporal scales governed by the ice-flow .

Gravitational Debris Fluxes
Gravitational debris movement (such as sliding and slumping) occurs on hillslopes and on glacier surfaces (Benn and Evans, 2014;Fyffe et al., 2020). Gravity-driven debris flux regulates the local thickness distribution of supraglacial debris during the ablation season when the surface topography is constantly changing under rapid melting, especially around melt hotspots such as supraglacial ponds and ice-cliffs. Field observations have identified sediment sliding or slumping off steepening ice-cliffs due to ice-cliff retreat and supraglacial lake expansion Miles et al., 2017).
Gravity, internal friction and basal resistance force govern the gravitational movement of debris. Based on the debris flow model by Chen and Lee (2000), the unit net force acting on a debris column, F, can be written as: where ρ d is debris density, g is the gravitational acceleration, z x and z y are the first derivatives of ice-surface elevation in the horizontal and vertical directions, respectively, k is the pressure ratio as defined by Chen and Lee (2000), u x and u y are debris velocity components, r u is the constant pore-pressure ratio, and φ is the dynamic internal friction angle of the debris. The first term on the right represents the gravitational force, the second term describes the internal friction of debris particles, and the third term represents the resistance force at the base of the debris column. Gravitational debris flux governs the local redistribution of debris thickness during the ablation season when the ice surfaces topography changes rapidly due to melt. The transport of a thin debris cover, however, is strongly controlled by more complex processes such as fluvial transport and rain wash (Fyffe et al., 2020).

Ice-Flow and Debris Advection
Ice-flow describes the glacier motion that governs the advective debris transport from the production zone to the terminus (Benn et al., 2012;Rowan et al., 2015;Anderson and Anderson, 2016;Wirbel et al., 2018). Therefore, ice-flow dynamics play a fundamental role in debris transport as demonstrated by recent numerical simulations (e.g., Anderson and Anderson, 2016;Wirbel et al., 2018). Ice-flow velocity is governed by multiple factors such as the ice thickness and basal water pressure, and many DCGs exhibit high flow velocities even though they have stable termini Quincey et al., , 2011. The surface velocity field can be determined from remote sensing analysis, estimating englacial ice-flow velocities can be challenging, although estimates are usually sufficient for computing particle trajectories if far enough away from the bed (about 80% of the ice depth). Therefore, most studies use ice-flow models (such as the shallow ice approximation) to solve the englacial velocity field (Herman and Braun, 2008;Bueler and Brown, 2009;Anderson and Anderson, 2016), which usually requires ice thickness measurements from geophysical surveys using seismic or radio sounding (McNabb et al., 2012). Modeled ice-flow velocities often suffer from high uncertainty (Farinotti et al., 2017), however, more accurate models have been developed to account for basal melt, temperature-adjusted ice properties and rugged bed topography (e.g., Bueler and Brown, 2009;Egholm et al., 2011).
Debris advection is governed by the ice-flow velocity field (Benn and Evans, 2014;Anderson and Anderson, 2016;Wirbel et al., 2018). A recent model was developed by Anderson and Anderson (2016), in which the englacial and supraglacial debris advection were modeled 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. Based on this , the rate of change of supraglacial debris thickness in the ablation zone can be written as (Anderson, 2000;Anderson, 2016, 2018): where h d is the surface debris thickness, t is time, c s is the nearsurface debris concentration, M s is surface ablation rate, and u s is surface ice 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. Figure 5D shows the simulated supraglacial debris advection based on Equation (13) given the initial debris thickness ( Figure 5A) and surface velocity ( Figure 5C). The transport of englacial debris and its contribution to the surface debris load is controlled by ice-flow dynamics, given that the vertical component of ice velocity usually points upward in the lower glacier where ice thickness increases toward the terminus. These processes have been described by Anderson (2000), and demonstrated in the simulations by Anderson and Anderson (2016) and Wirbel et al. (2018). The contribution of englacial debris load may explain the debris mass balance issue discussed by van Woerkom et al. (2019). The most recent englacial debris transport model was presented by Wirbel et al. (2018), in which both the advection and diffusion of englacial debris are accounted for: where c e is the concentration of englacial debris, D is the diffusion coefficient, u e is englacial ice velocity, and r represents englacial debris sources or sinks. Remote-sensing approaches have been used to map debris thickness distributions (Mihalcea et al., 2008;Juen et al., 2014;Gibson et al., 2017;Rounce et al., 2018); the results, however, suffer from high uncertainties due to the instantaneous nature of satellite imagery and the complex properties of debris. Further investigations on debris production and transport modeling are required to better address the spatial patterns of surface melt and supraglacial morphological conditions.

ABLATION DYNAMICS
Surface melt dominates the ablation of a DCG (Hambrey et al., 2008;Reid and Brock, 2014;Rowan et al., 2015). The melt rate is strongly affected by the properties of debris loads and the presence of supraglacial water bodies due to their strong impact on the energy transfer dynamics (Sakai et al., 2000;Nicholson and Benn, 2006;Reid and Brock, 2010;Fyffe et al., 2014;Miles et al., 2016).
Field measurement of ablation is usually conducted through monitoring ablation stakes at multiple locations on the glacier surface. Ablation rates at locations with different debris thickness are measured by several studies, and field data showed that debris thickness is a predominant factor that governs the sub-debris ablation rate (Mattson, 1993;Kayastha et al., 2000;Rounce et al., 2015). Given the complex debris effects, empirical approaches such as the degree-day method for estimating ablation may be highly unreliable for DCGs (Braithwaite and Olesen, 1990;Nicholson and Benn, 2006). Therefore, physics-based surface energy balance models are widely used for estimating sub-debris ablation (e.g., Nicholson and Benn, 2006;Brock, 2010, 2014;Fyffe et al., 2014;Miles et al., 2016), and several remotesensing approaches have been developed to estimate the key parameters, such as the thickness and thermal properties of debris (e.g., Mihalcea et al., 2008;Zhang et al., 2011;Foster et al., 2012;Juen et al., 2014).
Several recent studies have successfully used the surface energy balance model to estimate ablation on DCGs. For example, Fyffe et al. (2014) used a distributed model to compute melt rate on the Miage glacier in Alps. Rounce et al. (2015) compared the influences of debris thermal conductivity, albedo, and surface roughness on the surface energy balance of the Imja-Lhotse Shar glacier in Nepal. These results highlighted that: (1) Unlike the thick debris layer, a thin debris is sensitive to air temperature variation and water content, therefore areas with wet thin debris must be identified. (2) Subdebris ablation is more sensitive to the changes in thermal conductivity than surface roughness and albedo. (3) It is important to account for the realistic topographic shielding and wind speed over the glacier surface in the energy balance model. (4) Many model parameters need to be calibrated using field measurements or remote sensing data as they can be highly site-specific. The energy balance model at the debris/air interface can be written as (Nakawo and Young, 1981;Nicholson and Benn, 2006;Zhang et al., 2011;Rounce et al., 2015): 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. Net radiation fluxes dominate glacier surface energy balance and significantly control surface ablation rates (Nicholson and Benn, 2006;Cuffey and Paterson, 2010). The net shortwave radiation flux and the net long-wave 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 long-wave 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 is responsible for the long-wave adjacent-terrain irradiance.
The energy balance at the debris/ice interface can be written as (Nakawo and Young, 1981): 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 Frontiers in Earth Science | www.frontiersin.org heat flux toward the ice that is not used for ablation, which is often negligible under a temperate ice assumption (Cuffey and Paterson, 2010). Most models assume constant heat storage and a linear debristemperature gradient, then Q m during the ablation season can be computed using the simplified heat-flux approach based on daily mean values at a minimum 24 h timescale (Nakawo and Young, 1981;Nicholson and Benn, 2006): where k d is the bulk thermal conductivity of the debris layer accounting for water in pore-spaces, T i is the ice temperature. The sub-debris melt rate (M s ) can then be computed as (Nakawo and Young, 1981;Nicholson and Benn, 2006): where L f is the latent heat of fusion for ice, and ρ i is the density of ice. Precipitation can also contribute to surface melt (Fujita et al., 2014;Fyffe et al., 2014;Miles et al., 2016), and several methods have been used to account for this process (e.g., Reid and Brock, 2010;Fujita et al., 2014;Miles et al., 2016). Although more complicated processes have not been accounted for, such as the energy exchange between rainfall and heated debris that could eventually lead to additional ice melt. The latent and sensible heat fluxes have been found to be less significant with respect to melt variability compared to the shortwave and longwave fluxes in different climates (Sicart et al., 2008).
Using the aforementioned equations, Figure 6A shows the simulated ablation rates on the Baltoro Glacier in the Karakoram, using surface conditions presented by Mihalcea et al. (2008). Note that the modeled values are reasonably similar to field measurements on several Himalayan glaciers (Nicholson and Benn, 2006), including Barpu Glacier (Khan, 1989), Khumbu Glacier (Kayastha et al., 2000), and Rakhiot Glacier (Mattson, 1993). Figure 6B is an example of the temporal variation in ablation rate as a function of debris thickness over an ablation season, which highlights the important role of the dynamics involving debris redistribution in regulating local ablation rates on the glacier surface, especially for areas with water bodies and ice-cliffs. Figure 7 depicts simulated surface ablation rates compared to the remote-sensing-based estimates (Mihalcea et al., 2008). Both results show suppressed ablation in the terminus region and higher ablation around inter-moraine valleys corresponding to the difference in debris thickness ( Figure 5A).

SUPRAGLACIAL WATER BODIES AND ICE CLIFFS
Supraglacial water bodies and ice-cliffs elevate the ice loss of DCGs (Anderson, 2014;Reid and Brock, 2014;Steiner et al., 2015; FIGURE 6 | The debris-thickness control on surface ablation rates. (A) June-averaged ablation rate as a function of debris thickness on the Baltoro Glacier (red) compared to in-situ measurements in the ablation season on other Himalayan glaciers (replotted), including Barpu Glacier (Khan, 1989), Khumbu Glacier (Kayastha et al., 2000), and Rakhiot Glacier (Mattson, 1993). Note that the melt-enhancing effect of thin-debris is observed on all glaciers. (B) An example of the temporal variation in ablation rate vs. debris thickness around a supraglacial lake on the Baltoro Glacier. The debris layer is getting thinner as sediments are gradually sliding into the lake, which leads to an increase in ablation rate (the decrease of ablation rate in late ablation season is due to the decreasing irradiance). Buri et al., 2016;Miles et al., 2016Miles et al., , 2018Thompson et al., 2016;Dobreva et al., 2017;Mertes et al., 2017;Huang et al., 2018). The number of supraglacial ponds and lakes have been found to be increasing on many DCGs, which is related to the glacier morphological changes such as differential surface lowering and collapse of englacial channel roofs (Sakai et al., 2000(Sakai et al., , 2002Benn et al., 2001;Quincey and Glasser, 2009;Gibson et al., 2017). These water bodies also play an important role in regulating water storage and drainage in a glacier system (Cuffey and Paterson, 2010;Benn and Evans, 2014), and cause englacial ablation via efficient heat transfer (Sakai et al., 2000;Gulley and Benn, 2007;Benn et al., 2012;Miles et al., 2016;Mertes et al., 2017).
Supraglacial lakes can be relatively large to moderate-sized water bodies (lake azimuthal distances have been observed to be up to 0.6 km in length on the Baltoro Glacier) that can be relatively short lived or exist over multiple years, and supraglacial ponds are smaller but more common water bodies that usually only appear during the ablation season. These water bodies are typically abundant in the ablation zone where the surface slope is gentle and heterogeneous surface lowering occurs (Reynolds, 2000;Sakai et al., 2002;Reid and Brock, 2014). The evolution of many supraglacial lakes and ponds are governed by filling and draining cycles, and hydrological connectivity that are controlled by glacier surface and internal structures (Benn et al., 2001Wessels et al., 2002;Miles et al., 2017). Studies also found that proglacial lakes also have a significant control on the glacier dynamics (Sutherland et al., 2020). These processes have only been studied on a handful of glaciers in the field, and there are still many processes and feedbacks that need to be investigated.
Supraglacial ponds and lakes are efficient in transferring heat into glacier ice due to their low surface albedo and active convection (Lüthje et al., 2006;Miles et al., 2016). Studies have estimated that the ablation rates around the lakes can be much higher than that of most debris-covered areas (Sakai et al., 2000;Thompson et al., 2016). The simulations by Miles et al. (2018) revealed that ponds may be responsible for 1/8 of total ice loss in the Langtang valley, Nepal. Furthermore, supraglacial water bodies will most likely exhibit accelerated growth on Himalayan glaciers given current atmospheric temperature trends (Benn et al., 2001). For example, Gibson et al. (2017) 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 tripled, which suggests that supraglacial water bodies can expand quickly on DCGs.
One of the first energy-balance models for supraglacial ponds on DCGs was developed by Sakai et al. (2000), which was used to address field observations on the Lirung Glacier in the Langtang Valley, Nepal. The core of this model is the energy-balance of a supraglacial lake that accounts for heat flux input and output due to meltwater, heat storage of the lake, and the bare-ice vs. debris covered areas beneath the lake surface: where Q is the net radiation heat flux on the lake surface, Q in is the input heat flux from meltwater inflow into the lake, which in most cases, can be neglected (Sakai et al., 2000), 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 (Sakai et al., 2000), Q t is the change in heat storage of the lake, Q i is FIGURE 7 | Surface ablation rates for the Baltoro Glacier modeled with two different methods: (A) The remote-sensing based approach (Mihalcea et al., 2008), which uses surface temperature data acquired by satellite. (B) The surface energy balance based approach, which uses modeled surface irradiance. Both results are based on the same debris thickness distribution, assuming no water bodies and were produced to represent rates on August 14 2004. Note the high spatial variability in ablation rate caused by the debris thickness distribution ( Figure 5A).
the heat flux for ice melt at subaqueous ice-cliff, and Q d is the heat flux for ice melt beneath the subaqueous debris layer. A more recent model was presented by Miles et al. (2016), which improves the Sakai et al. (2000)'s model in several aspects, such as the energy flux due to rainfall, more accurate turbulent fluxes, and the free convection of water that may circulate within the pond. In addition to the energy-balance, Miles et al. (2016) used a mass balance model to account for the changes in pond volume ( V): where V in and V out are the inflows and outflows, V i is the subaqueous bare ice melt volume, V d is the subaqueous subdebris melt volume, V le is the volume of vapor exchanges, and V r is the volume gained from rain. All terms in the energy-balance model and the mass-balance model have been discussed. Although some parameters are based on empirical relationships, this work provides the most comprehensive physics-based pond model to date. Figure 8 depicts the large number of supraglacial water bodies identified in the terminus region of the Baltoro Glacier. Our simulation results (Figures 8C,D) based on the above model suggest that supraglacial water bodies make a significant contribution to the total ice loss and surface lowering in the ablation zone, and exhibit a non-linear seasonal trend toward the end of the ablation season, as the extent and amount of lakes and ponds gradually expand to yearly maximum. The dynamics of supraglacial lakes and ponds is often coupled with the evolution of ice-cliffs around them. Studies have suggested that this cliff/lake system plays an important role in governing the mass balance of DCGs, and most field observations and modeling results confirmed that ice-cliffs can make a significant contribution to the total ice mass loss on a DCG (Sakai et al., 2002;Anderson, 2014;Steiner et al., 2015;Buri et al., 2016), although many questions remain unresolved (Sakai et al., 2002;Buri et al., 2016;Miles et al., 2016). Several forms of interactions exist between ice-cliffs and ponds/lakes. For example, lake undercut can cause the steepening of some ice-cliffs  (Sakai et al., 2000(Sakai et al., , 2002, and the melt of ice-cliff below the water surface is strongly affected by water fluctuations . Sakai et al. (2002) found that the orientation and inclination of ice-cliffs significantly control the evolution of ice cliffs given the solar geometry and cast shadows. The geometry of the icecliff also controls the thickness distribution of debris cover which regulates meltwater production. Ice-cliffs are also controlled by englacial processes, and studies have attributed the formation of some ice-cliffs to the collapses of englacial conduits (Sakai et al., 2000(Sakai et al., , 2002. Some recent simulations (Steiner et al., 2015;Buri et al., 2016) confirmed that the complex terrain around the cliff has a non-negligible effect on ice-cliff evolution due to the local shading and the adjacent terrain irradiance. Furthermore, icecliff evolution can also be governed by pond-water ablation that undercuts adjacent lake slopes. In addition, the ice topographic load coupled with the debris load causes significant variations in the ice stress fields, that partially control the process of calving.
We have observed this process on a very large supraglacial lake near the terminus of the Liligo Glacier during the summer of 2005. Significant water depths coupled with wind and wave action can cause significant undercutting, thereby altering the ice stress fields which enables calving and the production of icebergs in glacial lakes. The ice cliffs that are generated from these processes typically have very steep slopes, and the debris loads above them are quickly transported. Collectively, various feedbacks promote ice-cliff retreat and rapid supraglacial lake expansion, that most likely represents a non-linear ablation response to climate forcing.
Understanding supraglacial water bodies and ice-cliff dynamics is critical to assessing DCG sensitivity to climate change, and valuable knowledge has been gained from the aforementioned studies. However, there are still many issues that need to be investigated including the filling and draining cycles that control the evolution of many supraglacial ponds, the connectivity between surface water bodies and englacial FIGURE 9 | Diagram illustrating the interactions between key controlling factors discussed in this paper that govern the dynamics of a debris-covered glacier. The feedbacks (dashed lines) are important to the system as they may increase a glacier's sensitivity to climate change over time.
conduits, englacial ablation due to warm water outflow, rockfall from ice-cliffs, and modeling the edge effects and temperature gradients in ponds (Benn et al., 2001Miles et al., 2016Miles et al., , 2017.

FEEDBACK AND SYSTEM COUPLINGS
The dynamics of a DCG system is controlled by the complex interactions between climate, debris load, ice dynamics, hydrologic conditions and topography (Scherler et al., 2011a,b;Dobreva et al., 2017;Gibson et al., 2017). Unfortunately, few studies have investigated the couplings and feedback mechanisms in a DCG system, which is a major issue in characterizing DCGs, because isolating a small number of processes out of the full system increases uncertainty in results, as system couplings and feedbacks can significantly affect glacier dynamics (Rowan et al., 2015;Anderson and Anderson, 2016). Studies are often limited in spatio-temporal scale (Anderson, 2014;Anderson et al., 2021), because they rely on instantaneous remote-sensing data or field measurements at limited sites that do not describe multi-year or longer-term variations in glacier conditions (such as ice thickness, debris distribution and topography). Furthermore, many system couplings and feedbacks operate at very different spatio-temporal scales, making them difficult to investigate in the field. Therefore, numerical modeling is required to address the feedbacks and non-linear relationships between major glacier components. Figure 9 is a conceptual diagram showing the major systems couplings we identified between climate, mass balance, supraglacial water bodies, debris load, ice flow, and topography that govern the dynamics of a DCG.

Feedbacks on Glacier Surface
Figure 10 is a conceptual diagram that illustrates the feedbacks we identified between ablation, surface ponding and debris thickness. The most obvious positive feedback describes the relationship between meltwater production and supraglacial pond expansion: as a pond grows, its surface area and FIGURE 10 | Conceptual diagram that illustrates the positive feedbacks between ablation, surface ponding and the spatial distribution of debris thickness.
surrounding ice-cliff area get larger, they absorb more solar energy due to the lower albedo, causing enhanced melt and calving, which in turn, leads to further expansion of the pond. Another feedback involves the differential lowering of the glacier surface, the distribution of debris and the increased ponding of surface water: the spatially heterogeneous distribution of supraglacial water bodies and ice-cliffs governs the differential elevation change, which in turn, facilities the expansion of existing ponds and the formation of new ponds as suggested by multiple studies (Benn and Owen, 2002;Benn et al., 2012;Reid and Brock, 2014;Rowan et al., 2015;Salerno et al., 2017;Miles et al., 2018). The presence of supraglacial ponds and lakes also increases the spatial heterogeneity of debris thickness. Rounce et al. (2018) found that debris thinning is usually associated with the formation of supraglacial lakes and ice-cliffs, then the heterogeneous distribution of debris thickness further enhances the differential surface lowering which creates more depression zones for new ponds.
Positive feedbacks also exist between the lowing of surface albedo and enhanced melt, as a wet surface often exhibits lower albedo. Pritchard et al. (2008) 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, i.e., the "dirty ice" that causes more ablation (Fyffe et al., 2020), 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, causing faster ice-cliff backwasting than previously thought (Sakai et al., 2000;Benn and Owen, 2002;Reid and Brock, 2014).

The Role of Glacier Surface Topography
Surface topography is an important factor that acts as a bridge between multiple components of a DCG, and studies have identified the following topographic effects on DCGs: 1. Surface topography has a direct impact 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 (Arnold et al., 2006;Garg et al., 2017;Olson and Rupper, 2019).
2. Surface topographic conditions control supraglacial hydrology, especially the development of ponds and lakes (Huo et al., 2021). Water transport is enhanced with steeper gradients, while ponding is enhanced with lower gradients, therefore supraglacial water bodies are more abundant where the glacier surface is flatter (Reynolds, 2000;Sakai et al., 2002;Immerzeel et al., 2014;Reid and Brock, 2014). Topographic depressions caused by differential ablation and debris loading are also important, as they provide topographic sinks for meltwater to accumulate.
3. Topography also affects the distribution and transport of supraglacial debris. Observations and simulations have shown that the morphometric conditions of the sidewalls and headwalls of the glacial valley control sediment input to the glacier (Benn and Evans, 2014;Anderson and Anderson, 2016).
4. The adjacent-terrain irradiance due to complex topography coupled with debris cover contributes to the surface melt, and especially the retreat of ice-cliffs (Sakai et al., 2002;Wessels et al., 2002).

Basal Processes
The glacier bed is slowly evolving due to differential glacier erosion caused by variations in ice-flow dynamics, such that erosion is most likely at a maximum near the ELA due to relatively thick ice and high flow velocities (Herman et al., 2011;Steer et al., 2012). Another complexity is the basal hydrological conditions that govern sliding velocity, as high basal water pressure may lubricate the ice-bedrock interface (Cuffey and Paterson, 2010;Quincey et al., 2011). The higher sliding velocity leads to more frictional heating that further increases basal water pressure, which forms a positive feedback (Benn and Evans, 2014). Many details about the interactions between basal sliding, erosion, subglacial ablation, basal hydrology, and subglacial debris transport need to be investigated in future studies.
Collectively, the DCG system is complicated by interactions and feedbacks between debris load, ablation, topography, ice-flow, basal processes, and supraglacial water bodies. An integrated modeling of these processes may help explain field observations over High-Mountain Asia, such as the rapid lowering on Himalayan DCGs Immerzeel et al., 2013;Fujita et al., 2014) and the advancement of many glaciers in the Karakoram Bishop et al., 2014;Dobreva et al., 2017;Farinotti et al., 2020). Furthermore, it is also important to account for ice-flow dynamics being a potential cause of glacier thinning, as the re-distribution of ice mass, especially the declining ice discharge plays a dominant role in setting the thinning patterns on debris-free glaciers (e.g., Cuffey and Paterson, 2010), and has also been shown to be important for DCGs (Vincent et al., 2016;Brun et al., 2018;Anderson et al., 2021).
Improved parameterization schemes (more parameters and processes) suggest that climate glacier dynamics for DCGs are more complex than debris-free glaciers and that DCGs may be more sensitive to climate change than previously thought.
We recommend that more energy input components (e.g., the adjacent-terrain irradiance) and multi-scale topographic effects be accounted for in improved parameterization schemes. Our simulations indicate that scale dependent processes such as the gravitational debris transport and supraglacial ponding must also be accounted for given the coupling of climate and surface process and the non-linear responses in ablation and ice-mass loss. In all likelihood, our simulations probably underestimate the magnitude of non-linear responses due to these coupled processes/systems, as other processes such as calving were not accounted for. Furthermore, it is essential that future research focus on addressing questions about the significance of missing and existing processes so that we have an improved understanding of the degree of dominance and partitioning of key parameters and processes/systems that need to be incorporated into DCG models. This will require important comparative analysis of parameterization schemes and parameter sensitivity analysis to account for uncertainty in so many parameters and processes that govern climate-glacier dynamics.
More remote sensing and field-based studies are required to validate numerical modeling results. Further research is also required to determine if DCGs in different geographic locations exhibit similar or different sensitivity to climate change given variations in internal and external forcing factors. We speculate that such complex non-linear systems have the potential to reach tipping points based upon changing climate and geomorphological conditions within glacierized basins.

CONCLUSIONS
Debris-covered glaciers represent complex non-linear systems that operate given the coupling of climate, geomorphological and glaciological processes and feedbacks. Many important processes and systems have not been adequately characterized. Consequently, our current understanding of DCG sensitivity to climate change is most likely oversimplified, as numerous processes and scale dependencies have not been accounted for. This has resulted in conflicting views about these glaciers, such as the Karakoram anomaly and the "debris-covered glacier anomaly, " where the lack of data and information about climate forcing and ablation dynamics force us to rely on empirical measurements that represent a static snapshot in space and time. This paper and our simulations have attempted to provide a more integrated synthesis into the complexity of DCG systems by addressing important concepts associated with the key external and internal forcing factors and glacial processes, highlighting the importance of feedback mechanisms, system couplings, and recognition of the high degree of uncertainty related to numerous properties and parameters that need to be taken into consideration. Our simulation results based on the Baltoro Glacier in the central Karakoram combined with a review of literature reveal that climate-DCG dynamics and responses can be significantly attributed to the following factors: 1. Climate forcing: solar radiation and precipitation are the main driving forces for ice loss or gain. An accurate surface energy balance model for debris-covered glaciers needs to account for short-wave and long-wave irradiance. Our simulations showed that multi-scale topographic effects on surface ablation can be significant for debris-covered glaciers, given glacier surface topographic evolution and differential sky-view conditions in the glacier valley. Other meso-scale topographic effects such as cast shadows and adjacent-terrain influences should also be considered.
2. Supraglacial debris: physical properties of debris, such as thickness, lithology, grain size, moisture content and thermal conductivity strongly govern the surface ablation rates on debriscovered glaciers. Many of these factors, however, have not been adequately characterized or studied. We showed that there are significant spatial variability in debris thickness, albedo, particle size and composition, which are also coupled with the dynamic movement of debris due to gravity and ice-flow. Such complexity needs to be accounted for in models and field studies to better assess debris-related effects on glacier mass balance.
3. Supraglacial water bodies and ice-cliffs: supraglacial ponds, lakes, and surrounding ice-cliffs are melt hotspots on debriscovered glaciers and are expanding in size on many glaciers. Our simulations based on the Baltoro Glacier suggest that supraglacial water bodies make a significant contribution to the total ice-mass loss over the ablation season. Further investigations are needed to address the evolution of supraglacial water bodies and icecliffs, and their interactions with debris load, topography and englacial processes.
4. System coupling and feedback: positive feedbacks on debriscovered glaciers have been identified in numerical simulations, such as the surface melt-lowering-ponding feedback, and the atmospheric convection feedback. These couplings and especially the positive feedbacks may govern the non-linearity of the glacier's response to climate forcing, which are represented as an acceleration in ice-mass loss and heterogeneous surface lowering. The combined effect of these processes may lead to the beginning of a critical transition of the glacier system that signifies an increasing level of glacier sensitivity to climate change.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
DH conducted the model simulations, made the figures, and wrote the sections 3.2-6. MB proposed and designed the research, developed the parameterization schemes, wrote the sections 2.1 and 3.1, and edited the manuscript. AB wrote the section 2.2 and edited the manuscript. All authors contributed to the writing of the introduction and conclusions and the final editing.