Controls on the relative melt rates of debris-covered glacier surfaces

Supraglacial debris covers 7% of mountain glacier area globally and generally reduces glacier surface melt. Enhanced energy absorption at ice cliffs and supraglacial ponds scattered across the debris surface leads these features to contribute disproportionately to glacier-wide ablation. However, the degree to which cliffs and ponds actually increase melt rates remains unclear, as these features have only been studied in a detailed manner for selected locations, almost exclusively in High Mountain Asia. In this study we model the surface energy balance for debris-covered ice, ice cliffs, and supraglacial ponds with a set of automatic weather station records representing the global prevalence of debris-covered glacier ice. We generate 5000 random sets of values for physical parameters using probability distributions derived from literature, which we use to investigate relative melt rates and to isolate the melt responses of debris, cliffs and ponds to the site-specific meteorological forcing. Modelled sub-debris melt rates are primarily controlled by debris thickness and thermal conductivity. At a reference thickness of 0.1 m, sub-debris melt rates vary considerably, differing by up to a factor of four between sites, mainly attributable to air temperature differences. We find that melt rates for ice cliffs are consistently 2–3× the melt rate for clean glacier ice, but this melt enhancement decays with increasing clean ice melt rates. Energy absorption at supraglacial ponds is dominated by latent heat exchange and is therefore highly sensitive to wind speed and relative humidity, but is generally less than for clean ice. Our results provide reference melt enhancement factors for melt modelling of debris-covered glacier sites, globally, while highlighting the need for direct measurement of debris-covered glacier surface characteristics, physical parameters, and local meteorological conditions at a variety of sites around the world.


Introduction
Rocky debris covers the surface of approximately 7% of mountain glacier ice globally, but over 10% for some regions and frequently over 20% for individual glaciers, where it is concentrated in lower ablation areas (Scherler et al 2018, Herreid andPellicciotti 2020). Due to difficulties in mapping debris-covered ice, as well as the melt-inhibiting effect of even moderate surface debris (e.g. Östrem 1959, Nicholson andBenn 2006), these glacier areas have been historically ignored in regional and global projections of glacier melt (e.g. Radić et al 2014, Huss and Hock 2015, Marzeion et al 2020. New debris thickness datasets offer promise for explicit representation of these areas in glacier models (Kraaijenbrink et al 2017, Compagno et al 2021, Rounce et al 2021, but the spatial variability of supraglacial debris thickness (McCarthy et al 2017, Nicholson et al 2018 and model parametric uncertainty pose considerable obstacles. Furthermore, bare ice exposures (i.e. ice cliffs) and supraglacial ponds are often scattered across the surface of debris-covered glaciers (Steiner et al 2019, Kneib et al 2021a; these features are disproportionately responsible for these glaciers' mass loss (e.g. Immerzeel et al 2014, Thompson et al 2016, Salerno et al 2017, Brun et al 2018, Mölg et al 2019 and are difficult to constrain due to their variability in space and time (Miles et al 2017a, Steiner et al 2019. Detailed assessments of ice cliff melt rates and volume losses have been performed for select glaciers, almost exclusively in High Mountain Asia (Sakai et al 1998, 2002, Han et al 2010, Reid and Brock 2014, Steiner et al 2015, Brun et al 2016, Watson et al 2017b, Anderson et al 2021, Mishra et al 2021, Stefaniak et al 2021. It is clear that ice cliff geometry is a key control on their melt rates: they experience differential melt rates based on their aspect (Sakai et al 2002, Buri and, and their steep surface complicates their radiative budget (Han et al 2010, Steiner et al 2015. In addition, ice cliffs' low albedo (Sakai et al 1998, Steiner et al 2015, Watson et al 2017b is likely to enhance their energy receipts and melt rates relative to debris-free glacier ice. Ice cliffs' steep slopes additionally increase their true surface area, often by up to 40% (Anderson et al 2021), further increasing their melt contributions.
The complex dynamics of supraglacial ponds (Röhl 2008, Benn et al 2012, Xin et al 2012, Watson et al 2016, Miles et al 2017a, 2017c, Steiner et al 2019 make it challenging to assess their role at the glacier scale (Salerno et al 2017. These features generally have a positive surface energy balance and contribute to local ablation at the glacier surface (Benn et al 2001, Röhl 2006, Salerno et al 2017, but their primary contribution to glacier ablation is the dissipation of absorbed energy to the glacier's interior when they drain (Sakai et al 2000, Benn et al 2012. Several cliff and pond studies have assessed a melt enhancement factor for these features relative to the general rate of glacier ablation beneath the surrounding surface debris (Sakai et al 2002, Immerzeel et al 2014, Brun et al 2018, Mishra et al 2021, enabling a basic representation of these features in regional glacier melt models (Kraaijenbrink et al 2017). This pragmatic approach is to-date hampered by the limited geographic representativeness of the few field observations and the use of sub-debris melt as a reference rate. To progress with understanding of mass losses for debris-covered glaciers, there is a need for an assessment of the melt rates of debris-covered glacier surfaces considering the geographic and climatic domains where surface debris is widespread. We hypothesize that the melt enhancement of ice cliffs and supraglacial ponds is highly variable due to differences in reference sub-debris melt rates, but that differences in melt rates can be related directly to local meteorological conditions. We aim to address two key questions: (a) How do the melt rates of ice cliffs and supraglacial ponds compare to sub-debris and bare ice melt rates across different climates? (b) Which physical parameters and meteorological conditions drive the heterogeneity of melt enhancement factors for debris, ice cliffs, and supraglacial ponds relative to clean ice?
To address these questions, we model the surface energy balance and associated melt rates for clean ice, debris-covered ice, ice cliffs, and supraglacial ponds using meteorological information from Automatic Weather Stations (AWSs) at debris-covered glacier sites around the world.

Methods
We compiled the most recent and complete AWS meteorological datasets available for each region. Our approach required a suite of high-quality meteorological observations encompassing much of an ablation season, which varies by site, in order to provide accurate forcing for our models and to ensure that the results are representative for real debris-covered glacier domains. Specifically, we required that stations include observations of downwelling measurements for both shortwave and longwave radiation, near-surface air temperature and relative humidity, and wind speed. These criteria enabled us to identify 13 sites covering a range of geographic and climatic conditions (figure 1; table 1). The available stations cover most regions of the Randolph Glacier Inventory (RGI, Pfeffer et al 2014) representing nearly all areas with considerable surface debris (Scherler et al 2018) and spanning the majority of the climatic envelope of debris-covered glacier sites (figure 1, right).
Using these data, our analysis aimed to calculate the melt rate for each surface type as if it were located exactly at the AWS location. We note that meteorological conditions and surface melt of course vary considerably across a single glacier, let alone regionally or globally. Our aim here is simply to ensure that the diversity of meteorological forcing where supraglacial debris is found globally is represented in our analysis, and to examine the effect of these conditions on glacier surface melt rates. We note that supraglacial ponds, ice cliffs, and debris-covered ice are indeed present at all the sites investigated.
We apply models of surface energy balance for debris-covered ice , ice cliffs , and supraglacial ponds , and also model the energy balance of bare glacier ice (Reid et al 2012) as a reference (see appendix). We model these distinct surface types for periods Figure 1. The 13 glaciers studied in this analysis cover nearly all regions of the Randolph Glacier Inventory (RGI) (Pfeffer et al 2014) which exhibit considerable surface debris (Herreid and Pellicciotti 2020) (left). These 13 sites also occupy very distinct climatic settings but are representative for the global climatic envelope of debris-covered ice (right), here evidenced by ablation season mean temperature and precipitation derived from WorldClim 2.1 (Fick and Hijmans 2017) for all pixels containing debris-covered ice. Glacier labels: 24K-24K, ARO-Arolla, CNU-Changri Nup, DJA-Djankuat, EXP-Exploradores, HAI-Hailuogou, KEN-Kennicott, LIR-Lirung, MIA-Miage, PIR-Piramide, SUL-Suldenferner, TAP-Tapado, TAS-Tasman.
where the WorldClim 2.1 (Fick and Hijmans 2017) monthly maximum air temperature is above 0 • C at the station's location, in order to compare only periods that are definitely within the local ablation season. It is worth noting that presence or absence of supraglacial debris alters lapse rates of near-surface air temperature Pellicciotti 2016, Shaw et al 2017) and that although turbulent fluxes differ over debris due to surface characteristics , similarity theory is equally applicable for clean and debris-covered ice (Nicholson and Stiperski 2020). The energy balance models for supraglacial ponds and ice cliffs have been developed based on meteorological data collected over debris, as in this study.
These four energy-balance models (debris, cliff, pond, clean) require specification of numerous physical parameters and characteristics (see appendix), which are constrained for only selected sites globally (e.g. Nicholson andBenn 2012, Rounce et al 2021). We thus prescribed conservative probability distributions for each parameter from the available literature, and determined 5000 independent, random parameter sets for each model (table 2). This Monte Carlo approach enabled us to ignore real differences in cliff, pond, and debris characteristics between sites, as we lack parameter estimates from each site, and to instead focus on the energy balance responses of each surface to distinct meteorological forcing. This allowed a comprehensive representation of all possible geometric and meteorological scenarios for each surface and site, therefore overcoming the potential bias of past studies that have investigated selected features in detail. Based on the 5000 Monte Carlo parameter sets, we determined functional relationships between debris thickness and sub-debris melt (Östrem 1959) at each site by fitting a first-order rational expression (Anderson and Anderson 2016). We also produced distributions of melt rates for ice cliffs, supraglacial ponds, and clean ice (figures 2 and 3).
We analysed the probabilistic melt rates to assess the relative melt rates compared to sub-debris melt rates at the 13 AWS locations. Specifically, we first examined the ratio of melt of cliff, pond, and clean ice surfaces to the sub-debris melt rate at each site. This ratio yielded a factor comparable to the melt 'enhancement factor' reported by past studies as a metric of the faster melt at ice cliffs compared to the surrounding debris (Sakai et al 2002, Immerzeel et al 2014, Brun et al 2018. We additionally calculated the ratio of cliff, pond, and sub-debris melt to clean glacier ice melt to provide an alternative enhancement factor unbiased by debris thickness, as in Rounce et al (2021). In these calculations we considered the full set of simulation results for each surface to produce probabilistic enhancement factors.
We then analysed the meteorological and parametric drivers of each surface's melt rate and enhancement factor variability. We first determined the mean value of key meteorological variables at each site for its model period, then related these meteorological characteristics to the variations of each surface's melt rates and enhancement factors across sites. We additionally correlated the melt rates at each site to each individual physical parameter (table 2) to assess which are most vital to constrain with field measurements, enabling us to compare the meteorological and parametric sensitivities for each surface.

Probabilistic melt rates of debris, cliffs, ponds, and clean ice
Our results show a wide variety of melt rates (figure 2, table 3). For sub-debris melt rates, the parametric variability appears much greater than the variability  between sites, leading to more than an order of magnitude of variation for each site, but mean sub-debris melt rates are below 0.02 m w.e. d −1 at all sites. Only a few sites exhibit melt rates up to 0.05 w.e. m d −1 for thin debris (24K, KEN, MIA, TAS).
Ice cliff melt rates also show considerable parametric variability, varying by a factor of two at most sites. This roughly parallels the magnitude of variation between sites, with mean melt rates falling between 0.09 (CNU, SUL) and 0.16 (DJA) m w.e. d −1 .
Ice cliffs exhibit the highest melt rates of any surface for all sites except CNU, where the distribution of ice cliff melt rates broadly overlaps with the distributions for supraglacial ponds and clean ice.
The mean supraglacial pond net energy balance (hereafter expressed as a mean melt rate) is much more variable between sites, and also varies in parametric sensitivity between sites. At PIR and TAP, the seasonal energy balance is negative for nearly all simulations, indicating overriding evaporation irrespective of parameter choice. At 24K and CNU, on the other hand, mean melt rates exceed 0.07 m w.e. d −1 . Pond mean melt rates exhibit very high parametric spread for some sites (CNU, EXP, KEN, MIA) and very tight distributions for others (LIR, TAS, SUL).
Modelled clean ice melt rates vary considerably between sites, from 0.02 m w.e. d −1 at TAP to 0.09 m d −1 at KEN and CNU. As the physicallymeaningful parameter space for this model is better constrained than for the other surfaces, modelled melt rates for most sites exhibit a fairly tight spread, with the exception of CNU. Conversely, TAS and EXP exhibit a considerably tighter spread of modelled melt rates, despite showing lower mean modelled melt rates than at many sites.

Sub-debris melt and controls
As has been demonstrated previously (Östrem 1959(Östrem , Rounce et al 2021, the presence of surface debris lowers ice melt rates relative to a clean glacier surface. Considering our hypothetical debris thickness distribution (which has a median thickness of 0.5 m), this leads to between 25× and 4× melt reduction on average. All sites show a strong reduction of melt rates with increasing debris thicknesses (figure 3), and our results indicate that debris thickness explains on average 67% of the variability in melt rates at a given site. The 95% confidence interval for our rational fit reaches 0 m d −1 for all sites between 1 and 2 m debris thickness, but modelled melt rates diverge sharply for thinner debris. Thermal conductivity (k d ) also exerts statistical control over sub-debris melt rates (mean pvalue < 0.001) and explains 8% of the melt rate variance for each site. z 0 d plays a role for sub-debris melt rates at some sites (mean p-values 0.03) but explains <1% of the melt rate variance (figure 6).
Considering a thickness of 0.1 m, our sub-debris melt rate results indicate a spread of 4× between sites based on meteorological conditions alone (figure 3). This inter-site variability dwarfs the parametric variability of melt rates for all debris parameters other than thickness. Our results indicate that the differences in sub-debris melt rates between sites are strongly correlated with differences in mean air temperature and mean downwelling shortwave radiation (figure 3 inset), which together account for 80% (50% for T a alone) of the difference in sub-debris melt rates between sites for 0.1 m debris thickness, as for clean ice (Braithwaite 1981). Notably, however, we did not find a clear association between site-specific meteorology and the sub-debris melt rate enhancement relative to a clean ice surface.
Since we determine the energy balance of debriscovered and clean ice glacier surfaces in parallel, we can leverage our results to consider, for each site, the critical debris thickness where sub-debris melt equates to the melt of a debris-free surface (Reznichenko et al 2010). In theory, at this thickness the increase in energy receipts due to low debris surface albedo is compensated by the reduced efficiency of conduction through thicker debris, but more recent work has highlighted that this effect may only occur in specific meteorological conditions (Evatt et al 2015) or not at all (Fyffe et al 2020). Our results indicate that at only six of 13 sites is the median clean ice melt rate equalled at a debris thickness on our Østrem curve, generally <4 cm but 8.5 cm at TAS (appendix). Three (ten) sites exhibited critical thicknesses for the 0.025 (0.975) confidence level. HAI, MIA and TAS exhibited critical thicknesses for the full confidence interval, while CNU, EXP and KEN did not exhibit a critical debris thickness above 0.01 m (the lower limit of our debris thickness parameter distribution) even at the 0.975 confidence level.

Melt enhancement factors for ice cliffs and supraglacial ponds
Directly comparing the model results, we determine enhancement factors normalised to sub-debris and clean-ice melt (figure 4). Cliffs and ponds exhibit high enhancement factors relative to debris-covered ice (2-30× higher) for all sites, with the exception of the pond energy balance at PIR and TAP. However, the enhancement factor spread in all instances is nearly an order of magnitude due to the variable debris thickness and its strong influence on sub-debris melt rates. Here, as we have maintained the theoretical parameter distributions between sites, the debris thicknesses do not match actual values at each site, so our enhancement factors do not entirely align with past studies leveraging, for example, geodetic mass balances (e.g. Immerzeel et al 2014, Brun et al 2018, Mishra et al 2021. However, our results are broadly in line with suggestions that cliffs and ponds enhance melt by a factor of roughly 10-20 relative to thick debris-covered areas (e.g. Sakai et al 2002, Immerzeel et al 2014, while cliffs and ponds enhance melt less where debris is generally thinner and sub-debris melt rates are higher (Brun et al 2018, Anderson et al 2021. We note that ice cliffs enhance melt more than ponds at all 13 sites considering equal planimetric areas.
Considering enhancement factors relative to a clean ice surface is more useful from a modelling perspective (Rounce et al 2021), and also reduces the apparent enhancement factor spread due to debris thickness variability. From this perspective (figure 4), it becomes clear that ice cliffs generally melt at rates 2-3 times higher than a clean glacier ice surface. Ponds, on the other hand, still show highly variable enhancement factors, and some sites are considerably more uncertain due to parametric factors. Notably, our results show that ponds are nearly as efficient as cliffs in absorbing atmospheric energy at sites in the Himalayas (LIR, CNU, 24K, HAI), but are much less efficient in all other locations. However, it is also crucial to note that not all energy absorbed by ponds will lead to glacier ablation; pond seasonal drainage and freezing must be accounted for (Miles et al 2017a.

Parametric and meteorological controls on ice cliff melt
As for sub-debris melt, differences in ice cliff melt rates between sites are most clearly associated to differences in air temperature (p = 0.024, Pearson's r = 0.62; figure 6). As for clean glacier ice, although shortwave radiation is the single largest energy flux, variations in melt rates between locations are primarily controlled by differences in turbulent fluxes, for which air temperature is a good proxy (Braithwaite 1981). Our results indicate that sky-view factors (V s,I and V s,L ) and surface slope (β i ) are controls of ice cliff melt rates at all sites (mean p-values < 0.001, ≪0.001, and 0.01, respectively). Of these, V s,L and β i help to explain the variance between model realizations (38% and 6%, respectively). Ice cliff aspect (ψ i ) varies in importance between sites but explains on average 22% of the melt rate variance for a given site ( figure 6). Investigations at individual sites have previously demonstrated that moderate aspect differences in fact dramatically alter ice cliff evolution (Sakai et al 2002, Buri and. The enhancement factor for ice cliffs relative to clean ice is also not directly associated to specific meteorological conditions, but is inversely correlated with clean ice melt rates themselves (figure 5). That is, clean ice melt rates are driven in different settings by distinct ablation regimes resulting from the interaction of air temperature and incident radiation (Hock 2005, Pellicciotti et al 2005, and all of these conditions affect ice cliff melt rates as well, with cliffs enhancing melt relative to clean glacier ice in all settings. The extra energy absorption of ice cliffs due to their low albedo represents a progressively smaller portion of the net energy balance as clean ice melt rates increase. Although the ice cliff melt enhancement appears to decrease linearly (figure 5, R 2 = 0.846), an exponential decay converging on parity with clean ice melt is more physically justified and explains slightly more of the variance in our results (R 2 = 0.881), also providing a straightforward implementation in glacier melt models.

Parametric and meteorological controls on supraglacial pond energy balance
Our results indicate that differences in supraglacial pond energy receipts are closely related to turbulent fluxes differences, both at individual sites and between sites. For any given site, only the water surface temperature estimation parameters (T ws,0 , T ws,m ) meaningfully control the supraglacial pond melt rates (mean p-values ≪0.001) explaining on average 51% and 32% of the net energy balance variance, respectively ( figure 6). Meanwhile, the variations in pond melt rates between sites are strongly associated with variations in relative humidity (p = 0.005, Pearson's r = 0.722) and wind speed (p = 0.005, Pearson's r = −0.722), and these meteorological variables also control variations in the ponds' enhancement factor relative to clean ice (figure 5). Turbulent fluxes have previously been shown to be an important energy flux for turbid supraglacial ponds in the Central Himalaya (Sakai et al 2000, and these results extend that finding to show that they are equally important globally, with contrasting outcomes in terms of pond energy balance.

Importance of cliffs and ponds
Studies have shown that ice cliffs and supraglacial ponds can considerably enhance local and glacierwide overall melt rates for debris-covered glaciers (Benn et  Our results provide several key pieces of information to this area of inquiry. First, we show that ice cliffs indeed melt at an enhanced rate for the entire climatic range of supraglacial debris, At the sites we have analysed, ice cliffs melt at 2-3 times the rate of clean glacier ice and up to 20 times the rate of debris-covered ice (depending on the debris thickness). These are higher values than previously accounted for but certainly insufficient to bring melt rates for debris-covered areas to parity with clean ice in the same setting (Brun et al 2018, Anderson et al 2021, Rounce et al 2021. We also note that the turbid supraglacial ponds typical of debris-covered glaciers absorb energy less quickly than clean glacier ice in most settings that we studied. Consequently our results support the conclusion that debris thickness patterns dominate the mass balance pattern of debriscovered glaciers, while surface features provide heterogeneity and amplify the overall melt rate in varying degrees (Anderson et al 2021).
Nonetheless, these features' ability to moderate local and glacier-scale melt underlines the importance of accounting for them carefully in glacier-and regional-scale models of mass balance and evolution. We note that melt due to cliffs and ponds is included in geodetic mass balance measurements, so glacier model calibration to these measurements implicitly represents these features (Compagno et

Priorities for new measurements
In this study we used on-debris meteorological measurements to infer the primary meteorological and parametric controls of distinct surfaces' melt rates. We sourced data exclusively from meteorological stations positioned over supraglacial debris, and our selection of sites roughly spans the climatic conditions of supraglacial debris, but for some regions no suitable data were available to us. We were unable to identify suitable on-debris AWS records for the Low Latitudes, Western Canada and US, Greenland, and Iceland, regions where debris covers a moderate portion of glacial area. Similarly, no suitable meteorological records were available for South Asia West, Central Asia, and Svalbard. Very little supraglacial debris is present for the Antarctic Periphery, Arctic Canada, the Russian Arctic, Scandinavia, and North Asia (Scherler et al 2018, Herreid andPellicciotti 2020); these regions are not represented in our study. Several other sites could be analysed for South Asia East and we opted to include four sites within this region. Although large-scale modelling studies have demonstrated the utility of distributed climate reanalysis datasets, the use of high precision on-glacier meteorological observations provides an important foundation to test and evaluate energy balance models. In this respect, the geographical distribution of sites in this study highlights the need for new observations in underrepresented areas where supraglacial debris is prevalent.
A key question relating to debris-covered glaciers relates to overall melt reduction or increase in the debris-covered area (Rounce et al 2021). At the local and glacier scale this depends heavily on debris thickness, but debris thickness remains poorly constrained globally. Rounce et al (2021) have produced the first global maps of debris thickness and sub-debris enhancement factors, opening the possibility of representing sub-debris melt within largescale models. However, the large parametric uncertainties underlying current data products may limit the utility of such datasets for projecting future melt. To reduce modelled debris thickness uncertainties, though, new methods are needed to constrain other debris properties. Furthermore, due to temporal discrepancies between the underlying global supraglacial debris area (Scherler et al 2018) and the Randolph Glacier Inventory (RGI, Pfeffer et al 2014), this dataset does not always represent debris-covered areas well (Herreid and Pellicciotti 2020). Regarding our sites, for example, the debris-covered area of Tapado Glacier is not represented within the RGI (or, therefore, the Scherler et al (2018) or Rounce et al (2021) products), although this complex debris-rich composite glacier landform is well-known (Pourrier et al 2014). Future efforts are needed to provide dateconsistent glacier extent, debris extent, and debris thickness datasets.
Finally, we note that the parametric variance is considerable for all four models, indicating that accurate knowledge of local properties is crucial for energy balance modelling. For subdebris melt, the average parametric variance (here indicated by the inter-site mean of the standard deviation of melt rates for each site) exceeds the meteorological variance (indicated by the standard deviation of site mean melt rates) by a factor of two, due to the strong control of debris thickness. For ice cliffs, the parametric variance is equal to the meteorological variance due to the importance of cliff geometry. For both clean ice and supraglacial ponds, parametric variance is half the meteorological variance. Consequently, measuring debris properties is more important for accurate sub-debris melt rate estimation than precise local meteorology. Measuring ice cliff geometry, for example leveraging aerial, terrestrial, or satellite photogrammetry to produce high-resolution digital elevation models , Watson et al 2017b, Brun et al 2018, is as important for estimating ice cliff melt as understanding local meteorology. On the other hand, an accurate representation of meteorological conditions is more important than local surface properties for estimating clean ice and supraglacial pond energy balance.

Conclusions
Our study used meteorological records collected over 13 debris-covered glacier sites around the world with four energy-balance models relevant to debris covered ice. Using Monte-Carlo simulations with meaningful distributions of physical parameters allowed us to disentangle the parametric and climatic drivers of melt rate and melt enhancement differences. Our key findings are that: • Knowledge of site-specific physical parameters (debris thickness and thermal conductivity) is crucial for estimating sub-debris melt. However, assuming equivalent parameter distributions, the spatial variability in sub-debris melt rates is best explained by air temperature. • Sub-debris melt rates are lower than clean ice and ice cliffs at all sites, and only six of thirteen sites exhibited a critical thickness. Given our empirical parameter distributions, sub-debris melt is generally one-tenth that of clean ice melt. • Ice cliff melt is 2-3 times higher than clean ice at our study sites, but this melt enhancement is small when clean ice melt rates are already high. For a given location, ice cliff melt rates are most strongly controlled by geometry: radiative view factors, slope and aspect. • The supraglacial pond net energy receipts are highly variable between sites, largely controlled by relative humidity and wind speed. The net energy receipts for supraglacial ponds are very sensitive to water surface temperature (derived parametrically in our model), but are generally lower than clean ice in the equivalent setting, and are in some cases negative (evaporation-dominated). • Climatic similarities between sites lead to distinctive regional patterns of melt rates and relative melt enhancement for each surface. • We identify that variability in supraglacial pond and ice cliff melt rates and melt enhancement factors can be effectively parametrized independent of debris (relative to clean ice melt), and are directly related to meteorological conditions, but that physical parameters can lead to considerable differences between locations at any site.
Consequently, our results are promising for effective representation of ice cliffs and supraglacial ponds in larger-scale models. In particular, provided that regional distributions of cliffs and ponds are better known in the future, our study opens efficient approaches to estimate their melt rates and contribution to glacier-scale ablation regionally. However, our results also emphasize the need to account for site-specific physical parameters, which have a strong influence on modelled melt rates.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Acknowledgments
We are supremely grateful for the provision of automated weather station data for our study sites. In particular, we would like to thank the following data providers: Patrick

Appendix
Here we provide additional details about the meteorological stations and distinct model setups, and provide a summary table of the key results for each site. Model code and results are available on request from the authors. Table 1

Energy balance models
We use established melt models to calculate the surface energy balance for a variety of glacier surface facies, always assuming that the glacier surface is snow-free (and in the ablation area) for each weather station. As the models have already been demonstrated, here we describe only the differences in model implementation from the studies in which they were established. For all models, we calculate the energy balance for a planimetric 1 m 2 patch of the surface as if directly situated beneath the meteorological station. For all model results, we report the resulting melt rate in water-equivalent (m w.e.).

Clean ice melt modelling
The melt rate for exposed, clean glacier ice is calculated using the clean-ice energy-balance model of Reid et al (2012), upon which the three other models herein have been developed. For consistency with those models, we assume that the surface is free of snow-cover and at the melting point. This is reasonable for most sites considering the relevant spatial and temporal domain (debris-covered ablation area during the ablation season), although ablation-season snowfall can modulate subsequent ablation (e.g. Fugger et al 2021). This model calculates the energy balance using in situ measurements of air temperature, relative humidity, wind speed, and incoming shortwave and longwave radiation, and requires parameters for ice albedo, roughness, and emissivity.

Sub-debris ice melt modelling
Sub-debris melt is calculated at the point scale with the model presented in detail in Steiner et al (2021). Incoming shortwave and longwave radiative fluxes are taken from measurements. Outgoing shortwave radiation is calculated with the prescribed albedo of the debris surface. Turbulent fluxes are calculated using the standard bulk approach. The model solves the flux through the debris iteratively at multiple layers of fixed thickness until surface temperature converges. This surface temperature is used for the calculation of outgoing longwave and sensible heat flux and is also used by the cliff model to derive longwave radiation emitted by the debris.

Ice cliff melt modelling
Ice cliff melt is calculated at the point scale, using measured station data of air temperature, relative humidity, wind speed, incoming short-and longwave radiation. Local shading is simulated based on slope and aspect of the cliff, assuming an even cliff shape and a constant horizon of 20 degrees for the far-topogaphy and determines the amount and timing of direct and diffuse shortwave radiation at the cliff surface. Viewing factors and emissivity values are taken from the parameter space and influence, together with the calculated surface temperature from the sub-debris ice melt model, the amount of longwave radiation emitted from the surrounding debris and received at the ice cliff surface. The model is described in detail in Buri et al (2016).

Supraglacial pond energy balance
The supraglacial pond energy balance is calculated at the point scale using the model of Miles et al (2018), which was developed specifically for meteorological measurements collected over debris , using measured station data of air temperature, relative humidity, wind speed, incoming shortwave and longwave radiation. The parametrization for water surface temperature has been established based on contemporaneous water surface and air temperature measurements . For this implementation, we do not perform a sky view correction based on nearby topography. Consequently, six parameters are needed for this model: three parameters to estimate water surface temperature, as well as an estimate of water surface roughness, broadband albedo, and emissivity.

Energy balance physical parameters
The key parameters specified for each model, including the probability distribution of values, is depicted in table 2. The probability distributions were derived based on available literature, and were used to derive 5000 randomly-sample model realizations.   Table 3.