Debris cover and the thinning of Kennicott Glacier, Alaska, Part B: ice cliff delineation and distributed melt estimates

The mass balance of many valley glaciers is enhanced by the presence of ice cliffs within otherwise continuous debris cover. We assess the effect of debris and ice cliffs on the thinning of Kennicott Glacier in three companion papers. In Part A we report in situ measurements from the debris-covered tongue. Here, in Part B, we develop a method to delineate ice cliffs using high-resolution imagery and use empirical relationships from Part A to produce distributed mass balance estimates. In Part C we describe feedbacks that contribute to rapid thinning under thick debris. Ice cliffs cover 11.7% of the debris-covered tongue, the most of any glacier studied to date, and they contribute 19% of total melt. Ice cliffs contribute an increasing percentage of melt the thicker the debris cover. In the lowest 4 km of the glacier, where debris thicknesses are greater than 20 cm, ice cliffs contribute 40 % of total melt. Surface lake coverage doubled between 1957 and 2009, but lakes do not occur across the full extent of the zone of maximum glacier thinning. Despite abundant ice cliffs and expanding surface lakes, average melt rates are suppressed by debris, the pattern of which appears to reflect the debris thickness-melt relationship (or Østrem’s curve). This suggests that, in addition to melt hotspots, the decline in ice discharge from upglacier is an important contributor to the thinning of Kennicott glacier under thick debris.

Greater surface elevation changes are documented from the Kennicott debris-covered tongue than from any portion of the largely debris-free Nabesna Glacier, north of Kennicott Glacier (Supplemental Figure 1; Das et al., 2014). It is not clear why the greatest thinning of Kennicott Glacier occurs under thick debris and at rates similar to nearby debris-free glaciers. We define a zone of maximum thinning or ZMT under debris cover where the glacier thinned at an average rate greater than 1.4 m yr -1 between 1957 and 2009 ( Fig. 1; Das et al., 2014). For Kennicott Glacier, thinning rates this high only occur near the terminus under thick debris (Fig. 2). The ZMT occupies a 2-km down-glacier by 3.5-km across-glacier portion of the debriscovered tongue.
The debris-cover anomaly can be explained by constraining the individual components of the continuity equation for ice from glacier ablation zones: where H is the ice thickness, t is time, ḃ is the annual specific ablation (or loosely ice melt in the ablation zone), and Q is the column integrated ice discharge (or loosely ice dynamics) (Fig. 3). Constraining ḃ is particularly difficult due to the presence of ice cliffs, lakes, and streams on debris-covered glacier surfaces. The annual specific balance in the ablation zone can be sub-divided,ḃ =ḃ s +ḃ e +ḃ b (2) where ḃ s is the annual surface ablation,ḃ e is the annual englacial ablation, and ḃ b is the annual basal ablation rate. Surface ablation typically dominates ḃ in most glacial settings and it is not yet possible, to quantifyḃ e or ḃ b within and under debris-covered glacier tongues. Building from Eq. (1),ḃ s is negative in the ablation zone, and therefore shifts dH dt towards negative values, thinning the glacier. In the ablation zone, the sum of −dQ dx − dQ dy tends to be positive, because more ice typically flows into a fixed planview area than leaves it leading to surface uplift. This ice emergence velocity counters surface lowering caused by melt.
Two common explanations for the debris-cover anomaly follow from Eq. (1) ( Immerzeel et al., 2014;Vincent et al., 2016;Brun et al., 2018). First, it is possible that ablation within debris-covered areas is higher than we expect, causing thinning. In this case, ḃ , averaged across the glacier, is more negative than what insulated, sub-debris ablation rates suggest. Local melt hotspots such as ice cliffs, lakes, streams, and thermokarst counter the insulating effects of debris (e.g., Kirkbride, 1993;Sakai et al., 2002;Reid and Brock, 2014;Miles et al., 2018). In addition, lakes can enhance melt on debris-covered glaciers by an order of magnitude compared to proximal sub-debris melt (e.g., Immerzeel et al., 2014). Second, thinning upglacier of the thick debris cover leads to reduced ice flow to the debris-covered tongue, leading to a less positive −dQ dx − dQ dy , reduced ice emergence rates, and glacier thinning (e.g., Nye, 1960;Vincent et al., 2016;Brun et al., 2018).
On Kennicott Glacier thousands of ice cliffs are scattered within the otherwise continuous debris cover. These abundant ice cliffs will increase area-averaged melt rates and counteract the insulating effect of debris. Debris thicknesses on Kennicott Glacier are typically less than 50 cm (Part A) reducing the insulating effect of debris relative to many other previouslystudied glaciers. If melt hotspots control the location of the ZMT then we should expect their effect to be maximized in the ZMT. But if melt hotspots are not dominant then we should expect the mass balance profile to be dictated by the debristhickness melt relationship (or Østrem's curve ) downglacier (e.g., Anderson and Anderson, 2018). We therefore address: What is the mass balance profile within the debris-covered tongue of Kennicott Glacier? And do melt hotspots within debris cover maximize glacier-wide melt in the location of maximum glacier-wide thinning (or ZMT)?
To address these questions, we quantify the role of ice cliffs and sub-debris melt across the debris-covered tongue. In Part A, we show that sub-debris melt rates tend to decrease downglacier as debris thickness increases. Mean ice cliff backwasting rates, on the other hand, increase downglacier. In order to determine the mass balance pattern within the debris cover, ice cliff distribution must be quantified. We therefore present 1) a new method for remotely delineating ice cliffs using high-resolution WorldView 1 images; 2) estimates of the ice cliff distribution; and 3) combined estimates of ice cliff and sub-debris melt across the study area. In order to assess if surface lakes control the location of the ZMT we also digitized surface lake extent in 1957 and 2009.

Study glacier
Kennicott Glacier is a broadly south-southeast facing glacier located in the Wrangell Mountains, Alaska ( Fig. 1; 42 km long; 387 km 2 area). As of 2015, 20% of Kennicott Glacier was debris-covered. Below the equilibrium-line altitude at about 1500 m (Armstrong et al., 2017), 11 medial moraines (e.g., Anderson, 2000) can be identified on the glacier surface. Above 700 m elevation debris is typically less than 5 cm thick, although, locally, areas with low surface velocities tend to have thicker debris (Anderson and Anderson, 2018). Medial moraines coalesce 7 km from the terminus to form a debris mantle with ice cliffs, streams, and lakes scattered within an otherwise continuous debris cover. The first sinkhole lakes and collapse features are documented in aerial imagery from 1937 within 400 meters of the terminus (Rickman and Rosenkrans, 1997). As surface slopes have lowered in the lowest four kilometers of the glacier, thermokarst collapse features and surface lakes have expanded upglacier (Rickman and Rosenkrans, 1997).

Methods
The presence of ice cliffs and surface lakes on debris-covered glaciers partially makes distributed estimates of mass balance difficult. The extent of melt hotspots must be defined to make distributed melt estimates. A new method for the detection of ice cliffs has been developed using high-resolution digital elevation models (DEMs) with 5-meter resolution (Herreid and Pellicciotti, 2018). Despite the efforts of projects like the ArcticDEM (Porter et al., 2018), glacier coverage with high resolution DEMs is still rarer than coverage with orthoimagery. Here we introduce a new method to delineate ice cliffs 3 70 75 80 85 90 95 using solely high-resolution satellite imagery. We use this method to delineate ice cliffs on the surface of Kennicott Glacier.
Using the delineated ice cliffs and the in situ measurements described in Part A and shown in Figure 4, we estimate ice cliff and sub-debris melt across the debris-covered tongue.

Remote sensing methods
We describe an automated algorithm to delineate ice cliffs from optical satellite imagery. We use 0.5 m resolution WorldView (WV) satellite imagery acquired on 13 July 2009 (catalog ID: 1020010008B20800) to delineate ice cliffs across the study area. The 2009 WV image was the closest high-resolution image available in time to the 2011 summer field campaign. We used WV stereoimagery from 2013 to produce glacier surface DEMs at 5 m spatial resolution using the Ames Stereo Pipeline (Shean et al., 2016), which we use to represent the glacier surface in 2011.

Automated ice cliff detection methods
Our method for detecting ice cliffs relies on the observation that ice cliffs are generally darker than the debris around them.
Ice cliffs are typically coated with a thin, wet debris film, which appears darker than the adjacent, dry debris in panchromatic optical imagery (Fig. 5). In addition, steep ice cliffs are often more shaded than nearby lower-sloped debriscovered surfaces.
The workflow we outline relies on open-source Python packages, which facilitates the method's replication and improvement by other researchers. Our workflow consists of three general steps: 1) stretching the image brightness histogram to a suitable range for our ice cliff detection methods; 2) applying an ice cliff detection method; and 3) morphologically filtering of the detected ice cliffs (Fig. 5).
The ABT approach runs a moving window over the image, calculates the mean brightness value within that window, and then uses a threshold to binarize the image. Because the brightness threshold varies across the image, the ABT approach is insensitive to changes in illumination and debris color.
The SED approach estimates spatial gradients in image brightness. The Sobel operator detects high contrasts between lightcolored debris and dark-colored ice cliffs. The saturation stretch applied on the orthoimage causes dark ice cliffs to appear as featureless black regions, which the Sobel operator returns as low gradient values. We apply a brightness gradient threshold to isolate ice cliffs.
Both detection methods (ABT and SED) produce false positives from shaded, over-exposed, or textureless debris cover (SED only). The SED approach produces many false positives, which generally have a characteristic speckled appearance, and often occur in small, isolated groups. We apply morphological opening (Dougherty, 1992) to remove these isolated, regions that have been over-exposed by the saturation stretch and therefore lack texture. We remove these SED false positives by masking pixels with the maximum brightness.
To maximize correct ice cliff identification and minimize false positives we compare our ice cliff estimates to handdigitized ice cliffs from twelve 90,000 m 2 regions. The cumulative area used in the error checks was 1.8 km 2 , approximately 7.4% of the 24.2 km 2 study area (Fig. 6). There is some operator subjectivity in delineating ice cliffs from satellite imagery, especially for smaller ice cliffs. To minimize this issue, two different human operators independently delineated ice cliffs.
As these independent delineations agreed within 3% in their ice cliff area, we consider operator misidentification to be a negligible source of error. Seven parameters determine the success of these ice cliff detection methods: i-ii) the low and high end brightness values used for the saturation stretch; iii-iv) the window size and offset from mean brightness in the ABT method, v) the high-end value to use for thresholding in the SED method, and; vi-vii) the kernel sizes used in morphological filtering of the SED and ABT results. To find the best parameter set we use a Monte Carlo approach for multi-objective optimization (Yapo et al., 1998). We ran the ice cliff detection algorithm 2500 times, while varying parameters sampled from uniform distributions (Duan et al., 1992). We evaluate algorithm performance by comparing ice cliff area from the automated routine against the hand-digitized ice cliff areas. Our optimization simultaneously seeks to maximize true positive ice cliff detection, while minimizing false positives and false negatives. We manually inspect the top-performing parameter sets, ranked by Euclidean distance from the origin, which defines perfect algorithm performance ( Fig. 7; Reed et al., 2013). We chose image processing parameters slightly off the set with the smallest euclidean distance to reduce false positives (Table 1).

Distributed estimates of melt
Previous studies have estimated ice cliff backwasting rates as they vary in space using DEM-differencing, models, and in situ measurements. These approaches have shown that 1) ice cliff survivability varies strongly with aspect at lower latitude (Sakai et al., 2002;Buri and Pellicciotti, 2018); 2) ice cliff backwasting is highly sensitive to cliff slope (Reid and Brock, 2014); 3) local topography plays an important role in local ice cliff backwasting rates (Steiner et al., 2015; Part A); and 4) lakes allow for the long-term persistence of ice cliffs (Watson et al., 2017). Although, we lack detailed on-glacier meteorological data, on the Kennicott Glacier, we take advantage of a rich dataset of in-situ backwasting measurements from 60 ice cliffs (Part A). We use an alternative approach to distributed melt estimation by extrapolating our in situ measurements across the debris-covered tongue.
We extrapolate the empirical measurements from Part A (Fig. 4) across the study area. We use empirical curve fits of debris thickness as it varies with elevation, sub-debris melt as it varies with debris thickness, and ice cliff backwasting as it varies with elevation to distribute our measurements across the study area. These estimates are meant to represent the period from 18 June to 16 August 2011 (Part A).
The summer specific mass balance ḃ s is divided into contributions from sub-debris and ice cliff melt: ḃ debris anḋ b icecliff . Each 1 m pixel is designated as debris or ice cliff using the ABT ice cliff delineation method. We use the ABT method because it consistently performs better than the SED method (see Results section). For the most-likely case we We extrapolate debris thickness across the study area by applying the elevation dependent curve to all debris-designated pixels. Debris thickness hdebris varies with elevation z according to: where the fitted parameters a, b, and c have values of 25.7 cm, 571 m, and -0.24 cm, respectively.
We apply sub-debris melt-debris thickness relationship to all debris-designated pixels. The relationship between specific sub-debris melt ḃ debris and debris thickness is: where ḃ ice , the bare-ice melt rate measured near the top of the study area, and hstar have values of 5.87 cm d -1 and 8.17 cm respectively. The hyperbolic fit between debris thickness and sub-debris melt assumes that energy is transferred through the debris by conduction (see Anderson and Anderson, 2016).
We then apply the ice cliff backwasting-elevation relationship to all ice cliff pixels. Based on the results from Part A, we ignore ice cliff backwasting variation with orientation. We fit a linear-relationship between elevation z and specific horizontal ice cliff retreat:ḃ backwasting =f * z + g , (5) where the fitted parameters f and g are -0.0123 cm (m d) -1 and 13.94 cm d -1 , respectively. In Part A we did not find a significant difference between backwasting for ice cliffs with and without lakes at their base. Because the backwasting rate is measured horizontally, we apply an average dip relative to the horizontal plane (θ) to estimate the melt perpendicular to the ice cliff surface:ḃ icecliff =ḃ backwasting cos( 90−θ) (6) In the most-likely case we assume a uniform ice cliff slope (θ) for all ice cliffs of 40° based on an analysis of 2m-ArcticDEMs (Porter et al., 2018) over several seasons.
In order to estimate the mass balance with elevation we integrate the contributions of ice cliff and sub-debris ablation across 50 meter elevation bands: where b i is the mean ablation rate within the elevation band i in units of m d -1 , A debris i is the total debris-covered area within the elevation band, A icecliff i is the total ice cliff area within the elevation band, A i is the total area within the elevation band ,and dx and dy are both 1 m.
We present one most-likely distributed empirical melt estimate, which we bound with two extreme cases. For the mostlikely case the curve fits are calculated using the median of data from the 50-m elevation bins (Fig. 4). See Table 2 for the extreme parameters used for the distributed melt estimates. In the extreme cases for the debris thickness and ice cliff backwasting, curve fits were made through the 25% and 75% data points in each elevation bin ( Fig. 4; Supplemental Figure   2).

Modeled bare-ice melt rates 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 melt factor from several ablation stakes in bare-ice in the northeastern portion of the study area (Part A). We calculate the melt factor for ice (e.g., Hock, 2003) MFice 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. Gates Glacier station is located off-glacier at 1240 m elevation and May Creek station is located at an elevation of 490 m located 15 km to the southwest of the town McCarthy (Fig. 1).

Digitization of surface lakes from 1957 and 2009
In order to compare the extent of surface lakes with the ZMT we hand-digitized lakes from the 1957 and 2009 summer images. Lakes were searched for using a fixed grid to insure complete coverage of the study area. Digitized lake extents were confirmed by independent operators. Depressions on the glacier surface with exposed ice and/or ice-cut shorelines were digitized and assumed to be former lakes that subsequently drained.

Performance of ice cliff detection methods
The adaptive binary threshold (ABT) method outperforms the Sobel edge detection (SED) method. Averaged across the error checks, the ABT method correctly identifies 58% of ice cliff area, with 21% false positives. Percentages are relative to the hand-digitized ice cliff areas. The SED method yields a lower percentage of correctly identified ice cliffs (45%), but also produces fewer false positives (14%). In regions where we do not have manually digitized ice cliffs, our estimates of ice cliff area represent both true and false positives. Assuming our success rate is consistent across the glacier, we expect the ABT and SED approaches to detect 79% and 69% of the true ice cliff area, respectively.
Some systematic errors are evident, as anomalously light and dark regions of the glacier produce higher error. Regions of thin debris are especially problematic when using the SED method ( Fig. 6; see also Herreid and Pellicciotti, 2018 ). To correct for this error in the SED results, where debris is very thin, we manually removed areas with highly erroneous ice cliff detections; these only occur at higher elevations in the study area (Fig. 6). Due to its poorer performance, we do not use the SED-defined ice cliff area for the distributed mass balance estimates.

Spatial distribution of ice cliffs
The two detection methods produce broadly similar ice cliff distributions. The SED method, specifically, overestimates ice cliff area at high elevation due to the thin, dark-colored debris. Over the 24.2 km 2 study area, we estimate that ice cliffs cover 2.14 km 2 (8.8%) and 2.32 km 2 (9.7%) of ice cliff planview area using the SED and ABT methods, respectively (Fig.   8). If we apply a bias correction to the SED (30%) and ABT (20%) estimates based upon under-detection rates in manually digitized areas, the ice cliffs cover 11.4% and 11.7% of the glacier respectively. Focusing on the ABT results, which provide the most accurate estimate, we find a "humped" profile in the elevational distribution of ice cliff area. Ice cliff area peaks between 600 and 620 m a.s.l. Below this elevation, ice cliff area decreases (Fig. 8).
In total, 11.6 % of the debris-covered tongue of Kennicott glacier is occupied by ice cliffs. This is 60% more coverage by percentage than on the Changri Nup Glacier, the glacier with the second highest coverage of ice cliffs studied to date (Table   3). The Kennicott Glacier has the lowest mean debris thickness (13 cm) of glaciers with reported ice cliff coverage percentages and supports the highest percentage of ice cliffs. This implies that ice cliff coverage varies with debris thickness.
We normalized ice cliff area by glacierized area within each elevation band, which we refer to as ice cliff fractional area. Ice cliff fractional area is relatively uniform at 7-8% except for a broad peak between 500-660 m elevation within which fractional area reaches 13% at 540-560 m. The lower edge of this peak overlaps with the upper end of the ZMT (see Part C for further discussion).

Distributed estimates of melt
In Figure 9 we show the distributed estimates of melt split into sub-debris and ice cliff contributions across the study area.
While sub-debris melt decreases toward the terminus due to thickening debris, ice cliff backwasting rates increase toward the terminus due to increasing energy available for melt. When averaged across the entire study area, 81% of mass loss is expected to come from sub-debris melt and 19% from ice cliff melt. Maximum bounds for the total contribution of ice cliffs to mass loss are 12% and 28%. Figure 10 shows that the insulating effects of debris cover dominates the melt-enhancing effect of ice cliffs when averaged across elevation bands.
Modeled bare ice melt rates, which are meant to represent the hypothetical melt rate if debris were absent from the study area, increase towards lower elevations and range from -5.9 to -7 cm d -1 (Fig. 10). The dominance of decreasing sub-debris melt downglacier, due to thickening debris, results in a sharp deviation from the bare-ice melt rate near 700 m elevation (relative to the 2015 glacier surface). Elevation-band averaged sub-debris melt rates decline from 5.7 cm d -1 at the top of the study area to 1.6 cm d -1 near the terminus.
Ice cliffs produce mean elevation-band averaged melt rates of 0.33 cm d -1 at the top of the study area and 0.78 cm d -1 near the terminus. The maximum contribution of ice cliffs to band-averaged melt occurs near 500 m and has a value of 1.1 cm d -1 . Ice cliffs contribute most to mass loss in the 500 to 520 m elevation band, close to where the ice cliff fractional area also maximizes. Across all of the elevation bands, the ice cliffs between 500 and 520 m generate a maximum of 40% of the total mass loss due to ice cliffs and sub-debris melt.

Surface lake coverage
In 1957

Ice cliff detection methods
Our automated methods provide an accurate estimate of ice cliff area, though both the ABT and SED ice cliff detection methods underpredict ice cliff area, without bias corrections. These methods require that ice cliffs are dark relative to surrounding debris cover. Ice cliffs may be brighter than the surrounding debris if the ice cliffs are not covered with thin debris films or if they are strongly illuminated. Our method will therefore likely underpredict south-facing ice cliffs, although we observe many correct detections.
Future improvements to these detection methods may be achieved using more advanced image segmentation techniques (e.g., Leyk and Boesch, 2010), by utilizing image texture analysis, or by adaptively changing image processing parameters within a window moving across the image and mosaicing the results. Using multispectral imagery would also likely improve detection, although such imagery is less readily available. The detection methods presented here could be compared to the cliff delineation algorithm of Herreid and Pellicciotti (2018) using existing high-resolution DEMs.

Distributed estimates of melt
On Kennicott Glacier, ice cliffs most likely contribute 19% of volume loss of the debris-covered tongue. This percentage is more than twice the percentages reported from other glaciers with mean debris thicknesses less than 50 cm (Table 3) likely due to the high fractional coverage of ice cliffs on the Kennicott Glacier. For glaciers with mean debris thicknesses much larger than 50 cm, ice cliff contributions are larger than 19% and reach 40%. These high ice cliff contributions occur despite much lower ice cliff coverage when compared to Kennicott Glacier (Table 3).
Ice cliffs tend to contribute a higher fraction of mass loss as debris thickness increases. This trend is visible on Kennicott Glacier as debris thickens toward the terminus (Fig. 10). This relationship also appears to hold when considering debriscovered glaciers from different regions ( Table 3). As debris thickens the contribution of ice cliff melt also tends to increase.
This appears to occur even though the fractional coverage of ice cliffs tends to decrease as mean debris thicknesses increase.
Ice cliffs do not counteract the insulating effects of debris on Kennicott Glacier (Figs. 10 and 11). Thin debris leads to melt rates closer to bare-ice melt rates. Kennicott Glacier also has the highest fractional coverage of ice cliffs, relative to other studied glaciers, which serves to increase melt rates (Table 3). Despite this, ice cliffs do not compensate for the insulating effects of debris. This strongly suggests that the presence of ice cliffs is unlikely to completely counter the insulating effects of debris on other glaciers with thicker debris or lower ice cliff coverage.

Do melt hotspots maximize melt in the ZMT?
In order to discuss the relationship between debris, mass balance, and ice dynamics to the thinning of Kennicott Glacier we must address the assumptions inherent in using the melt pattern from the summer of 2011 to represent the melt pattern over the last 60 years. At least from 2011, the surface mass balance within the zone of maximum thinning (ZMT) is strongly suppressed by thick debris cover. We assess what changes to our melt estimates would be required to produce the highest glacier-wide melt rates within the ZMT.
Debris cover: Debris in the ZMT would have to decrease to 10% of its current thickness to produce maximum glacier-wide melt rates. In Part A, we noted that most of our debris thickness measurements were derived from the top of ice cliffs and topographic highs. Because debris tends to concentrate in topographic lows our debris thickness measurements may be biased toward thinner debris, making the required reduction in debris thickness even more extreme.
The ZMT has been continuously debris covered since at least 1957 (Fig. 12). The presence of thermokarst features, ice cliffs, and lakes in the lower 4 km of the glacier surface imply that debris greater than 5 cm was present in the ZMT between 1957 and 2009 as is the case for other debris-covered glaciers with abundant surface ponds and thermokarst features (e.g., Sakai et al., 2000;Benn et al., 2001;Wessels et al., 2002;Thompson et al., 2016). Debris cover expanded upglacier by 1.6 km between 1957 and 2009, but upglacier from the ZMT (Figs. 2 and 12). This suggests that sub-debris melt did not control the location of the ZMT.
Ice cliffs: In order for ice cliffs, in the ZMT, to enhance melt and produce maximum glacier-wide melt rates in the ZMT, backwasting rates would need to be 7.5 times higher than those measured in the summer of 2011. This required backwasting rate is well beyond potential biases introduced due to the summer of 2011 having anomalously low air temperatures. Our backwasting estimates are based on repeated measurements at a single location at the top of each ice cliff. But as described in Part A maximum backwasting rates across each ice cliff are likely to occur near the top (Buri et al., 2016). Applying our measurements across single ice cliffs or the entire ice cliff population may therefore overestimate ice cliff melt. The hypothetical backwasting rates required to maximize melt in the ZMT are therefore unreasonable; a compilation of previously published backwasting rates shown in Part A and Table 3 support this.
In order for ice cliffs in the ZMT to compensate for the insulating effects of debris and enhance melt in the ZMT beyond bare ice melt rates, ice cliff area would need to increase from 11.7% to 90% of the glacier surface. While aerial photos from the summers of 1957 and 1978 reveal an anecdotal increase in exposed ice (Fig. 12), at the upper end of the ZMT, far less than 90% of the glacier surface was occupied by ice or ice cliffs. This suggests that ice cliff melt did not solely control the location of the ZMT.
Surface lakes: Figure 11 shows that in both 1957 and 2009 lakes are most abundant near the terminus and immediately downglacier from the ZMT. Surface lakes and thermokarst depressions do not coincide with the full extent of the ZMT and are most notably absent in its upper reaches ( Fig. 11; Rickman and Rosenkrans,1997). Mass loss beneath lakes surfaces in the ZMT would have to be 190-times the local sub-debris melt rate for sub-aqueous melt to compensate for the insulating effects of debris in the ZMT. This further suggests that sub-aqueous melt under lake surfaces is unlikely to compensate for the insulating effect of thick debris.
But the mottled nature of the thinning pattern in the ZMT highlights the melt-enhancing effects of melt hotspots. The patches of most rapid thinning occur near the 2015 glacier margin which are bordered by alluvial-bedded streams and ice cliffs. While the local effect of these processes is apparent in the thinning map, it is unclear how thermokarst, surface lakes, and ice-marginal streams enhance melt outside of the hotspots themselves.

Østrem's curve expressed in the mass balance profile
Debris tends to thicken towards debris-covered glacier termini as is the case for Kennicott Glacier (e.g., Anderson and Anderson, 2018). This leads to the expectation that sub-debris melt rates will decline towards the terminus. At least on Kennicott Glacier it appears that melt hotspots do not compensate for the melt-insulating effects of thick debris. Figure 10 shows this Østrem's curve like pattern. The low melt, low melt gradient portion of the mass balance profile extends further than the high melt high-melt gradient portion of the mass balance profile. Similar theoretical expressions of this mass balance profile can be seen in the numerical simulations presented in Anderson and Anderson (2016).

Conclusions
Our new ice cliff delineation method using high-resolution satellite imagery reveals that Kennicott Glacier supports the highest percentage of ice cliffs (11.7%) of any debris-covered glacier studied to date. Ice cliffs within the debris-covered portion of Kennicott Glacier partly counters the insulating effect of debris. Approximately 19% of melt in the study area is attributable to ice cliffs. In the lowest 4 kilometers of Kennicott Glacier, debris is thick, ice cliff coverage is low, but ice cliffs still contribute up to 40 % of mass loss in this area. Ice cliffs contribute a larger percentage of mass loss within thicker debris covers, a trend that can be seen across the Kennicott Glacier and for other studied glaciers.
The upglacier and almost doubled their areal coverage in the debris-covered portion of the glacier (see Rickman and Rosenkrans, 1997). But lakes are not extensive enough to solely control the location of maximum thinning.
It appears that melt hotspots (ice cliffs and surface lakes) are unable to compensate for the insulating effects of thick debris cover on Kennicott Glacier. This suggests that, in addition to melt hotspots, ice dynamics and the decline in ice discharge from upglacier has played an important role in the thinning of the glacier under thick debris. The mass balance profile within the debris covered portion of the glacier appears to follow Østrem's curve (the debris thickness-melt relationship). In Part C we explore feedbacks that help define the rapid thinning of Kennicott Glacier under thick debris.

Data availability
Datasets and results are available upon request.    Note the substantial thinning that has occurred upglacier from the continuous debris cover.  Table 2 for curve fit parameters). The double-headed arrow represents the zone of maximum thinning (ZMT).