Morphology and evolution of supraglacial hummocks on debris‐covered Himalayan glaciers

Thick supraglacial debris layers often have an undulating, hummocky topography that influences the lateral transport of debris and meltwater and provides basins for supraglacial ponds. The role of ablation and other processes associated with supraglacial debris in giving rise to this hummocky topography is poorly understood. Characterizing hummocky topography is a first step towards understanding the feedbacks driving the evolution of debris‐covered glacier surfaces and their potential impacts on mass balance, hydrology and glacier dynamics. Here we undertake a geomorphological assessment of the hummocky topography on five debris‐covered glaciers in the Everest region of the central Himalaya. We characterize supraglacial hummocks through statistical analyses of their vertical relief and horizontal geometry. Our results establish supraglacial hummocks as a distinct landform. We find that a typical hummock has an elongation ratio of 1.1:1 in the direction of ice flow, length of 214 ± 109 m and width of 192 ± 88 m. Hummocky topography has a greater amplitude across‐glacier (15.4 ± 10.9 m) compared to along the glacier flow line (12.6 ± 8.3 m). Consequently, hummock slopes are steeper in the across‐glacier direction (8.7 ± 4.3°) than in the direction of ice flow (5.6 ± 4.0°). Longer, wider and higher‐amplitude hummocks are found on larger glaciers. We postulate that directional anisotropy in the hummock topography arises because, while the pattern of differential ablation driving topography evolution is moderated by processes including the gravitational redistribution of debris across the glacier surface, it also inherits an orientation preference from the distribution of englacial debris in the underlying ice. Our morphometric data inform future efforts to model these interactions, which should account for additional factors such as the genesis of supraglacial ponds and ice cliffs and their impact on differential ablation.


| INTRODUCTION
Debris-covered glaciers are widespread throughout tectonically active mountain ranges (Herreid & Pellicciotti, 2020) where high rates of rock uplift provide sediment to glacier surfaces by landsliding (Schomacker, 2008). In High Mountain Asia, 30% of the ice mass in glacier ablation areas is covered with rock debris (Kraaijenbrink et al., 2017). Supraglacial debris occupies approximately 36% of the glacierized area in the Everest region, and the area of debris-covered ice in this region has doubled between 1962 and 2011 (Thakuri et al., 2014). Spatial variability in debris thickness results in differential ablation and heterogeneous melt production across the debris-covered area (e.g. Nicholson & Benn, 2006). Supraglacial debris layers between about 0.01 and 0.1 m thick enhance sub-debris melt by reducing albedo, whereas debris layers that exceed this critical thickness insulate the underlying ice from solar radiation and atmospheric warming to reduce sub-debris melt (Østrem, 1959;Reznichenko et al., 2010).
The debris layers on glacier tongues in the Himalaya exhibit highly undulating topography punctuated by ice cliffs and supraglacial ponds ( Figure 1) that further contribute to spatially variable and non-linear ablation regimes (Benn et al., 2012).
Recent work has quantified the contribution of differential ablation to debris-covered glacier mass loss (e.g. Brun et al., 2016;Buri et al., 2016;Fyffe et al., 2020;Miles et al., 2016Miles et al., , 2018Thompson et al., 2016;Watson et al., 2017). However, the processes that form hummocky topography and control differential ablation are not yet fully understood . At the glacier scale, a feedback exists whereby ablation promotes exhumation of englacial debris towards the terminus and the formation of medial and transverse moraines (Anderson, 2000), resulting in reduced melt towards the terminus, enhanced ablation up-glacier and declining ice flow downglacier (Anderson & Anderson, 2018;Nicholson & Benn, 2006;Rowan et al., 2015). The resulting low-angle or inverted mass balance gradient causes debris-covered glaciers to lose mass by surface lowering rather than terminus recession (Benn & Lehmkuhl, 2000;Hambrey et al., 2008;Quincey et al., 2009). Thick debris cover at the terminus and decreasing ice flow in the lower ablation area of such glaciers encourages toward-surface migration of englacial debris-rich septa and subsequently debris emergence, leading to thickening of supraglacial debris layers (Kirkbride & Deline, 2013). Steep-sided lateral moraines and thick terminal moraines inhibit the removal of debris from the glacier surface, causing the zone of debris emergence to expand up-glacier (Anderson & Anderson, 2016).
Supraglacial debris thickness in the Everest region can vary by several metres within horizontal distances of a few decimetres (Nicholson & Mertes, 2017). Differential ablation resulting from these variations in debris thickness causes non-uniform surface lowering, topographic inversion (Iwata et al., 2000) and modifies the supraglacial landscape, manifesting as a transient hummocky topography (Benn et al., 2012;Juen et al., 2014). Depressions in this topography enable the collection of meltwater and development of supraglacial ponds (Qiao et al., 2015;Thakuri et al., 2016), which enhance ablation by an order of magnitude compared to surrounding areas of debris-covered ice (Miles et al., 2018). Topographic highs may become drastically reworked where sufficiently steep slopes develop so that debris is avalanched off to expose bare ice , and the ice cliffs thus produced typically melt six times faster than the mean ablation rate for the debris-covered area (Brun et al., 2018;Reid & Brock, 2014). Lateral thermal erosion at supraglacial ponds may also form ice cliffs by steepening surrounding slopes and promoting the removal of debris (Röhl, 2008). The volume of ice contained within individual hummocks, as well as their three-dimensional geometry, limits the amount of ice-cliff backwasting (Brun et al., 2016;Miles et al., 2016).
These observations demonstrate that differential ablation is itself influenced by hummock formation and the subsequent gravitational redistribution of debris across local slopes Moore, 2017;Tonkin et al., 2016;Westoby et al., 2020). Thus, generally, hummocks may be regarded as the outcome of complex spatial dynamics operating at short, sub-glacier length scales that arise from feedbacks between ablation, topographic evolution, sediment transport and supraglacial hydrology.
Whilst concepts regarding the establishment of debris cover and its subsequent evolution are reported in the literature (e.g. Anderson, 2000;Anderson & Anderson, 2018;Iwata et al., 2000;Kirkbride & Deline, 2013), the genesis of hummocky topography has received limited study (e.g. Mölg et al., 2020;Thompson et al., 2016;Westoby et al., 2020), and how this topography is governed by the spatial dynamics outlined above and by glaciological factors is far from well understood (King et al., 2020). Classification of landforms through a full assessment of their morphometry can yield valuable information about their formation (Smith, 2014). Here, we take this approach to derive insights into the origin and evolution of hummocky topography on debris-covered glaciers. Focusing on a sample of five debriscovered glaciers in the Everest region, we undertook extensive statistical characterization of the topography of their debris-covered sections and developed a first-order classification of hummocks as a supraglacial landform. We used these results to discuss the physical factors and processes that determine the range of observed hummock forms, the glacier-scale controls on their morphometry and the conditions for their formation and evolution.

| Study area
Glaciers in the Everest region exhibit extensive hummocks across their debris-covered surfaces. We sampled the topography of five  (Watson et al., 2016). Although their proximity to each other suggests that these glaciers experience similar synoptic climatic forcing, they do vary in elevation, aspect and size (Table 1). This diversity of form presents an opportunity to characterize glacier-scale influences on hummock formation and evolution.
2.2 | Topographic data SETSM digital elevation models (DEMs) were sourced from the Polar Geospatial Centre, University of Minnesota, at 8-m grid resolution (Noh & Howat, 2015). These DEMs were assessed for geolocation errors and biases following the three-step approach of Nuth and Kääb (2011), and were co-registered to the 30-m Shuttle Radar Topography Mission (SRTM) DEM. Over stable (i.e. nonglacierized) terrain, the mean difference in the vertical direction between the SETSM and SRTM DEMs was -0.416 m with a standard deviation of 10.12 m. Following co-registration, the SETSM tiles covering the study area were mosaicked to create a single, region-wide DEM for analyses. Our measurements of hummock geometries are therefore subject to 3.8 m horizontal and 2.0 m vertical uncertainties resulting from the accuracy of the Worldview scenes that remain after correction of the SETSM DEM (Shean et al., 2016). The SETSM data were gathered in 2014 (Noh & Howat, 2015), so the topography sampled in our study represents the 2014 surface elevations of the five glaciers.  (Figure 2b). The horizontal point interval on each transect was 8 m. Bilinear interpolation was used to estimate elevation where transects did not intersect DEM grid points. The distance between transects was chosen to be small enough to ensure dense sampling of the hummock morphometry necessary for their statistical characterization but large enough to limit redundancy in the detection of peaks and troughs. The Y profiles were linearly detrended to remove the prevailing slope of the glacier tongue as required for our topographic roughness analysis.
Successive peaks and troughs on each transect were identified using the 'findpeaks' function in MATLAB® R2015b. The minimum distance for peak detection was set to 32 m (i.e. four times the horizontal point distance), because initial tests showed that peak detection was heavily influenced by low-amplitude (<1 m) fluctuations in surface elevation between adjacent points. Hummock peak spacing (peak-to-peak distance), length (trough-to-trough distance), amplitude and mean flank slope were calculated for all X transects and Y profiles (see Figure 3a for the definitions of these variables). Any uncertainty associated with our measurements arises only from relative uncertainties due to the internal consistency of the DEM, which is considered to be negligible (e.g. Quincey et al., 2007). However, the 8-m F I G U R E 2 Location of study area in the Everest region of the Himalaya. (a) Topographic data were extracted for analysis from areas highlighted in green. (b) Surface elevation data for the debris-covered tongue of Khumbu Glacier showing examples of an X transect and a Y profile in our sampling scheme. Glacier names and outlines were obtained from the Randolph Glacier Inventory version 5.0 (Pfeffer et al., 2014). Topographic imagery is from ESRI World Imagery Basemap [Colour figure can be viewed at wileyonlinelibrary.com] resolution of the DEM means that the true elevations of hummock peaks and troughs may not always be recorded, and so amplitude is considered a minimum estimate of hummock relief and peak spacing incorporates an inherent uncertainty.

| Statistical measures of hummocky topography
To classify hummocks as a landform, we examined their morphometry by compiling probability distributions of hummock peak spacing, length, amplitude and mean flank slope, and interpreting these distributions and the associated descriptive statistics. The results were used to conceptualize a typical hummock and inform the assessment of glacier-scale controls on the hummock expression. For each hummock shape variable, probability density functions (PDFs) were computed for each glacier as well as for the entire population of hummocks. Peak spacing, length, amplitude and flank slope measurements were analysed in both the X and Y directions. Skewness and kurtosis were calculated for every distribution as these statistics can inform interpretations of the conditions and controls of hummock formation (Chandler et al., 2016). Unequal variance t-tests were conducted to highlight statistically significant differences in morphometric characteristics along-glacier and across-glacier. Any directional biases identified were evaluated against anisotropy detected independently from surface roughness calculations.

| Surface roughness calculations
To complement the statistical characterization, which treats individual hummocks as countable, a surface roughness method was employed to analyse the continuous topography of each sampled area with no a priori assumptions of what constitutes a recognizable landform (Gadelmawla et al., 2002). This approach quantifies the landscape so that information on processes can be inferred over a length-scale continuum. We used it to identify the horizontal scale at which natural background roughness 'transitioned' to the topography of interest (Smith, 2014), in this case, the hummocky topography. Specifically, scaling behaviour is quantifiable by the Hurst exponent, H, in the where v is the root mean square deviation (RMSD) of detrended surface elevation measured at the horizontal scale or 'lag step', Δ, and v 0 is the RMSD at the unit scale (Shepard et al., 2001;Smith, 2014). The v and Δ data compiled for a transect constitute its variogram, and T A B L E 1 Morphology of the five sampled glaciers from the Randolph Glacier Inventory version 5.0 (Pfeffer et al., 2014) Glacier name Area (km 2 ) Elevation range (m a.s.l.) Aspect ( ) Debris-covered area (km 2 ) Slope of the sampled area ( )  Rosenburg et al. (2011). Our perception of hummocks as nonscale-invariant (thus non-self-similar) forms suggests that their topographic scaling is characterized by low H values. The range of Δ over which such scaling is found would identify their scale of occurring as distinct landforms. We expect the corresponding b (at the lower end of this range) to distinguish hummocks from background roughness that arise from short-scale supraglacial processes not intrinsic to hummock formation.
Two deviograms were computed for each glacier, one for X transects and one for Y profiles. Each deviogram was examined visually for the presence of a transition and b. Then, more precisely, sequential Bayesian change point detection was conducted to determine b by pinpointing where H changes significantly, using the algorithm of Ruggieri and Antonellis (2016). We recorded H for the data left of b as H1, and H between b and the maximum RMSD as H2. 3 | RESULTS

| Statistical characterization of hummock morphometry
In total, we recorded 1528 hummock lengths and 1116 hummock peak spacings in X, and 1088 lengths and 1067 peak spacings in Y across the five glaciers. Hummock lengths are more numerous than hummock peak spacings because some transects are so short that they sample just a single hummock, yielding a measurement of length but not spacing. The minimum, median and maximum hummock length, amplitude and flank slope measurements (Table 2) provide a first indication of a systematic difference in hummock morphometry between the two directions (recall that X and Y are defined with respect to each glacier's flow direction, while the aspects of the glaciers differ substantially). These results suggest that hummocks tend to be elongated along the glacier flow line, with a mild exception being Khumbu Glacier, and more subdued in amplitude along-glacier than across-glacier (Table 2). Consistent with these findings, hummock flanks are steeper in X than in Y. We affirm these tendencies by examining the PDFs of each morphometric variable in X and Y ( Figure 4) and their associated statistics (Table 3) and make comparisons between glaciers as well as considering the data as a whole. ing and length are greater in Y than in X, whereas the reverse tendency is found for mean hummock amplitude and flank slope (Table 3), in agreement with our preliminary interpretations (Table 2).
In considering the planform geometries of our population of hummocks, the mean peak spacing and mean length of all hummocks are by about 3 m in X than in Y ( Figure 4; Table 3). This directional excess in mean amplitude varies between glaciers from 1.4 to 5.9 m and is only significant for Jiuda Glacier. More strikingly, mean hummock flank slopes are significantly greater in X than in Y on every glacier by several degrees (Table 3). The greatest mean slope in X is found on Jiuda Glacier, and the greatest mean slope in Y is found on Ngozumpa Glacier.
Several relationships between glacier attributes (Table 1) and hummock morphometry (Table 2) are noteworthy, despite the variability between glaciers and their differing levels of anisotropy reported above. Mean peak spacing and length in both axes increase with glacier size; broadly in the order from smallest to largest glacier of Jiuda, Khumbu and Gezhongkang, Ngozumpa and Rongbuk Glaciers (Table 3). For each glacier, we discovered positive correlations in X of mean peak spacing (r 2 = 0.97) and mean length (r 2 = 0.60) with glacier width. Correlations in Y with glacier width, which increases with glacier length, are strongly positive with mean peak spacing (r 2 = 0.88) and mean length (r 2 = 0.95). The mean amplitude in Y also correlates positively with glacier width (r 2 = 0.67), although in X the correlation is much weaker (r 2 = 0.10), apart from for Jiuda Glacier. Correlation between mean flank slope and the slope of the sampled area (Table 1) was found to be weakly positive in X (r 2 = 0.42) and Y (r 2 =0.08).
Skewness and kurtosis of the frequency distribution of slopes varies across glaciers more so than the other morphometric characteristics (Table 3), hence it is difficult to observe a particular trend.
From the results presented above, we identify a typical hummock morphology. Although the hummocks sampled are not regular in shape or size, nor always symmetrical or co-aligned in orientation, the idealized form in Figure 3b summarizes key aspects of their shape: (1) their elongation along-glacier; (2) their higher amplitude acrossglacier than along-glacier; and (3)     Directional anisotropy in surface roughness is indicated by marked differences between along-glacier roughness and across-glacier roughness ( Figure 5). The signs of these differences vary between glaciers. Roughness in X exceeds roughness in Y on Ngozumpa, Rongbuk and Khumbu Glaciers, whereas the reverse is true on Gezhongkang and Jiuda Glaciers. For all glaciers, mean Max(Δ) in X exceeds mean Max(Δ) in Y by 37 m (Figure 5f). Since we find strong positive correlations between Max(Δ) and mean peak spacing in both X (r 2 = 0.84) and Y (r 2 = 0.89), this excess can be interpreted as a directional anisotropy in peak spacing.
The deviograms of the five glaciers (Figures 5a-e)

| Glacier-scale controls on hummock morphology
Across the five study glaciers, larger hummocks generally occurred on larger glaciers. This observation is based on the strong positive correlation between hummock peak spacing and length and glacier widths and lengths, which is interpreted as the result of the volume of ice and englacial debris available to form large-scale surface undulations, and the longer ice-flow histories of larger glaciers that allow more time for larger landforms to develop (Herreid & Pellicciotti, 2020;Mölg et al., 2020). The deviograms computed for the glaciers reinforce this assessment, as both the measured maximum RMSD and the lag distance at which it occurs (Table 4) increase with glacier size. Positive skewness in all of our morphometric data indicates a higher probability of the occurrence of smaller hummocks, although hummocks of greater length and peak spacing may form on larger glaciers. Hummocks were largest on Ngozumpa Glacier, which is one of the largest glaciers in the region and has the lowest elevation in our study (Table 1). Glaciers at relatively low elevations will be subject to warmer mean annual air temperatures and so are likely to experience greater ablation (Oerlemans, 1989) and thus greater rates of debris rearrangement, recorded by higher-relief hummocks, than those at higher elevations.
Hummock amplitude is also likely to be influenced by supraglacial and englacial hydrology. The skewness of across-glacier (in X) amplitudes on Gezhongkang, Rongbuk and Ngozumpa Glaciers was greater than for the other glaciers, exceeding 0.9 in each case and reflecting a greater probability of small hummocks. Jiuda Glacier has a greater probability of larger hummocks, indicated by the skewness of amplitudes in X of 0.09, and Khumbu Glacier is in the middle of the range, with skewness of amplitudes in X of 0.6. Gezhongkang Glacier is observed to have fewer supraglacial ponds than the other study glaciers. For comparison, nearby Jiuda Glacier occupies the same elevation range as Gezhongkang Glacier but exhibits a larger number of supraglacial ponds. The area of the debris-covered tongue occupied by supraglacial ponds is 2.0% for Ngozumpa and Rongbuk Glaciers (excluding the large terminal lakes) and 3.2% for Khumbu Glacier (which does not host a terminal lake) (Watson et al., 2016). While these relationships should be interpreted cautiously, our results support the suggestion of greater relief on glaciers with a greater proportion of supraglacial ponds (King et al., 2020;Mölg et al., 2020;Thompson et al., 2016;Vincent et al., 2016), but with the caveat that such relief may be dampened by the formation of large terminal lakes, which are associated with more efficient supraglacial stream networks that remove water from the glacier surface where ponds may otherwise develop (Watson et al., 2016).
The positive relationship between hummock amplitude and glacier size is supported by our surface roughness calculations. Roughness was greater along-glacier (in Y) than in X on Gezhongkang and Jiuda Glaciers, where the slopes of the sampled areas are 3.7 and 4.1 compared to <2.0 on the other three glaciers (Table 1). This observation suggests that more steeply sloping debris-covered glaciers experience slightly greater gravitational redistribution of debris along- glacier. An issue that arises from the detrending Y profiles for glacier slope is that where the slope of the sampled area is removed, the possible influence of glacier gradient on slopes across the study area is negated. However, the observation of twice as steep slopes in X would hold true without detrending in Y, as the hummock slopes would remain subdued but more tilted in the direction of the glacier slope.

| Processes and conditions for hummock formation and evolution
Here we discuss the most important processes that are expected to influence hummock morphology both along-and across-glacier. Hummock morphometry exhibits consistent directional anisotropy in several ways: the along-glacier profiles undulate on longer length scales and with lower amplitude and thus are more gently sloping than the across-glacier transects (Figure 4). The systematic difference found between the RMSD measurements in the across-and along-glacier directions confirms that hummocky topography is anisotropic ( Figure 5; Table 4), and its creation and/or evolution involve(s) directional processes (Smith, 2014). We expect that variations in hummock geometry along-glacier could result from variations in debris supply to the glacier from hillslopes, variations in ice flow, redistribution of debris on the glacier surface by gravitational processes and differential ablation. Variations in hummock geometry across-glacier are likely to be influenced by the same factors, and additionally by the erosion of debris from the inner faces of lateral moraines, which further promotes differential ablation and may explain the greater across-glacier relief we observed from our population of hummocks.
Ice flow is an obvious candidate in controlling the formation and evolution of hummocks that would account for their universal elongation along-glacier. Rockfalls from hillslope failures and sediment incorporated into snow avalanches contribute debris directly to the glacier surface. In the accumulation area, this debris is submerged and transported englacially to emerge in the ablation area.
Emergent debris and sediment deposited directly on the glacier surface in the ablation area will be transported supraglacially and become elongated along-glacier. Observations from glaciers with less mature debris layers than those in the Everest region indicate that hummocky topography initially forms from pronounced medial moraines (Mölg et al., 2020), and these longitudinal forms appear to persist as the debris layer thickens to cover the entire glacier width. These observations support the interpretation that the distribution and transport of englacial debris by ice flow is an important control in the formation of supraglacial hummocks (Figure 7a). we know of only one such field study in the Himalaya (Miles et al., 2021).
As hillslope erosion rates exhibit temporal and spatial variability that is difficult to constrain or infer from observations (e.g. Banerjee & Wani, 2018;Scherler, 2014), the rate of debris input to glaciers is often assumed to be uniform over both time and space in models of debris-covered glacier evolution (e.g. Anderson & Anderson, 2018;Rowan et al., 2015). However, pronounced spatial non-uniformity in englacial debris concentrations (Miles et al., 2021) and distinct englacial debris structures occur within such glaciers (Kirkbride & Deline, 2013). Spatially discrete debris inputs, such as rockfalls from hillslopes and the inner surfaces of lateral moraines, constitute a leading candidate reason behind such non-uniformity, and could lead to the formation of highly variable debris surface undulations. Furthermore, reduction in ice flow of the lower reaches of the glacier poses the possibility of englacial thrusting (Moore et al., 2013), which provides an additional mechanism for altering englacial debris structures and inducing variability in the spatial distribution of debris. Consequently, the spatial distribution and geometry of englacial debris septa may be fundamental to the evolution of hummocky topography. To verify these ideas, observations are required of englacially entrained debris volumes and debris concentrations within the ice and their spatial variability (e.g. Miles et al., 2021).
The glaciers in our study area are enclosed by prominent icemarginal moraines and corresponding lateral morainic troughs that prevent the transport of debris from hillslopes to the glacier tongues (Hambrey et al., 2008) and inhibit the removal of debris from the glacier surface (Anderson & Anderson, 2016). Therefore, at present, direct deposition of sediment from surrounding slopes to the lower ablation area is limited to that sourced from the collapse of the inner faces of lateral moraines, but the runout distance of such deposits is expected to be much less than the glacier width (van Woerkom et al., 2019). Given the predominant tendency for the medial moraines to have their axes oriented along-glacier, we generally anticipate that ice is uniformly thermally insulated by medial moraines, whereas acrossglacier insulation is more variable (due to alternation between successive tracts of cleaner ice and more sediment-laden ice at the position of medial moraines), and the areas of cleaner ice between debris septa will ablate more rapidly to form depressions and create gradients for the gravitational redistribution of debris (Figures 7b and c). Quite probably such a configuration helps precondition the observed anisotropy. Where medial moraines intersect underlying transverse (acrossglacier) englacial debris septa, there is potential for localized differential ablation near the intersection (Figure 7b), particularly if ice flow is reduced so that rearrangement of debris occurs mostly by gravity (Kirkbride & Deline, 2013). As multiple undulations in the debris surface also occur across-glacier, redistribution of debris away from hummock flanks must also occur in this direction.
Greater differential ablation gradients across-glacier promote the evolution of deeper surface depressions and greater slopes (Figures 7c and d), which is reflected in the form of surface hummocks, where amplitudes and mean slopes are greater across-glacier than along-glacier. In addition, the occurrence of maximum RMSD at greater horizontal intervals across-glacier than along-glacier highlights the prevalence of larger topographic variations in the across-glacier direction, which is also evident from the differences in H2 scaling behaviour in X and Y (Figures 4 and 5). Once hummocks are established, the debris surface will continue to be modified by differential ablation due to spatial disparities in debris thickness. Hence, the genesis of hummocky topography is likely the combined result of the effect of ice flow on the spatial distribution and concentration of englacially entrained debris, and subsequent feedbacks between the evolving debris topography and differential ablation . Gravitational redistribution of debris, driven by differential ablation, acts to damp surface roughness ( Figure 7e) (Capt et al., 2016). Therefore, for the observed hummocky topography to exist and persist, additional processes are required to maintain relatively high levels of surface roughness, particularly across-glacier. These processes likely comprise enhanced ablation at supraglacial ponds and ice cliffs ( Figure 7f) Miles et al., 2016Miles et al., , 2018Reid & Brock, 2014;Westoby et al., 2020). As hummocks provide the conditions for cliff and pond formation, feedbacks may exist whereby hummocks enable supraglacial pond and ice cliff development, which in turn deepens and steepens local areas on the glacier surface to sustain their topographic relief (King et al., 2020;Miles et al., 2016). Additionally, the development of glacier karst may enhance hummocky topography by the collapse of englacial conduits and the development of large supraglacial ponds (Benn et al., 2012;Watanabe et al., 1994). Therefore, supraglacial debris hummocks, and the inherent spatial non-uniformity in englacial debris volume that causes them to form and evolve, have important implications for supraglacial hydrology and the hazards posed by the evolution of supraglacial lakes on debris-covered glaciers (Benn et al., 2012;Thompson et al., 2016).
Significant variability occurred across our observed hummocks, which limits the strength of assessments made in this study. We could not, for example, convincingly identify changes in hummock geometry in the down-glacier direction along each glacier tongue; not only did cursory examination of the X and Y transects reveal no trends, it was unclear at what spatial scale(s) the X transects should be grouped and the Y profiles segmented to enable robust analysis of down-glacier changes-the present work is thus restricted to the scale of each glacier. Future investigations should evaluate and test our inferences about the glacier-scale controls on hummock morphometry with observational data from many more debris-covered glaciers, by examining their surface topography and morphology with the methods and classification developed here.
There is also an opportunity to investigate the time evolution of hummock morphology on some glaciers by using high-resolution DEMs from multiple years (e.g. Westoby et al., 2020), such as the recent 8-m High Mountain Asia DEM (Shean, 2017), which would allow identification of the processes operating as areas of highrelief expand up-glacier (King et al., 2020).

| CONCLUSIONS
We present the first detailed quantification and statistical analyses of supraglacial hummocky topography using over 1000 individual hummocks on five debris-covered glaciers in the Everest region. We Despite much geometrical variation in the population of hummocks studied, we discovered universal patterns in their form across the five glaciers. We propose that relatively large and spatially nonuniform debris concentrations within debris-covered glaciers, which may result from changes in ice flow and headwall erosion, are the dominant control on the origin of hummocky topography. The observed anisotropy in hummock geometry supports the notion that the orientation of ice flow causes the landform to originate and/or affects its evolution. Hummocky topography is further modified by gravitational redistribution of debris across the glacier surface, differential ablation at ice cliffs and supraglacial ponds, and reworking of the glacier surface by supraglacial streams. We therefore expect that improved understanding of the processes that form and maintain hummocky topography on debris-covered glaciers is important for understanding how glacier mass balance varies over time. Moreover, if hummocky topography does indeed reflect past variations in ice flow and headwall erosion, then these landforms could be considered a record of climatically driven changes in glacier behaviour over a longer period than can be investigated directly (e.g. by using satellite imagery).