Debris cover and the thinning of Kennicott Glacier, Alaska: in situ measurements, automated ice cliff delineation and distributed melt estimates

Abstract. Many glaciers are thinning rapidly beneath melt-reducing
debris cover, including Kennicott Glacier in Alaska where glacier-wide
maximum thinning also occurs under debris. This contradiction has been
explained by melt hotspots, such as ice cliffs, scattered within the debris
cover. However, melt hotspots alone cannot account for the rapid thinning at
Kennicott Glacier. We consider the significance of ice cliffs, debris, and
ice dynamics in addressing this outstanding problem. We collected abundant in situ measurements of debris thickness, sub-debris
melt, and ice cliff backwasting, allowing for extrapolation across the
debris-covered tongue (the study area and the lower 24.2 km2 of the
387 km2 glacier). A newly developed automatic ice cliff delineation
method is the first to use only optical satellite imagery. The adaptive
binary threshold method accurately estimates ice cliff coverage even where
ice cliffs are small and debris color varies. Kennicott Glacier exhibits the highest fractional area of ice cliffs (11.7 %) documented to date. Ice cliffs contribute 26 % of total melt across
the glacier tongue. Although the relative importance of ice cliffs to area-average
melt is significant, the absolute area-averaged melt is dominated by debris. At Kennicott Glacier, glacier-wide melt rates are not maximized in the zone
of maximum thinning. Declining ice discharge through time therefore explains
the rapid thinning. There is more debris-covered ice in Alaska than in any
other region on Earth. Through this study, Kennicott Glacier is the first
glacier in Alaska, and the largest glacier globally, where melt across its
debris-covered tongue has been rigorously quantified.



S1 In situ measurements
) Image with medial moraines delineated, labeled, and with elevation contours from 2013.69 % of debris thickness measurements are from medial moraines 5-7.A trail up the east side of the glacier allowed easy access to medial moraines 1, 2, and 3.More data is collected there than on the west side of the glacier.Note the light color debris in medial moraine 9, which is debris in excess of ~1 m thick.The ice extents are derived from WV imagery and aerial photos.The 'dead' ice portion of the debris-covered tongue is shown in light grey.It is defined by areas that only have measured daily mean surface velocities above 5 cm d -1 (Armstrong et al., 2016).This definition coincides with the definition of 'dead' ice from Rickman and Rosenkrans (1997).c) Location of air temperature poles and internal debris temperature profiles.The letters label the air temperature poles: Upper weather station (UWS); Middle weather station (MWS); Lower weather station (LWS).

Fig. S4
Error in debris thickness measurements plotted with sub-debris melt.Error bars for debris thickness represent changes upon repeated measurement.An exponential fit is shown here though we do not use this curve fit elsewhere.
Because melt measurements were made over different time periods we corrected each measurement to represent the full measurement period between 18 June and 16 August 2011.A degree-day factor for sub-debris melt DDF debris was therefore calculated for each sub-debris melt measurement: where T ✛ is the positive degree-days defined as the mean daily air temperature when above 0° C and Δt is one day (e.g., Hock, 2003).Air temperatures did not drop below freezing during the study period.We use hourly air temperature data from the Gates Glacier and May Creek meteorological stations to estimate the T ✛ at each site.Gates Glacier meteorological station is located at 1240 m a.s.l. and May Creek meteorological station is located 15 km to the southwest of McCarthy at 490 m a.s.l.(Figure 1, main text).The average period between measurements was 24 days, the minimum was 8 days, and the maximum was 56 days.Sub-debris melt rates presented in the main paper have had this correction applied, but ultimately this correction exceptionally small.The correction has a negligible effect on the curve fit applied to the distributed estimates of melt, but we correct for the measurement period nonetheless (Fig. S5).Measured sub-debris melt rates follow the typical shape of Østrem's curve from other glaciers at similar latitudes (Fig. S6).But sub-debris melt rates are higher for debris thicknesses greater than 15 cm, than on other debris-covered glaciers at similar latitudes.This discrepancy could be explained by the variability of air temperature within the study area: we measured sub-debris melt across a 250 meter elevation range, and most melt measurements under thick debris were derived from the lowest portion of the study area.This difference could also reflect the proportion of fine material making up the thicker debris cover downglacier (e.g., Owen et al., 2003).Juen et al. (2013) demonstrated that increasing the proportion of fine material in debris covers increases water retention, ultimately increasing debris thermal conductivity and hence melt rates.

S1.3 Mean ice cliff slopes
In establishing the uncertainty in our distributed melt estimates we need to determine the plausible uncertainty in mean ice cliff slope across the study area.In order to estimate this uncertainty we compiled all studies that present a measured mean ice cliff slope (for the population of ice cliffs analyzed in the specific study).

S1.4 Near-surface air temperature
Boundary layer conditions differ between debris-free and debris-covered glaciers (e.g., Brock et al., 2010).This is primarily due to differences in surface temperature between debris-free and debris-covered glacier surfaces.During the melt season, debris-free surface temperatures are typically near the melting point.But debris surface temperatures vary strongly with debris thickness and surface topography (e.g., Shaw et al., 2016).
We installed three air temperature poles (Fig. 2; at 513, 600, and 704 m a.s.l.) to document the interplay between debris, air temperature, and surface melt across the study area.Air temperatures were measured hourly between July 21 and August 9, using iButton thermistors in standard HOBO radiation shields.Poles were partially drilled into the ice and then debris 30 cm thick was accumulated around the pole in a 1 meter radius.Thermistors were initially placed 1.8 m above the debris surface.Thermistors were moved back to their original elevation above the debris surface in the middle of station deployment.The maximum deviation of thermistors from the original height above the surface was 14 cm.The thermistors have a measurement uncertainty of 0.5 °C based on ice bath calibration.We refer to the individual temperature poles as the Lower, Middle, and Upper weather stations (LWS, MWS, UWS).Unperturbed debris thicknesses averaged 2 cm at the UWS, 10 cm at the MWS, and 22 cm at the LWS.The UWS and MWS were placed in areas with lower surface relief while the LWS was located near the base of a topographic bulge.
Fig. S9 shows data from the LWS as well as near-surface air temperature lapse rates.Lapse rates between LWS and MWS were extremely steep, up to -74 C km -1 , with a mean lapse rate of -25 C km -1 over the measurement period (Table S2).Lapse rates were more shallow between the MWS and UWS, with a mean of -7.3 C km -1 , more typical of lapse rates above debris cover (-7.5 C km -1 , Mihalcea et al., 2006; -8.0 and -6.7 C km -1 , Brock et al., 2010; -5 to -7.8 C km -1 , Steiner and Pellicciotti,  2016).The study area is spanned by off-ice meteorological stations, which are typically used for estimating melt on glacier surfaces (Figure 1, main text).For comparison, lapse rates calculated between these off-ice meteorological stations were a more typical -5.6 C km -1 (e.g., Minder et al., 2010).
Lapse rates tended to be steeper during the afternoon and evening than the night and morning (Table S2 and Fig. S9).An observation similar to those of Fujita and Sakai (2000) and Steiner and Pellicciotti (2016) on Lirung Glacier, Nepal and Brock et al. (2010) andShaw et al. (2016) on Miage Glacier, Italy.The steepest lapse rates occurred between 5 and 7 pm.
Steep lapse rates in the afternoon and evening are likely caused by higher debris surface temperatures, and convective heating at lower elevations.The daily-lapse rate cycle was amplified on clear-sky days (Fig. S9; e.g., day-of-the-year (DOY) 203 and 204 compared to DOY 212-215).
Variable, steep lapse rates occurred on clear-sky days, while lapse rates were less variable and shallower on cooler, cloudy days.Amplified debris surface temperatures, as well as a decreased sky-view (see Steiner and Pellicciotti, 2016) at LWS may explain the extreme temperature differences between LWS and MWS.UWS and MWS are more likely to be influenced by cold-air drainage from up valley.Because debris is thicker at LWS debris surface temperatures will also tend to be higher.These extreme air temperature differences highlight the need for further study of the relationship between debris, surface temperature, and air temperature variations on debris-covered glaciers.
On sunny days, a steep air temperature lapse rate is observed above the debris.The extreme lapse rates observed in the lower portion of the study area are likely caused by feedbacks related to debris thickness.On clear days, high direct solar radiation fluxes increase debris surface temperatures for debris thicker than ~20 cm much more than for thinner debris.Debris thicknesses at LWS are in the range where debris surface temperatures rapidly increase with debris thickness.In contrast, at MWS and UWS, debris thicknesses were only 10 and 2 cm, respectively.This suggests that strong gradients in air temperature may be related to feedbacks between debris thickness and surface temperature.The steep lapse rates could also be related the large areas of exposed ice above the debris-covered tongue.The upper stations are more likely to be affected by descending katabatic winds.Air temperatures at the LWS may also be affected by the local debris-covered topography (e.g., Shaw et al., 2016).

S1.5 Debris thermal properties
In order to constrain debris thermal properties, thermistors were placed within the debris at ten sites.At each site, three to four iButton thermistors were placed in vertical profile within the debris.All profiles were paired with nearby ablation stake measurements.Temperature was measured every 30 minutes and each profile was in debris for at least one week, allowing for an assumption of a linear temperature gradient through the debris (e.g., Conway and Rasmussen, 2000; Fig. S10).
Thermistors were placed in debris between 8 and 46 cm thick, in four medial moraines composed of distinct lithological combinations.The debris cover is composed of sedimentary and volcanic rocks, including: andesite, pyroclastics, dacite, limestone, greenstone, and shale (MacKevett, 1972;MacKevett and Smith, 1972).We estimated debris conductivity using simultaneous measurements of sub-debris melt, internal debris temperature, and debris thickness (e.g., Nakawo andYoung, 1981, 1982;Kayastha et al., 2000;Mihalcea et al., 2006).The heat flux (Qm) is proportional to the effective thermal conductivity (Ke) of the debris for a given temperature gradient, which is set largely by surface temperature (Ts) due to the nearly constant 0° C summertime temperature of melting ice (Ti): where hdebris is debris thickness.The modifier effective is used to emphasize that heat is also transferred by advection of air and water (e.g., Juen et al., 2013).The sub-debris melt rate (per unit area) is related to the energy available for ablation (Qm) through: where Lf is the latent heat of fusion of ice (334 x 10 3 J kg -1 ), ρi is the density of ice (900 kg m -3 ), and a is the sub-debris melt rate (m d -1 ).We estimated mean debris surface temperatures at each site by linearly extrapolating mean internal debris temperatures.Using Eq. (S2 and S3), we then calculate Ke for each temperature profile over the entire duration of the temperature measurements (Fig. S11).Debris thickness appears to be a primary control of debris conductivity, even though the thermistor profiles were placed in debris composed of distinct lithological mixes.Perhaps the increase in thermal properties with debris thickness reflects a decrease in porosity as debris thickens.Thicker debris covers tend to have increased fine material at depth (e.g., Owen et al., 2003), which tends to reduce porosity.If little air convection occurs in the pores of the thinner, higher porosity debris then the bulk Ke would be greatly reduced due to the low thermal conductivity of air relative to rock (as described by Juen et al., 2013).

S1.6 Bare-ice melt rates extrapolated across the study area
For reference we also estimate the bare-ice melt rate through the study area for the summer of 2011, in the hypothetical case that no debris was present on the glacier.We calculate the bare-ice degree-day factor from several ablation stakes in bareice in the northeastern portion of the study area.We calculate the degree-day factor for ice (e.g., Hock, 2003) DDFice using measured bare-ice ablation ( ḃice ) and air temperatures interpolated across the glacier: where T ✛ is the positive degree-days defined as the mean daily air temperature when above 0° C and Δt is one day.Air temperatures did not drop below freezing during the study period.We use hourly air temperature data from the Gates Glacier and May Creek meteorological stations to estimate the T ✛ at each measurement location.We do not use air temperature measurements from the glacier surface because they are not necessarily representative of the ~1 km of the atmosphere above the glacier (e.g., Ohmura, 2001), our near-surface air temperature measurements do not span the full study period.Our intention is to show what bare-ice melt rate we might expect if surface debris did not perturb boundary layer air temperature.For all these reasons off-glacier meteorological stations are better for estimating hypothetical bare-ice melt rates within the study area.On-glacier positive degree-days were calculated based on off-glacier meteorological stations.The mean bare-ice degree-day factor was then used to calculate bare-ice melt rates across the study area.

S3 Additional distributed melt estimates
For all explorations in the main text and supplementary materials the summer 2011 melt rate is not maximized in the zone of maximum thinning (ZMT) (Table S4).Below we show how distributed melt estimates vary depending on the debris thickness assumed for medial moraine #9; we show more extreme uncertainty bounds for the best case shown in the main text; 3) we apply a linear backwasting rate, that varies with elevation; 4) we apply linear curve fits through the debris thickness from each medial moraine; and 5) we apply a single sigmoidal curve fit through all applicable debris thickness data.

S3.1.1 'Best case' in main text (Piecewise sigmoidal and linear curve fits)
In the main text we extrapolate debris thickness across the 9 medial moraines in a piecewise fashion (Fig. S14).For the 5 central medial moraines we apply the elevation dependent sigmoidal curve fit to all debris-designated pixels across these 5 medial moraines (Fig. 5).For medial moraines with debris thicknesses that do not show a sigmoidal debris thickness-melt relationship we apply linear curve-fit: where m is the slope, z is the elevation, and n is the y-intercept.
For medial moraine #9 on the far western portion of the glacier, where access was difficult, we apply a uniform debris thickness across the medial moraine area (Fig. S14).Some debris thickness measurements were made at the very lowest portion of the medial moraine but they are not representative of the upper part of the moraine.North and west of the study area debris thicknesses greater than 1 m were observed on medial moraine 9 in a portion of the glacier of similar character to the light colored debris presented in Fig. S2.This medial moraine also shows no ice cliffs in its upper reaches within the study area (Fig. S2).We set the mean debris thickness to 20 cm and and the interquartile range to 20 cm.
For the five central medial moraines debris thickness uncertainties are based on the interquartile range as it varies in 50 m elevation bins.For the other medial moraines, the uncertainty is based on the interquartile range for all debris thickness measurements made within the medial moraine.Note that we anecdotally saw similar debris thicknesses between medial moraine # 8 and the other central 4 medial moraines.For this reason we include it with the other central medial moraines.
The 'best case' distributed melt rate estimates are presented in Figure 12 of the main text.We also present uncertainty estimates with ±2-sigma bounds applied to backwasting rate, mean ice cliff slope, and sub-debris melt (Fig. S15).In this case for the sub-debris + ice cliff melt 99.996% of estimates are inside of the red band in Figure S15.
Our approach assumes that small-scale debris thickness variability has a negligible effect on area-averaged melt rates, despite the non-linear debris-melt rate relationship.Including small-scale variability in debris thickness cannot decrease the effective debris thickness (therefore increasing local melt) below the minimum in situ measured debris thickness in any given elevation band.We see from Figure 5 that the minimum measured debris thicknesses are much larger in the ZMT than immediately above the ZMT.Small-scale variability in debris thickness therefore cannot to lead to a peak in melt rates in the ZMT.It is also highly likely that the range of plausible effects from small-scale debris variability lie within extreme sensitivity tests as described in this and sections 4.4, and S3.1, further suggesting that this assumption is unlikely to affect our conclusions.

Fig. S14
Debris thickness measurements for each medial moraine in the study area.Shaded brown areas are off glacier such that the white in the figure indicates that the medial moraine exists at that elevation.For moraine # 9 the light area has very thick debris likely over 1 m thick.A survey found that debris was greater than 1 m just above the study area.We assign a most likely debris thickness to this moraine of 20 ±10 cm for the base case, shown in the main text.We test this assumptions' affect on the uncertainty estimates below.The middle five panels are used for the sigmoidal curve fit shown in the main paper.The red lines in the panels near the glacier margin show the linear curve fits applied to those medial moraines.Note that the lower two panels are on Root Glacier, where it joins Kennicott Glacier.Because of its accessibility from the glacier margin measurements were also taken from these medial moraines.Here, we explore reducing the assumed mean debris thickness in medial moraine #9 by 50%.The mean debris thickness applied is 10 cm instead of 20 cm as it is in the best case.The uncertainty bounds use debris thicknesses of 5 and 15 cm as compared to 10 and 30 cm uncertainties for the 'best case.'The difference between this case and the original 'best case' is small (Table S3).Here, we apply a linearly varying ice cliff backwasting rate with elevation.We apply the best fit linear curve through the median values of ice cliff backwasting rate in Fig. S7c.For this curve fit, the y-intercept is 13.94 cm d -1 and slope -0.0123 cm (d m) -1 .For the ice cliff backwasting uncertainty we perturbed the linear fit ±1 standard deviations of 2.46 cm d -1 .Except for the change in how ice cliff backwasting is represented all other parameters are the same as the 'best case.'The difference between this case and the original 'best case' is small (Table S3) and melt is not maximized in the ZMT.

S3.2 Linear curve fits applied to each medial moraine individually
For completeness, we also explore the effect of applying a linear curve fit through debris thickness measurements for each medial moraine individually.The most likely outcome for this is shown in Figure S18.
As for the case above: For medial moraine #9 on the far western portion of the glacier, where access was difficult we apply a uniform debris thickness across the medial moraines (20±10 cm; Fig. S14).Uncertainty bounds were not made in this case as linear curve fits through the interquartile range bounds produce negative debris thicknesses.The difference between this case and the original 'best case' is small (Table S3) and melt is not maximized in the ZMT.

S3.3 A single, sigmoidal debris thickness-elevation curve across all medial moraines
We also explore the application of a single, sigmoidal debris thickness-elevation curve across all nine medial moraines.Note that the estimates presented in this section of the supplementary materials produce the highest sub-debris melt rates of any of the cases explored.But the area-averaged melt rates are still not maximized within the zone of maximum thinning (ZMT), nor do they approach what bare ice melt rates that are expected to be at a similar elevation, even in the extreme uncertainty case.
We apply a sigmoidal curve fit because it best matches the pattern of debris thickness with elevation using a root mean squared error (RMSE) metric (Fig. S19).This curve fit was produced using the median debris thickness from 50 m elevation bins producing an RMSE of 0.9 cm.The sigmoidal shape is able to represent the clear changes in slope of the median debris thickness (in 50 m bins) with elevation, while a linear curve fit is unable to capture this effect.The best fit linear curve produces an RMSE to the median debris thickness in 50 m elevation bins of 1.9 cm, twice as high as the sigmoidal curve fit.

Fig
Fig. S1 Schematic of sub-debris melt rate and ice cliff backwasting rate measurements.

Fig.
Fig. S2 High-resolution orthoimage with medial moraines delineated.WorldView image from 13 July 2009 (catalog ID: 1020010008B20800) is shown in both panels.a) Image with medial moraines delineated and labeled.b) Image with medial moraines delineated, labeled, and with elevation contours from 2013.69 % of debris thickness measurements are from medial moraines 5-7.A trail up the east side of the glacier allowed easy access to medial moraines 1, 2, and 3.More data is collected there than on the west side of the glacier.Note the light color debris in medial moraine 9, which is debris in excess of ~1 m thick.

Fig
Fig. S3 Map of the general study area with in situ measurement locations.a) Locations of sub-debris melt rate measurements.Map of the study area with contours derived from 2013.b) Locations where ice cliff backwasting was measured.The ice extents are derived from WV imagery and aerial photos.The 'dead' ice portion of the debris-covered tongue is shown in light grey.It is defined by areas that only have measured daily mean surface velocities above 5 cm d -1(Armstrong et al., 2016).This definition coincides with the definition of 'dead' ice fromRickman and Rosenkrans (1997).c) Location of air temperature poles and internal debris temperature profiles.The letters label the air temperature poles: Upper weather station (UWS); Middle weather station (MWS); Lower weather station (LWS).

Fig. S5
Fig. S5 Difference between measured melt rates and expected melt rates (those corrected to the full study period) with debris thickness.The root mean squared difference between the measured and expected melt rates is 0.5 [cm d -1 ].

Fig. S7
Fig. S7 Pattern of ice cliff backwasting as it varies with elevation.a) Measured ice cliff backwasting over different time intervals.b) Ice cliff backwasting corrected to the full measurement period using the individual degree-day factor for each ice cliff using equation S1 but for backwasting instead of sub-debris melt.Note the minor effect of the correction.The mean of these data is 7.05 cm d -1 and the standard deviation is 2.46 cm d -1 .c) Expected backwasting boxplots in 50 meter elevation bins.Outliers are represented as +'s.

Fig. S8
Fig. S8 Ice cliff backwasting rate as it varies with elevation and medial moraine.a) Medial moraine #9 is the medial moraine furthest to the west in the study area and 1 is the eastern most medial moraine at the end of Root Glacier.Note that no backwasting measurements were made on medial moraine 2. There are not clear trends in backwasting rate with medial moraine.b) Ice cliff backwasting rates with error bars based on a measurement uncertainty of ±20 cm per distance measurement.

Fig. S9
Fig. S9 Air temperature and near-surface air temperatures lapse rate within the study area.a) Screen-level air temperature data from Lower weather station (LWS), Middle weather station (MWS), and incoming shortwave radiation data from the May Creek meteorological station.b) Lapse rate between LWS and MWS with May Creek shortwave radiation.c) Screen-level air temperature from MWS, Upper weather station (UWS), and incoming shortwave radiation data recorded from May Creek meteorological station.d) Lapse rate between MWS and UWS with May Creek shortwave radiation.

Fig. S10
Fig. S10 Debris internal temperatures as they vary in time.The elevation where each profile was located is shown in the title.The uncertainty of each temperature measurement is ±0.5 deg.C. Uncertainties associated with the measurement of thermistor depth are shown in Fig. S12.Debris thermal conductivity ranged from 0.525 to 2.16 W (m C) -1 with a mean of 1.06 W (m C) -1 (Figs.S11 and S12).Our estimates broadly agree with the range of previous direct measurements of Ke from debris-covered glaciers (0.85 to 2.6 W (m C) -1 ;Nakawo and Young, 1982;Conway and Rasmussen, 2000;Juen et al., 2013).The measurements are also within the range of Ke estimated from physical constants (0.47 to 1.97 W (m C) -1 ;Nicholson and Benn, 2006).

Fig. S11
Fig. S11 Thermal conductivity of the debris cover.Effective thermal conductivity (Ke) as it varies with debris thickness.Estimates are based on internal debris temperatures and sub-debris melt measurements over at least a week.

Fig. S12
Fig. S12 Thermal conductivity of the debris cover with error bars based on thermistor measurement uncertainty and debris thickness uncertainty.

Fig. S13
Fig. S13 Parameter optimization for a) adaptive binary threshold (ABT) and b) sobel edge detection (SED) ice cliff detection algorithms.On each plot, every point represents algorithm ice cliff detection error averaged across twelve manually digitized zones.The horizontal axes show true positives (i.e., automated and hand-digitized ice cliffs agree) and vertical axes show false positives (i.e., automated method predicts ice cliff where none exists).Perfect algorithm performance would plot on the origin.The coordinates for our chosen parameter sets provide an estimate of the error associated with the method.Note the differing axis limits on the horizontal and vertical axes.Markers are colored by Euclidean distance from the origin.

Fig. S15
Fig. S15 The best distributed melt rate estimate with ±2 st.dev.bounds applied to backwasting rate, mean ice cliff slope, and sub-debris melt.a) For the sub-debris melt estimates, 96.8 % of estimates are within the grey band.For subdebris + ice cliff melt 99.996% of estimates are inside of the red band.b) The percent contribution of ice cliffs to mass loss (sub-debris + ice cliff) with elevation.The red band contains the extreme range of melt contributions from ice cliffs.c) The fractional area of ice cliffs.

Fig. S16
Fig. S16 The best distributed melt rate estimate with the debris thickness medial moraine #9 halved.±1 st.dev.bounds are applied to backwasting rate, mean ice cliff slope, and sub-debris melt.The red bad contains an extreme range of sub-debris plus ice cliff melt based on compounding parameter choices such that 98.4 % of estimates lie within it (see section 2.3.1).84.1% of estimates for sub-debris melt are within the grey shaded band.b) The percent contribution of ice cliffs to mass loss (sub-debris + ice cliff) with elevation.The red band contains the extreme range of melt contributions from ice cliffs.c) The fractional area of ice cliffs.

Fig. S17
Fig. S17 The best distributed melt rate estimate with ice cliff backwasting rates that vary with elevation.±1 st.dev.bounds are applied to backwasting rate, mean ice cliff slope, and sub-debris melt.The red bad contains an extreme range of sub-debris plus ice cliff melt based on compounding parameter choices such that 98.4 % of estimates lie within it (see section 2.3.1).84.1% of estimates for sub-debris melt are within the grey shaded band.b) The percent contribution of ice cliffs to mass loss (sub-debris + ice cliff) with elevation.The red band contains the extreme range of melt contributions from ice cliffs.c) The fractional area of ice cliffs.

Fig. S18
Fig. S18 Distributed melt estimates if a linear curve fit is applied to each individual medial moraine.

Fig. S19
Fig. S19 Pattern of debris thickness with elevation across all medial moraines with sigmoidal and linear fits through the data.a) In situ debris thickness measurements with curve fits.With the sigmoidal curve fit to the 25% (grey), median (black), and 75% (grey) quartile debris thickness in 50 m bins.The points plotted are the mean debris thicknesses.The red line is the best linear fit.b) Debris thickness boxplots in 50 meter elevation bins.Outliers are represented as +'s.The median measured debris thickness shows a non-linear shape when binned with elevation.For that reason we use a debris thickness curve fit that is also non-linear.

Fig. S20
Fig. S20 The distributed melt rate estimate with a single sigmoidal curve fit applied across all medial moraines ±1 st.Dev.bounds are applied to a uniform backwasting rate, the mean ice cliff slope, and the sub-debris melt.The red bad contains an extreme range of sub-debris plus ice cliff melt based on compounding parameter choices such that 98.4 % of estimates lie within it (see section 2.3.1).84.1% of estimates for sub-debris melt are within the grey shaded band.b) The percent contribution of ice cliffs to mass loss (sub-debris + ice cliff) with elevation.The red band contains the extreme range of melt contributions from ice cliffs.c) The fractional area of ice cliffs.

Table S2 .
Near-surface air-temperature lapse rates within the study area *Lower, Middle, and Upper air temperature poles are on the glacier, May Creek and Gates Glacier stations are off-glacier.**Estimated debris surface temperatures based on clear sky days using a Fluke Infrared Thermometer.

Table S3
Optimal image processing parameters for ice cliff delineation in the July 2009 WorldView scene, as determined by a Monte Carlo method.
*DN = digital number, a 0-255 representation of relative radiance of each pixel.

Table S4 . Comparison of addition uncertainty cases
* (best case) means that in this case all parameters and assumptions are the same for the 'best case' except for the change indicated in the table.**Extremebounds are in the parentheses.