Modeling Supraglacial Ponding and Drainage Dynamics: Responses to Glacier Surface Topography and Debris Flux Conditions

Supraglacial ponds play a significant role in the mass loss of many debris-covered glaciers in the Himalaya. Glacier surface topography and debris flux conditions are thought to govern supraglacial ponding and drainage. Existing studies, however, have not adequately investigated the relationships and feedbacks between meltwater production, debris transport, topographic evolution and ponding, because field measurements are limited in time and space, and most existing models either neglect these processes or use oversimplified assumptions. Such limitations restrict our understanding of supraglacial 5 hydrology and introduce uncertainties in our assessments of glacier sensitivity to climate forcing. This study develops a more comprehensive numerical model to provide insights into the couplings between topographically-controlled surface ablation, meltwater drainage, ponding, and gravitational debris transport under radiative forcing. We investigate supraglacial ponding and drainage dynamics in response to different topographic and debris flux conditions through numerical simulations based on Baltoro Glacier in the Karakoram and several hypothetical scenarios. Results suggest that: 1) Supraglacial ponds make a 10 significant contribution to the total ice loss (more than 20%) in the lower-mid ablation zone over one ablation season, which elevates the glacier’s nonlinear response to radiative forcing. 2) Gravitational debris transport has a non-negligible control on the growth rate of supraglacial ponds by governing debris thickness and ablation rates on the ice-cliffs around ponds. 3) Glacier surface gradient and local topographic depressions control pond formation by affecting supraglacial water storage and drainage. Our simulations provide a possible explanation to the abundance of ponds in the mid ablation zone where slope is 15 gentle and more local depressions are present. These findings may contribute to more accurate predictions of future glacier changes in response to climate change.


Introduction
Supraglacial ponds play an important role in glacier mass-balance and glacial hydrology (Fountain and Walder, 1998;Sakai et al., 2000;Wessels et al., 2002;Cuffey and Paterson, 2010;Miles et al., 2016Miles et al., , 2018Huo et al., 2021a), and they are likely 20 to grow rapidly on debris-covered glaciers (DCGs) given projections of atmospheric warming (Benn et al., 2001;Gibson et al., 2017), which will have a significant impact on regional water resources and hydro-power supply (Wessels et al., 2002;Dobreva et al., 2017;Bush et al., 2020). A typical lifespan for a supraglacial pond on Himalayan glaciers ranges from a few months to a few years, during which they can form, grow, merge and be completely drained when they intersect with englacial conduits (Benn and Lehmkuhl, 2000;Gulley and Benn, 2007). Many supraglacial water bodies are hydrologically connected 1 https://doi.org/10.5194/esurf-2021-53 Preprint. Discussion started: 19 July 2021 c Author(s) 2021. CC BY 4.0 License. (Watson et al., 2016;Benn et al., 2017;Miles et al., 2017), and most of them exhibit periodic filling and drainage processes (Wessels et al., 2002). Large supraglacial ponds can also be destructive, outburst floods originate from DCGs in the Himalayas have caused injuries, deaths and property damage to downstream villages (Reynolds, 2000;Richardson and Reynolds, 2000). Therefore, a better understanding of the processes and factors that govern supraglacial pond formation and evolution is critical for a more accurate assessment of regional water resources and geohazard conditions. 30 Supraglacial ponds and adjacent ice-cliffs are considered zones of rapid melting on DCGs (Sakai et al., 2000(Sakai et al., , 2002Reid and Brock, 2014;Miles et al., 2016;Thompson et al., 2016;Mertes et al., 2017;Huang et al., 2018;Miles et al., 2018). A better understanding of pond dynamics can provide valuable insights into the nature of glacier sensitivity to climate change (Reynolds, 2000;Gibson et al., 2017;Benn et al., 2017;Miles et al., 2018;Huo, 2020). Field studies in the Himalaya have suggested that the degree of glacier surface lowering is related to the spatial density of supraglacial ponds in the ablation zone due to more 35 efficient heat absorption and ablation caused by ponds (Benn et al., 2012;Miles et al., 2017). Most existing knowledge of supraglacial ponds is gained from field measurements or remote-sensing-based analysis, these field studies provide valuable insights into pond dynamics, however, they are usually limited in time and space, and do not permit a quantitative understanding of the relationship between ponding and multiple surface processes. For example, we do not know what percentage of icemass loss on a DCG is due to the presence of supraglacial ponds, and what are the factors and processes that control pond 40 expansion. Supraglacial ponding is also neglected in most existing glacier models due to the complexity of characterizing the coupling between meltwater production, morphological evolution and hydrological processes. These limitations often lead to an underestimation of ice loss on DCGs. This could be, for example, one of the reasons why current models cannot explain the increasing number of supraglacial water bodies and accelerated surface lowering observed on some debris-covered Himalayan glaciers (Kääb et al., 2012;Fujita et al., 2014;Huang et al., 2018). Therefore, an integrated model that accounts for the 45 evolution of supraglacial ponds and surface ablation dynamics is sorely needed to investigate the mass balance of DCGs and their sensitivity to climate change.
The development of supraglacial ponds on a DCG is not only governed by meltwater drainage and filling, but is also strongly controlled by topographic conditions and debris-flux. Studies have suggested that the morphological changes on a glacier surface contributes to the increasing number of supraglacial ponds (Quincey and Glasser, 2009;Benn et al., 2012), 50 For example, the undulating surface and the gentle slope of the ablation zone on a DCG encourages the formation of ponds in depression areas (Sakai et al., 2000;Benn et al., 2001;Wessels et al., 2002;Benn et al., 2017;Miles et al., 2017), and the expansion of ponds further lowers the slope that encourages meltwater accumulation. The distribution and movement of supraglacial debris also controls pond expansion. Studies have identified significantly reduced debris thickness due to ponding (Rounce et al., 2018), and sediment flow on ice-cliffs during their retreat (Benn and Owen, 2002;Sakai et al., 2002) , which 55 explains why the pond-facing slopes usually have a thinner debris cover that enhances melting, creating steeper slopes, and causing further retreat of ice-cliff (Sakai et al., 2002;Reid and Brock, 2014;Mertes et al., 2017). Unfortunately, the impacts of topographic condition and debris transport on supraglacial ponding dynamics have not been adequately characterized in existing models, which limits our understanding of multi-scale surface morphological evolution on DCGs.
to better understand supraglacial ponding and the associated controlling factors. Specific research objectives are: 1) Develop a numerical model for DCGs that integrates meltwater production, surface drainage, ponding, topographic evolution, debris transport, and the feedback mechanisms between them; 2) Quantify the potential impact of supraglacial ponds on ice-mass loss based on a case study simulation of the Baltoro Glacier in the central Karakoram; 3) Understand how gravitational debris flux controls supraglacial ponding; and 4) Investigate surface topographic conditions that govern supraglacial pond formation using 65 numerical simulations.

Data
The lower-mid ablation zone of the Baltoro Glacier is simulated to evaluate the model and to answer our research questions.
The Baltoro Glacier is a notable debris-covered glacier located in the central Karakoram in Pakistan, it is one of the largest 70 temperate glaciers in the world outside the polar regions. Baltoro Glacier is an ideal glacier for studying supraglacial ponding because quantitative estimates of surface ablation rates, debris-thickness distribution and surface morphological maps (including supraglacial ponds) over the Baltoro Glacier have been produced from field measurements and remote-sensing analysis (Mihalcea et al., 2006(Mihalcea et al., , 2008, and a recent study has revealed an increasing number of ponds on the Baltoro Glacier (Gibson et al., 2017).

75
Initial surface topography conditions and land-surface parameters were derived from a 30m digital elevation model (DEM), which was generated based on the use of stereo-correlation from ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) satellite imagery using SILCAST software. The initial surface temperature was acquired from the ASTER surface-kinetic temperature product derived from the same ASTER imagery (acquisition date: August 14, 2004, ID: 00308142004054614) in order to maintain temporal consistency and facilitate calibration with the field measurements ac-80 quired at that time by Mihalcea et al. (2006). The air temperature measurements in 2004 on the Baltoro Glacier (measured at 2m above the surface) have been published (Mihalcea et al., 2006), we extrapolate this data to cover the entire ablation season using the modeled mean diurnal trend of air temperature variation (Collier et al., 2013) and the monthly trend based on the NCEP-CFRS reanalysis 2-m air temperature product (Reggiani et al., 2016) in the central Karakoram region. A vertical lapse rate of 0.0065Km −1 is used following Reid and Brock (2010) to account for the decrease in air temperature with increasing 85 elevation.

Surface Ablation Under Radiative Forcing
Solar radiation and the thermal irradiant flux from the atmosphere are the main radiative-forcing components for melting ice (Cuffey and Paterson, 2010). We use a radiation-driven glacier ablation model because the traditional degree-day factor approach is highly site-specific and cannot adequately account for debris and topographic effects (Braithwaite and Olesen, 90 1990;Nicholson and Benn, 2006;Huo et al., 2020). Based on the transfer of solar radiation, the two energy fluxes received by the glacier surface are the net shortwave radiation flux (Q s ) and the net longwave radiation flux (Q l ). The total shortwave surface irradiant flux in our model consists of two components: 1) the direct solar-beam irradiance; and 2) the diffuse-skylight irradiance. They are computed following Huo et al. (2021a). In this distributed model, we classify the glacier surface into three types: debris-covered ice, bare-ice or supraglacial pond. A pixel is classified as supraglacial pond if the water depth above 95 debris surface does not fall below 0.5m. Ablation rates are calculated differently for these three surface types based upon surface energy balance.
The energy balance at the debris/air (or the ice/air) interface can be written as (Nakawo and Young, 1981;Nicholson and Benn, 2006):

100
where 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 or ice which governs ablation. Note that in this 2-D distributed model, heat flux terms are computed at the beginning of each iteration due to the changing topography, albedo, and surface temperature. The typographically-controlled net radiation fluxes (Q s , Q l ) at the glacier surface and surface albedo for debris-covered area are computed following Huo et al. (2021a), Q h and Q e are computed following Nicholson and Benn (2006).The surface albedo of supraglacial pond as a function of water depth 105 are estimated using the empirical equation reported by Taylor and Feltham (2004). Based on the energy balance at the ice/air interface, the melt rate for bare ice under temperate conditions can be calculated as: where ρ i is the density of ice, and L f is the latent heat of fusion for ice.
The sub-debris ice melting is governed by the energy balance at the debris/ice interface (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 conductive heat flux towards the ice that is not used for ablation. We assume the temperature at ice-debris interface is constant at 0 • C during the ablation season after Nicholson and Benn (2006), causing Q c to be negligible (Cuffey and Paterson, 2010).
Under the assumptions of constant heat storage and a linear debris-temperature gradient, Q m during the ablation season can 115 be estimated as (Nakawo and Young, 1981;Nicholson and Benn, 2006): and Benn, 2006): The ablation process beneath the pond surface is different from sub-debris or ice melting because pond water can store and transport energy. The energy-balance model for a closed supraglacial pond can be written as (Sakai et al., 2000):

125
where Q cw is the conductive heat flux into the ice beneath the water surface that causes ablation, and ∆Q t is the change in heat storage of the pond (Sakai et al., 2000): where c w is the specific heat of pond water, ρ w is the density of pond water, and they are held as constants in our simulations.
∆T w and ∆V w represent the changes in water temperature and pond-volume over the time step ∆t respectively. We assume 130 a turbulent well-mixed water for the entire pond, such that the water temperature does not change with depth (Thomsen and Reeh, 1986), also note that this model neglects englacial energy transfer because we only focus on surface processes in this study.
Because ice is at melting temperature while in contact with pond water, the melt rate beneath the pond surface can be computed as (Nakawo and Young, 1981;Nicholson and Benn, 2006): Collectively, the surface ablation rate M is computed differently for different surface conditions: Surface topography evolves by a denudation rate of ∆z i after each iteration (with a time step of ∆t) based upon the ablation rate and surface slope angle (θ t ): Note that this ablation model only accounts for the radiation-driven surface ice melt, and neglects other forms of iceloss, such as calving around the perimeter of a supraglacial pond caused by lateral ablation, which is a process that also contributes to the expansion of ponds, because calving is governed by other processes (e.g., gravitational ice movement, icemass stress caused by debris load, structure of the ice-cliff and ice velocity) that are beyond the scope of this study, and cannot 145 be characterized using the available data.

Debris Transport Around Supraglacial Ponds
Supraglacial sediment (debris) is mobile under gravity (Benn and Evans, 2014), especially during the ablation season, when the surface topography around supraglacial ponds is constantly changing (usually steepening) due to rapid melting (Sakai et al., 2002;Reid and Brock, 2014). This parameterization scheme neglects debris advection in ice-flow velocity field because 150 we focus on the local redistribution of debris over a relatively short period of time (one ablation season), during which, the movement of debris is governed by gravity and local topographic change. From a mass flux perspective, the change of debris thickness can be estimated as: where q in and q out are the debris fluxes into and out of a grid cell, respectively, and A c is the planimetric area of a grid cell.

155
The flux out of a grid cell can be represented as: where A s is the cross-sectional surface area which the flux passes through, and u is the movement speed of a debris column.
The momentum equation is used for solving u from force analysis:

160
where ∆u represents the change in the speed of a debris column, F is the net force applied on the debris column. A tuning parameter (β) is used here to account for surface conditions that control the flux rate, and different values are evaluated in this study. The bulk density of debris (ρ d ) is assumed to be a constant, and V is the volume of the debris block. Very thin debris cover is governed by more complex dynamics: the debris is strongly affected by meltwater and ice, where cohesion and refreezing effects can significantly restrict debris particle movement. For simplicity, the bare-ice condition is determined as Mihalcea et al. (2008).
To estimate the net force applied on each debris column, we adopt the force analysis for unsteady gravity-driven debris flow characterized by Chen and Lee (2000), which accounts for gravity, internal friction and basal resistance for each finite moving debris column. Based on this model, the unit net force acting on a debris column, F , can be written as (Huo et al., 2021a): 170 where ρ d is debris density, g is the gravitational acceleration constant, z x and z y are the first derivatives of ice-surface elevation in the E-W and N-S directions, respectively, k is the lateral Earth-pressure ratios in active state as defined by Chen and Lee (2000), u x and u y are debris velocity components (from previous iteration) in the E-W and N-S directions, respectively, r u is the constant pore-pressure ratio, and ϕ is the dynamic internal friction angle of the debris. Figure 1 shows a sensitivity test for 175 the β value, which controls the debris flux rate withing a reasonable range based on field observations (Huo et al., 2021b). In our numerical experiments, we vary the value of β to understand how different debris flux rates affect supraglacial ponding. Figure 2 is a diagram highlights the coupling between the local debris transport under gravity, topography and ponding, which is the main process investigated in this study.

180
In the ablation zone, numerous water channels form on glacier surfaces that allow efficient meltwater drainage (Cuffey and Paterson, 2010;Benn et al., 2017). A supraglacial drainage model that characterizes water-depth variation is necessary in order to simulate the formation of supraglacial ponds and to estimate meltwater runoff from the glacier surface. Based on the ablation model and the DEM of the glacier surface, we numerically simulate the supraglacial drainage and ponding conditions in the ablation season using the model presented by Lüthje et al. (2006a) and Lüthje et al. (2006b). This model is a second-order differential solution accounting for topography-controlled flow path, surface lowering due to melting, and variations in water level, and has been successfully utilized to model the evolution of supraglacial ponds and sea-ice melt ponds (Lüthje et al., 2006a,b). This model can be written as: where h w is water depth, t is time, He(h) is a Heaviside function to prevent negative water level, ρ w is the density of water, η 190 is the dynamic viscosity of water, Π h is the horizontal permeability of the debris layer that controls water-flow speed, we use a value for well-sorted gravel and sand mixture that is recommended by Bear (2013), and z i is ice-surface elevation which is updated at each iteration based on equation (10). This model is adopted because it accounts for the resistance effect caused by a porous top layer, e.g., the barriers of porous sea ice (Lüthje et al., 2006b), which is a good analogy of modeling meltwater passing through supraglacial debris.

195
The timestep for the simulations is determined using an adaptive timestep approach for finite differencing implementation to ensure stability, such that the timestep for each iteration is a dynamic function of the surface gradient (∇z i ), cell size (ds) and the remaining time (t r ) for the simulation:

Simulation Scenarios and Initial Conditions
Three groups of numerical experiments (Table 2) are simulated over the ablation season in 2004 to answer our research questions.

205
The first group of simulations (S1-S2) investigates the contribution of supraglacial ponds to the overall surface ablation and differential thinning. To quantify the ablation caused by ponds, we compare and contrast two scenarios: S1: a debris-covered glacier with supraglacial ponds, and S2: a debris-covered glacier without supraglacial ponds.
The second group of simulations (S3-S8) addresses the influences of gravitational debris flux and surface topography on pond expansion. Gravitational debris flux plays an important role in pond expansion because it controls the redistribution of 210 debris thickness on a changing topography, and supraglacial ponds evolution could be more sensitive to the debris-flux due to the steep ice-cliffs around them. Specifically, we simulate the expansion of supraglacial ponds on different topography and with different debris-flux rates.
The third group of scenarios (S9-S11) investigates how the glacier surface gradient influences supraglacial meltwater drainage and pond formation. These simulations were designed to provide modeling evidence to support the hypothesis that a 215 low gradient glacier surface facilitates the formation of supraglacial ponds (Reynolds, 2000). Specifically, we simulate surface drainage condition and estimate discharge for glacier surfaces with 3 different gradients at 2º, 5º, and 10º (by altering the slope of the best-fit-plane over the study glacier surface).
The initial surface topography, debris thickness distribution and meteorological measurement over the Baltoro Glacier were discussed in the data section. All scenarios start with a dry surface and we assume there is no precipitation, tributary glacier 220 input, or englacial processes throughout the simulations, such that all supraglacial ponds are created via the accumulation of meltwater generated on the glacier surface given radiation forcing. For scenarios S3, S4 and S5, the initial topography is an inclined surface with artificial topographic variations (a hummocky surface) and a 0.5m thick debris load. For scenarios S3', S4' and S5', the initial topography is changed to an inclined surface without any undulations. For scenarios S6, S7 and S8, a subset of the actual topography on the Baltoro Glacier is used.
225 Table 2. List of simulation scenarios used to provide insight into the impacts of glacier ponds on ice-mass loss and surface morphometry.
The β values for low, moderate, and high debris-fluxes are 0.01, 0.05, and 0.10 respectively. HU represents the hypothetical topography with undulations, and HP represents the hypothetical inclined plane.

Supraglacial Ponding and Drainage
Supraglacial ponding and drainage conditions control the melt water depth distribution on the glacier surface. We simulate the supraglacial water depth variations in the ablation zone over the ablation season in 2004 (pre-existing lakes and ponds are not accounted for). Figure 3 shows the initial ( Figure 3A) and final ( Figure 3B) water depth distribution on the glacier surface. The  Figure   3C), as reported by previous studies (Taylor and Feltham, 2004;Lüthje et al., 2006a). Figure 3D depicts the decrease of the overall surface albedo for the study area (spatially averaged over each 30m altitude bin) from June 1st to September 28th in 2004, we noticed that the mid ablation zone where ponds are more abundant exhibits the most significant decrease in surface albedo over ablation season. Simulation results produce higher estimates of the total area of water bodies as compared to the published remote-sensing analysis over the same area by Gibson et al. (2017), (1.74km 2 versus 1.07km 2 ). Many factors could contribute to this difference, both process and methodology related, for example, the quality of satellite data, and the use of thresholding can influence remote-sensing-based estimates. The dominant factor most likely is the neglection of englacial conduits in the model, as englacial conduits can rapidly reduce the volume of surface water through englacial drainage, but it is beyond the scope of this 245 study as we only focus on surface processes. Although there is a discrepancy between the simulation and the remote-sensing results, we consider our simulations represent a meaningful first attempt to model the spatio-temporal dynamics of meltwater ponding and drainage over the ablation zone of a DCG.

Surface Ablation Due to Supraglacial Ponds
To understand how supraglacial ponds govern surface ablation, we computed surface ablation rates over the ablation season.
250 Figure 4 compares the mean surface ablation rates and the temporal standard deviation between pond-present and no-pond scenarios (simulation S1 and S2). The mean surface ablation rates represent temporally averaged values over the ablation season (from June 1st 2004 to September 28th 2004), and the standard deviations are computed based on monthly mean ablation rates (temporal variability over the ablation season for each grid cell). The simulated sub-debris ablation rates have been validated using the field measurements, which was detailed in our previous studies (Huo et al., 2021a,b). This comparison 255 indicates that ponds create numerous melting hotspots on the glacier surface and increased the spatial heterogeneity of surface ablation. The temporal standard deviations reveal that the ablation rates for ponds also exhibit very high temporal variability over the ablation season, which suggests that supraglacial ponds are more sensitive to the variation in radiative forcing than other areas on the glacier surface.
To better highlight ponding and its impact on the glacier, we show the increase in total surface water volume over the ablation 260 season ( Figure 4E), and temporal comparisons between the two scenarios from two perspectives: The diurnal mean ablation rate ( Figure 4F), and the cumulative ice mass loss ( Figure 4G). Simulation results suggest that the presence of meltwater and pond expansion dramatically accelerate the overall ice loss in the ablation zone, and the magnitude progressively increases towards the end of the ablation season because of pond expansion, although the magnitude of ablation decreases after early September, which is most likely caused by the decrease in the magnitude of irradiance. We notice that the nature of cumulative 265 ice-mass loss is nonlinear, and we would expect that the nonlinearity can cause dramatic change over longer periods of time if the ponds remain on the surface over multi-year time frames. These results strongly suggest that the formation and evolution of supraglacial ponds represent a key indicator of a glacier's overall sensitivity to radiative forcing. Table 3. Spatial statistical comparison between simulations S1 (with ponds) and S2 (without ponds).Ms represents mean surface ablation rate, σM s represents the standard deviation of surface ablation rate,M l represents the mean ablation rate beneath pond surface, A l % represents the area fraction of ponds at the end of the 120-day simulation, and M l % represents the fraction of ponds' contribution to total surface ablation. The ablation rates are temporally averaged over the ablation season first, and then the statistics are computed over the glacier surface.  Table 3 compares the statistics of the two simulated scenarios, which demonstrates that the presence of supraglacial ponds significantly increased the magnitude and spatial variability of surface ablation. The ponds in late September (when they reach 270 their maximum extend) accounts for 7.57% of the surface area, but they contributed 22.78% of the total ice-mass loss over the ablation season, and the modeled ablation rates beneath the pond water are on average 5.6 times higher than the average for non-pond areas. This strongly suggests that supraglacial ponds can significantly increase ice-mass loss on a DCG, and contribute to the high-magnitude heterogeneous thinning observed on many DCGs, such as the Khumbu Glacier (Sakai et al., 2000), the Baltoro Glacier (Mihalcea et al., 2006), and the Koxkar glacier (Huang et al., 2018).

Debris Flux and Topographic Controls on Ponding
Supraglacial ponds expand during the ablation season, and the growth rate is thought to be linked with debris thickness distributions and topographic conditions around the ponds (Benn et al., 2001). To understand how debris flux and surface topography control pond expansion, we conduct simulations with three different gravitational debris flux rate and surface topography conditions (S3-S8).
280 Figure 5 depicts simulated ponding for simulations S3-S8, based on 1) a hypothetical topography consisting of equallyspaced periodic variations ( Figure 5A); 2) an inclined surface without any topographic undulations ( Figure 5B); and 3) a subset of the modern-day surface topography on the Baltoro Glacier ( Figure 5C). The comparison between S3 (low flux rate scenario), S4 (moderate flux rate scenario), and S5 (high flux rate scenario) indicates a positive correlation between gravitational debrisflux rate and pond size. A higher flux rate increases pond size and water volume, and also generates higher connectivity 285 between the ponds. This is because a faster debris-flux can dramatically decrease debris thickness on adjacent slopes by transporting debris into depression locations (e.g., ponds), which exposes large thin-debris-covered area that enhances ablation and facilitates meltwater production and pond expansion. To achieve a more realistic simulation, we performed numerical experiments based on modern-day surface topography, and focused on a subarea on the Baltoro Glacier ( Figure 5C) where ponds are surrounded by steep slopes (early stage ice-cliffs).

295
Similarly, we simulated the corresponding pond evolution under three different debris-flux rates (S6-S8). The bottom row in Figure 5 shows that the more intensive meltwater ponding on the glacier surface is associated with higher debris-flux rate.
These results also revealed a positive correlation between the gravitational debris-flux rate and pond size, a faster sediment flux generating larger areas of thin debris on the ice-cliffs and around the pond, thereby enhancing ice melt, and producing additional meltwater that accelerates pond expansion. Given that the differential topographic lowering caused by these ponds 300 further encourage water ponding, a positive feedback may exist between pond expansion, debris-flux, and topography that could lead to nonlinear increases in ice-mass loss on DCGs.
Note that over annual or decadal time frames, englacial drainage is a controlling mechanism on pond size, as most supraglacial ponds start to drain rapidly when they intersect englacial conduits (Benn et al., 2001). Therefore, a pond cannot infinitely expand, and our simulations only serves as a demonstration of supraglacial pond expansion before it encounters an englacial 305 drainage channel.
Supraglacial pond evolution is also controlled by the retreat of ice-cliffs, which is a form of lateral ice-mass loss. The debris is usually thinner on the pond-facing side of the ice-cliff due to the steep slope ( Figure 6A), a thin layer of debris significantly lowers the albedo when mixed with meltwater, while not providing enough insulation, which potentially leads to faster retreat if the orientation of the ice-cliff is facing the solar azimuth direction . We then compare the mean 310 ablation rate at pond boundaries (ice-cliffs as indicated by the black line in Figure 5C) as a function of debris thickness for scenarios S6, S7 and S8 ( Figures 6B, 6C, 6D). Results show the thinning of debris at pond boundaries caused by an increase in debris-flux rate, and the accelerated melt as debris gets thinner, which facilitates pond expansion. The different flux rates could be due to the different slope angles of the ice-cliff, which is a topographic control on the debris flux rate.  Bush). Note that the debris cover on the ice-cliff is thinner than adjacent areas and exhibits a relatively high moisture content (low albedo "dirty ice") due to meltwater. Simulated mean ablation rate versus debris thickness at pond boundaries are plotted for (B) A low debris flux rate scenario (S6) , (C) A moderate debris flux rate scenario (S7); and (D) A high debris flux rate scenario (S8). Note that the degree of debris thinning and ablation rate increase is regulated by the debris flux rate, which can be influenced by the topography.

315
Field studies have suggested that the overall glacier slope gradient in the ablation zone controls the formation of supraglacial ponds, such that most ponds exist over areas exhibiting relative gentle slope (Reynolds, 2000;Sakai et al., 2002;Immerzeel et al., 2014;Reid and Brock, 2014;Thompson et al., 2016). To understand and quantify the effect of glacier surface gradient variations on supraglacial pond formation and drainage, we performed three numerical experiments on glacier surfaces with three different gradients (scenarios S9: slope angle = 2º, S10: slope angle = 5º and S11: slope angle = 10º) and investigated 320 the variations in discharge and surface water storage (Figure 7, note that we define slope here as the overall gradient of the entire studied ablation zone, not the local slope). The 2º (S9) scenario represents a realistic topographic condition of the study area, and scenarios S10 and S11 represent hypothetical high slope conditions. Figure 7B shows the variations in surface water volume under different slopes, and Figure 7C depicts the increase in surface discharge as the glacier surface gets steeper. Note that the simulations start with an initial dry surface, so the filling process took about 25 days before the water volume and 325 discharge stabilize. The total amount of meltwater produced from the glacier for each scenario is adjusted to be equal to the meltwater production in S9, such that surface slope is the only factor that is responsible for the difference in discharge. The increasing discharge and decreasing water storage on a steepening glacier surface suggests that surface gradient has a significant control on supraglacial drainage efficiency. Supraglacial discharge is indicative of supraglacial pond formation because higher discharge means lower water storage on the glacier surface, which inhibits pond expansion based on our 330 assumption of no englacial filling and drainage. The lower drainage efficiency on lower gradient indicates a higher water storage capacity on the glacier surface, which is a favored condition for pond formation. These numerical results provide a quantitative explanation for the observed abundance of supraglacial ponds in ablation zones that exhibit relatively low gradients.

335
Although supraglacial ponds and surrounding ice-cliffs only make up a very small portion of the glacier surface area, they are responsible for a significant amount of ablation that is disproportionate to their size (Table 3). Previous studies have estimated that the ablation rates around ponds can be one-or two-orders of magnitude higher than that of most debris-covered areas (Sakai et al., 2000;Benn et al., 2012;Thompson et al., 2016). Our simulation over the ablation zone of the Baltoro Glacier suggests that at least 20% of ice loss can be attributed to ponding-related ablation. This is in agreement with many field observations 340 (e.g., Sakai et al., 2000Sakai et al., , 2002Reid and Brock, 2014;Anderson, 2014;Miles et al., 2016;Thompson et al., 2016;Mertes et al., 2017;Huang et al., 2018). For example, Anderson (2014) estimated that about 30% of net mass loss from the ablation zone of the Kennicott Glacier in Alaska is due to the retreat of ice-cliffs. The simulations by Miles et al. (2018) revealed that supraglacial ponds may be responsible for 1/8 of total ice loss in the Langtang valley, Nepal. Our simulated percentage based on the Baltoro Glacier fall into the range of these two studies. Supraglacial ponds can cause high-magnitude ablation 345 and ice-mass loss for two reasons: 1) ponds efficiently absorb solar energy due to their low albedo (Sakai et al., 2000;Lüthje et al., 2006a). The energy can be effectively used for ice melt because heated pond water can easily percolate through porous debris. Sakai et al. (2000) estimated that the absorbed heat per unit area for supraglacial ponds is about 7 times higher than the average for the entire debris-covered area, and our simulations indicate that the ablation rates beneath the pond surface are on average 5-6 times higher than debris-covered ice in the ablation zone.
2) The deepening of ice surface depressions due to 350 ponding can potentially cause faster debris transport around ponds, as more and more debris translocate and accumulate in the pond. Pond-facing ice-cliffs are usually covered by a thin layer of debris given their steep slopes ( Figure 6A), which lowers the surface albedo when mixed with meltwater. Threrefore, supraglacial ponds not only enhance the vertical lowering of the glacier surface, but also encourage the lateral retreat of surrounding ice-cliffs. These processes partially explain why the glacier surface tends to be more sensitive to radiative forcing, and exhibits higher non-linearity in ice-mass loss when ponds are present 355 ( Figures 4F,4G, Table 3). We also note that supraglacial ponds significantly increased the spatio-temporal variability in surface ablation ( Figures 4B,4D), given that they act as melting hotspots and vary in size and spatial density at smaller timescale than other supraglacial features. Field study has inferred a similar pattern (Huang et al., 2018), and our results provide the first modeling support to this observation.
supraglacial ponds are accounted for, which can be characterized as an acceleration in ablation rate and total ice-mass loss. This is potentially caused by a pond's status or "stage of development", such that more rapid change is likely to occur in the later portions of the ablation season, as it reaches its yearly maximum extend and depth. As such, a nonlinear response to radiative forcing may represent the beginning of a critical transition of the glacier system that signifies a critical nonlinear response that further signifies the surface ablation-pond expansion feedback. The collective feedback mechanisms associated 365 with supraglacial pond formation and evolution may be a sensitivity metric that indicates the rapid surface lowering of alpine debris-covered glaciers.

Topographic Control on Supraglacial Ponding and Drainage
Glacier surface topography plays an important role in meltwater ponding and drainage (Reynolds, 2000;Sakai et al., 2002). Our simulations show that glacier surface gradient controls pond formation by governing the magnitude of surface discharge ( Figure   370 7). A higher surface gradient promotes discharge that leads to lower water storage on the glacier surface, then less supraglacial ponds could form. Therefore, the low drainage efficiency on a lower gradient surface creates the favored condition for pond formation. This phenomenon has been observed on several DCGs, where supraglacial ponds are more abundant in zones that exhibit relatively low gradients (Reynolds, 2000;Sakai et al., 2002;Immerzeel et al., 2014;Reid and Brock, 2014). Many factor contribute to the formation of supraglacial ponds, and our simulation results partially explained this phenomenon based 375 on the gradient control on surface drainage. In addition to gradient, topographic depressions caused by differential ablation and surface water flow are also necessary for pond formation, as indicated by our simulations S3'-S5', because topographic depressions provide sinks for meltwater to accumulate. The topographic depressions are ubiquitous on DCGs given the high spatial variability in ice topography, debris thickness and ablation rate (Sakai et al., 2002;Mihalcea et al., 2008;Zhang et al., 2011;Bishop et al., 2020).

380
Collectively, our simulations suggest that the favorable surface conditions for supraglacial ponding on a DCG include relatively low surface gradient, presence of multi-scale topographic depressions, and heterogeneous debris thickness distribution, facilitated by the presence of meltwater and rapidly changing ice topography. Under these conditions, supraglacial ponds can form and expand rapidly if no englacial drainage occurs.

385
Multiple processes control supraglacial ponding (Benn and Lehmkuhl, 2000;Benn et al., 2012;Huo et al., 2021c). Our model characterizes two important positive feedback mechanisms that govern pond expansion on a DCG: 1) The pond size-ablation rate feedback. As the pond grows, its surface area and surrounding ice-cliff area get larger, they absorb more solar energy which in turn accelerates ablation rate, and the retreat of these ice-cliffs leads to further expansion of the pond.
2) The ponding-surface lowering-debris flux feedback. The decrease in debris thickness is strongly related to the development of supraglacial ponds 390 and ice cliffs (Rounce et al., 2018), and a dense spatial distribution of ponds often increases surface lowering, as they form a large number of englacial channels and the collapse of the channel roofs creates new ponds and further accelerates ablation. The steepening of ice-cliffs at the edge of ponds encourages debris thinning, which contributes to the further expansion of ponds. This process has been described in previous studies (Reid and Brock, 2014;Miles et al., 2019), and this model provides the a numerical solution to partially simulate this positive feedback. Future work is required to include more complicated processes 395 that also contribute to this feedback, for example, the relatively high surface temperature on ice-cliffs can produce significant emission of long-wave radiation that is part of the adjacent-terrain irradiance that also facilitates ice-cliff retreat, and ultimately further pond expansion (Sakai et al., 2002;Miles et al., 2016;Buri et al., 2016;Buri and Pellicciotti, 2018).

Assumptions and Limitations
Our simulations focus on the relationship between supraglacial ponding and glacier surface processes over one ablation season.

400
This model does not account for precipitation forcing and ice-flow, which are known to regulate glacier mass balance, sediment fluxes and surface topography over a multi-year time scale (Cuffey and Paterson, 2010;Benn and Evans, 2014). Therefore, this model is not suitable for characterizing ablation dynamics, ponding or pond-related englacial processes over longer time periods, as other processes must be accounted for.
This model neglects adjacent-terrain irradiance, which could cause significant ablation on steep ice-cliffs (Sakai et al., 2002;Buri et al., 2016). The fluvial debris transport that could potentially decrease debris thickness and enhance ablation is also not accounted for. In addition, this model assumes no undercut on pond boundaries that causes calving and pond expansion.
This model also neglects englacial filling and drainage, which are considered as an important controlling factors for many large supraglacial ponds (Benn and Lehmkuhl, 2000;Benn et al., 2001Benn et al., , 2017. Furthermore, this model does not account for the formation of supraglacial ponds due to the collapse of water channel roofs, which has been described in previous studies 410 (Sakai et al., 2000). The ablation caused by warm pond water outflow is also neglected and the water temperature is assumed to be constant. Collectively, although this model accounts for more processes and feedback mechanisms than any existing models we are aware of, it may still underestimate the amount of ice-mass loss and the glacier's non-linear response to radiative forcing.

415
Supraglacial ponds are melt hot-spots on debris-covered glaciers. A better understanding of the ponding process provides valuable insights into the state and fate of many glaciers under a changing climate. Currently, little quantitative knowledge is known about how debris flux and surface topography control the ponding dynamics, as multiple processes and feedback mechanisms associated with pond formation and evolution have not been accounted for in existing glacier models. In this study, we investigated supraglacial ponding and drainage in response to different debris flux and topographic conditions, and 420 quantified ponds' impact on glacier mass loss using a relatively comprehensive numerical model that characterizes surface energy balance, topographic effects, gravitational debris-flux, meltwater production, drainage and ponding. Simulation results based on the Baltoro Glacier and various hypothetical scenarios indicate that: