Englacial Drainage Drives Positive Feedback Depression Growth on the Debris‐Covered Ngozumpa Glacier, Nepal

The development of hummocky topography is a poorly understood aspect of down‐wasting on debris‐covered glaciers that is often attributed to variable debris thickness. Thousands of enclosed depressions pit the hummocky topography. To better understand depression growth, we examined the size distribution and geometry of depressions on the Ngozumpa Glacier, in the Everest Region of Nepal. The depressions exhibited a power‐law size distribution, fractal perimeters, and power‐law depth‐area scaling, which suggest positive feedback growth. With a simple model, we showed that positive feedback growth produces similar power‐law size distributions. Based on these findings, we propose a “sinkhole” hypothesis for the development of depressions. Drainage into englacial sink points removes debris from the depressions and inhibits ponds from overflowing, thereby enabling positive feedback growth via incision, increased sub‐debris melt rates, and ice cliff retreat. By facilitating sustained depression growth, englacial drainage preconditions the ablation zone for the rapid growth of glacial lakes.

The proliferation of melt hotspots coincides with the development of complex, hummocky topography on debris-covered glaciers (Bartlett et al., 2021;Mölg et al., 2020;Westoby et al., 2020). Understanding feedbacks between the development of hummocky topography and melt hotspots represents a key challenge to predicting melt rates from debris-covered glaciers. The development of hummocky topography is often attributed to variable debris thickness (Benn & Evans, 2014;Sharp, 1949). Topographic evolution models (Moore, 2021) show that variable debris thickness produces moderate topographic relief, but hillslope diffusion homogenizes debris thickness and causes relief to decay (Fowler & Mayer, 2017). This contrasts with observations that show relief increased as hummocky topography developed (King et al., 2020) and debris thickness remained highly heterogeneous (Nicholson et al., 2018;Rounce et al., 2018). These disparate findings demonstrate the need to reframe how we think about the development of hummocky topography.
Hummocky topography is usually associated with its positive expression, that is, the hummocks (Bartlett et al., 2021). However, it is the negative topographic expression-the depressions-that directly represent ablation. An enigmatic quality of the depressions is that many do not permanently contain ponds (Röhl, 2008;E. S. Miles et al., 2016E. S. Miles et al., , 2017Thompson et al., 2016;Watson et al., 2016;Benn et al., 2017). Karst landscapes also develop an abundance of dry depressions (sinkholes) that form when water dissolves soluble bedrock and evacuates sediment as it drains into fractures. Much like karst, meltwater drainage into crevasses and permeable seams of englacial debris forms caves within debris-covered glaciers (Gulley & Benn, 2007). Surveys showed that caves often originate or terminate within depressions Gulley, Benn, Screaton, & Martin, 2009). Based on these observations, we hypothesize that depressions form similarly to sinkholes in karst. Although hummocky topography on debris-covered glaciers has been compared to karst topography (Clayton, 1964), these similarities have not been quantitatively examined.
Analyses of size distributions and geometric scaling relationships have proven useful for revealing the origins of depressions in karst landscapes (Pardo-Igúzquiza et al., 2018;Raines, 2020). Karst bedrock dissolution is a positive feedback process that can produce power-law size distributions of sinkholes (Pardo-Igúzquiza et al., 2020;Yizhaq et al., 2017). Meltwater drainage, pond development, and ice cliffs initiate similar positive feedbacks on debris-covered glaciers (Anderson et al., 2021;Huo et al., 2021;Watson et al., 2016) that could drive depression growth. In contrast, negative feedback produces geomorphic feature size distributions with characteristic scales (Dong et al., 2019;Evans, 2003). Models suggest that sub-aqueous, sub-debris melt rates are >20 times slower than melt rates along the pond margins (E. S. Miles et al., 2016). Thus, without englacial drainage, we hypothesize that most depressions will have ponds, that debris accumulation and slow sub-aqueous melt rates in ponds will limit deepening, and pond overflow will limit depression areas.
To better understand depression development, we examined the size distribution and geometric scaling of depressions on the Ngozumpa Glacier, in the Everest Region of Nepal (Figure 1). Using a geometric model, we explored how growth rates, depression densities, and mergers influence the size distribution of depressions. Collectively, these analyses provide new insights to the development of depressions and hummocky topography on debris-covered glaciers.

Measuring Topographic Depressions
Our analyses of depressions on the Ngozumpa Glacier are based on a digital elevation model (DEM) generated using Structure-from-Motion photogrammetry. Similar to Wigmore and Mark (2017), we conducted an aerial drone survey of the lower ∼6 km of the Ngozumpa Glacier in December 2019 using a DJI Mavic 2 Pro with a 24 mm camera. Fifteen survey flights of the ablation zone were launched from the west lateral moraine, comprising 3,791 nadir images, each recorded with the image capture latitude, longitude, and elevation. Flights were conducted at an elevation 50 m above the moraine take-off sites. Additionally, we recorded the locations of 12 ground control points (Table S1 in Supporting Information S1) on the glacier surface using a pair of Emlid Reach dGPS receivers to aid model construction.
Following the 2017 USGS Agisoft PhotoScan workflow ( Figure S1 in Supporting Information S1), we constructed a high resolution (<0.1 m/pixel) DEM of the lower ablation zone of the Ngozumpa Glacier using Agisoft Photo-Scan v1.0. With QGIS 3.22, the high resolution DEM was rescaled to 1-m resolution and reprojected to UTM Zone 45N.
We applied the Wang and Liu (2006) sink-filling algorithm to identify topographic depressions. This method does not identify sub-depressions that may exist within a depression. The depressions were converted to shapefiles to calculate surface areas, perimeters, and depths. Since ponds may influence the development of depressions, we identified them within the DEM. We classified connected areas ≥10 m 2 with surface slopes <2° as ponds (Text S2 in Supporting Information S1).

The Size Distribution of Depressions
We examined the complementary cumulative distribution function (CCDF) of depression areas. Distinguishing power-law distributions from other heavy-tailed distributions is non-trivial (Avnir et al., 1998;Clauset et al., 2009;Newman, 2005;Voitalov et al., 2019) and has led to controversy (Corral & González, 2019;Dong et al., 2019;Raines, 2020;Stumpf & Porter, 2012). Crucially, commonly applied hypothesis testing methods (e.g., R 2 , KS-tests) are statistically invalid for identifying power-law distributions (Newman, 2005;Voitalov et al., 2019). Voitalov et al. (2019) presented an appropriate method to identify power laws from empirical data. In this approach, the power-law CCDF takes the form where ( ) represents probability of a feature of size A and γ is the power-law tail exponent. A slowly varying function, ℓ(A), accounts for the stochasticity of natural phenomena and reclassifies the distribution as a regularly varying distribution. This reclassification makes it possible to calculate the tail exponent using three independent estimators. Voitalov et al. (2019) showed that a distribution can be conservatively considered power-law if the extreme value index (ξ) of the Hill, Moments, and Kernel estimators converge to values greater than 1/4. Power-laws denote scale invariance, thus true power-law scaling should span the full size range of the features.
The value of the power-law tail exponent determines the number of defined statistical moments 〈A m 〉 (Newman, 2005). Tail exponents γ ≤ 2 denote distributions with undefined moments. For all tail exponents, the statistical moment 〈A m 〉 of the distribution exists if m < γ − 1.

Depression Geometric Scaling Relationships
If the geometric properties of depressions, such as perimeter, area, and depth, follow tight scaling relationships with each other, then it can be presumed that the depressions evolve along these relationships. Otherwise, the relationships would be quickly erased with scatter. Therefore, to quantify how depression morphology changes as the depressions grow, we examine perimeter-area and depth-area scaling relationships. Geometric scaling exponents were calculated by fitting a least-squares regression to the data in logarithmic space (Korvin, 1992;Mandelbrot, 1982;Turcotte, 1997).
Perimeter-area scaling describes how the perimeter P of a shape increases as area A increases. Shapes that obey the scaling relationship P ∝ A D/2 are self-similar, where D is the perimeter dimension. Euclidean shapes have dimension D = 1, whereas fractal shapes have perimeter dimensions 1 < D < 2 (Mandelbrot, 1967).
Similarly, the depth-area scaling relationship describes how the depression depths d scale with area. Euclidean volumes (e.g., hemispheres) follow the depth-area scaling relationship d ∝ A 1/2 . Self-similar, non-Euclidean volumes follow a more general power law, d ∝ A δ .

The Distribution and Geometry of Depressions
The sink-filling algorithm identified 31,324 depressions with areas ranging 10 0 -10 5 m 2 . We assumed that many small depressions result from debris roughness instead of depressions in the ice surface. To minimize the influence of debris surface roughness on the depression scaling relationships, we excluded depressions smaller than 10 m 2 from analysis, reducing the total number to 3,073. Only ∼4.3% of these depressions contained ponds (Text S3 in Supporting Information S1).
The distribution of depression areas ( Figure 2a) was power-law, with tail exponent γ ≈ 1.68 (γ Hill = 1.667, γ moments = 1.667, γ kernel = 1.714). The extreme value indices of three consistent estimators (ξ Hill = 1.5, ξ moments = 1.499, ξ kernel = 1.403) of the power-law tail exponent converged to similar values >1/4, indicating robust power-law scaling. The power-law scaling spanned the full range of depression areas (10 1 -10 5 m 2 ). Empirical power laws typically deviate slightly from the expected scaling for the largest features because of the finite sample size (Newman, 2005;Voitalov et al., 2019). The glacier size limits the number of large depressions, causing the power-law tail to "roll-off" for the largest depressions.
The depressions followed power-law perimeter-area scaling relationship. The perimeter-area scaling (Figure 2b) showed that the depressions have fractal perimeters, D = 1.16 (R 2 = 0.965), spanning four orders of magnitude in area (10 1 -10 5 m 2 ). The presence of ponds within the depressions did not significantly influence the fractal perimeter.
Depressions deepened (Figure 2c) by more than two orders of magnitude (∼10 −1 -10 1 m) as area increased, with depth-area scaling exponent δ = 0.49 (R 2 = 0.684). Substantial depth variability existed for the smallest depressions, but the scatter diminished as depression area increased. Most depressions with ponds were shallower than predicted by the scaling relationship. This difference is most apparent for the smallest depressions with ponds.

Model Setup
Both positive feedback growth and critical phenomena produce power-law distributions (Newman, 2005). Percolation models describe critical phenomena (Newman, 2005). This framework helped explain the power-law size distribution and fractal scaling of Earth's terrestrial lakes (Cael & Seekell, 2016). Percolation models predict a size distribution power-law tail exponent of γ ≈ 2.05 (Isichenko, 1992) and fractal dimension of D ≈ 1.33. These values are inconsistent with the size distribution tail exponent (γ ≈ 1.68) and fractal dimension (D = 1.16) of depressions on the Ngozumpa Glacier. Positive feedback growth places fewer constraints on tail exponents and fractal dimensions (Newman, 2005).
We developed a geometric model to test the hypothesis that the power-law size distribution of depressions arose through positive feedback growth. Although our model simulates no physical processes, it directly examines how feedback and mergers between depressions influence the size distributions.
Depressions are initiated as circles placed randomly within a square model space, with side length ℓ = 1,000 m. Initial depression areas were selected from a log-normal distribution (ln μ = ln 0.1 m 2 , ln σ = ln 0.5 m 2 ). Functionally, the log-normal distribution provides a broad range of initial depression areas while also ensuring that extremely small values remain uncommon (Strahler, 1954) and approximates a "random" topographic surface better than other distributions (Huang & Qin, 2014;Landy et al., 2020).
The area, A, of each depression increased following the growth rate function: where r is the growth rate exponent and k is a proportionality constant. We discretized the growth equation and explicitly approximated solutions using the fourth order Runge-Kutta finite difference method, with time steps Δt = 0.01 and k = 0.25. We simulated four growth rate exponents ( Figure 3a): constant radial growth (r = 0.5), exponential growth (r = 1), and two cases of hyperbolic growth (r = 1.25, 1.5).
Mergers occur when two or more depressions join to become a single depression. Thus, mergers potentially influence the size distributions. To examine this, we varied the initial number of depressions n = [500, 750, 1,000, 1,250] and ran simulations with and without mergers. For simulations with mergers, increasing n increases the likelihood of mergers as the depressions grow. When a merger occurs, the newly merged depression maintains its irregular shape and its area increases following Equation 2 by expanding its perimeter.
For each parameter set, we combined the results from 11 simulation runs to produce distributions comprising ∼10 3 -10 4 simulated depressions. Running 11 simulations ensured the number of simulated depressions was approximately equivalent to the number of depressions on the Ngozumpa Glacier (∼3,000). The distributions from each set of 11 simulations were then tested for power-law tails (Section 2.2). Each simulation was stopped when the largest depression exceeded 10 5 m 2 , corresponding to the scale of the largest depressions on the Ngozumpa Glacier.

Model Results
Size distributions were deemed similar to the depressions on the Ngozumpa Glacier if they exhibited a power-law tail exponent γ ≤ 2 and scaling over at least four orders of magnitude in area. Figure 3c visualizes a simulation with mergers where n = 500. Representative distributions from simulations (n = 500) without mergers ( Figure 3b) and with mergers ( Figure 3d) illustrate the influence of the growth rate exponent on the resulting CCDFs. Results from all simulations are summarized in Text S5 and Table S2 in Supporting Information S1.
Of the parameters varied, the growth rate exponent exerted the greatest influence on the distributions. When r ≤ 1 without mergers, no distributions met the criteria to be classified as power-law, nor did the depression areas span the range observed on the Ngozumpa Glacier. Only hyperbolic growth rates (r = 1.25, 1.5) produced power-law size distributions similar to the Ngozumpa Glacier. These distributions spanned areas ranging ∼10 1 -10 5 m 2 and produced tail exponents γ ≤ 2.
Allowing mergers between depressions did not significantly influence the distributions for the range of depression densities we tested. We noted that for hyperbolic growth with mergers, increasing the initial number of depressions (n) slightly decreased the value of the power-law tail exponent. Exponential growth with mergers produced power-law size distributions when n ≥ 750, but power-law scaling only spanned ∼10 3 -10 5 m 2 .

The Scale Invariance of Depressions
Topographic depressions on the Ngozumpa Glacier were scale invariant. The size distribution exhibited robust power-law scaling (Figure 2a) over more than four orders of magnitude in area (∼10 1 -10 5 m 2 ) and the depressions tail exponent, γ ≈ 1.68, was within the range of tail exponents measured for sinkholes with power-law size distributions, 1.51 ≤ γ ≤ 1.98 (Pardo-Igúzquiza et al., 2020;Yizhaq et al., 2017). Positive feedback growth is the most widely applicable generative mechanism for power-law distributions (Newman, 2005), and power-law distributions cannot arise in the presence of strong negative feedbacks, which tend to produce characteristic scales (Evans, 2003). If differential melt drives depression growth, then we should expect the progressive deposition of debris within the depressions to act as a negative feedback (Figure 4a). The power-law distribution suggests negative feedbacks exert minimal influence on depression development and challenges the differential melt hypothesis for depression genesis.
Scale invariant size distributions are associated with fractal geometry (Mandelbrot, 1982). The depressions on the Ngozumpa Glacier exhibited uniform fractal scaling across the full range of depression areas (Figure 2b). Scaling breaks suggest a significant transition in growth processes. For example, the abrupt increase in the fractal dimension of Arctic melt ponds corresponded with a critical pond size where pond mergers increased (Hohenegger et al., 2012). We might expect a break in scaling to appear if variable debris thickness produced nascent depressions, but some other process (e.g., ice cliff retreat) took over as the depressions grew larger. If such a transition existed, these different processes would influence the perimeter-area scaling. The uniform fractal scaling of depressions on the Ngozumpa Glacier suggests that there is no distinct transition in the mode(s) of depression growth over a broad range of depression sizes.
The power-law depth-area scaling shows that depressions consistently deepen as they grow (Figure 2c), which is difficult to explain using differential melt, whereby debris would ultimately accumulate in depression bottoms. Instead, the observed depth-area scaling suggests the importance of other processes that deepen depressions, such as debris removal from the depressions via englacial meltwater drainage (Benn & Evans, 2014), sub-aqueous melt/calving in ponds , meltwater drainage incision (Anderson et al., 2021), turbulent atmospheric heat fluxes (Bonekamp et al., 2020), and/or collapse of englacial voids exposing new ice cliffs (Sakai et al., 2000;E. S. Miles et al., 2016;Yizhaq et al., 2017).

Positive Feedback Depression Growth
Only hyperbolic positive feedback growth (dA/dt ∝ A r , r > 1) produced size distributions of depressions similar to the Ngozumpa Glacier. Although mergers between depressions occur on the Ngozumpa Glacier, our model shows that mergers are not required to generate power-law distributions. Only strong positive feedback depression growth was needed.
The solution to Equation 2 for hyperbolic growth rates approaches a singularity. Although an infinitely large depression is unrealistic, Imja Tsho, a large proglacial lake on the nearby Imja-Lhotse-Shar Glaciers, provides an example of how positive feedback depression growth could progress on the Ngozumpa Glacier. From 1957 to 1976, a collection of supraglacial ponds and depressions rapidly grew to form Imja Tsho (Watanabe et al., 2009). By 1978, Imja Tsho spanned the width of the glacier, and continued to expand up the ablation zone (Watanabe et al., 2009). On the Ngozumpa Glacier from 2001 to 2010, the area of Spillway Lake more than doubled, growing to 2.7 × 10 5 m 2 Thompson et al., 2012). However, by the end of 2019, the area of Spillway Lake decreased to 2.3 × 10 5 m 2 . During this same period, approximately 1 km up-glacier of Spillway Lake, several large depressions merged to form a depression exceeding 2.5 × 10 5 m 2 . This suggests that mergers could be an important growth mechanism for large depressions and that Spillway Lake may undergo further expansion.

The Development of Hummocky Topography
Hummocky topography on debris-covered glaciers has been compared to karst topography (Clayton, 1964). Our analyses deepened this analogy by showing that the size distribution of depressions on the Ngozumpa Glacier is quantitatively similar to some types of karst sinkholes (Pardo-Igúzquiza et al., 2020;Yizhaq et al., 2017). We infer that depression genesis on debris-covered glaciers is analogous to karst sinkhole formation (Figure 4b). Sinkholes form via the focusing of bedrock dissolution at locations where surface water drains through cracks, often mobilizing soil through suffosion. Sinkholes also form via cave collapse. Similar to karst bedrock, only preexisting crevasses and permeable seams of englacial debris can route meltwater through glacier ice (Gulley, Benn, Screaton, & Martin, 2009). Human-traversable conduits exist within debris-covered glaciers and develop when meltwater drains through preexisting lines of permeability Gulley & Benn, 2007;Gulley, Benn, Screaton, & Martin, 2009). Surface observations (Benn et al., 2001;Benn et al., 2012;E. S. Miles et al., 2017;K. E. Miles et al., 2020), borehole investigations (K. E. Miles et al., 2022), and englacial cave surveys documented Gulley & Benn, 2007;Gulley, Benn, Müller, & Luckman, 2009, Gulley, Benn, Screaton, & Martin, 2009) small scale (∼10 −3 -10 −1 m thick) permeable structures within debris-covered glaciers. We infer that these structures are widespread and act as sink points for supraglacial meltwater and debris. Supraglacial meltwater drainage into these sink points produces nascent depressions by focusing incision and evacuating sediment at the drainage point.
The power-law size distribution of depressions on the Ngozumpa Glacier suggests that depressions undergo positive feedback growth along the observed geometric scaling relationships. Our depression growth model demonstrated that hyperbolic positive feedback growth produced similar distributions. In karst, sinkhole growth is a positive feedback process driven by the increase in drainage area (Pardo-Igúzquiza et al., 2020;Yizhaq et al., 2017). Similar dynamics may apply for depressions on debris-covered glaciers. We infer that the increase in drainage area increases meltwater available for channel incision, sediment transport, pond development, ice cliff production, and the growth of topographic relief. Englacial drainage is an essential aspect of positive feedback growth, because it limits negative feedback through enabling pond drainage and debris removal from depression bottoms.
Supraglacial ponds profoundly influence the development of depressions and glacier melt rates (K. E. Miles et al., 2020). Ponds grow the depressions through sub-aqueous melt and calving (E. S. Miles et al., 2016;Mertes et al., 2017), destabilizing debris slopes (Moore, 2018), and ice cliff retreat along the pond margins (Steiner et al., 2019). Importantly, hydrologically connected ponds create more ice cliffs and undergo faster deepening than isolated ponds (Watson et al., 2018). Thus, depressions with hydrologically connected ponds should grow faster than isolated depressions. During our survey, ponds were present in only ∼4.3% of the depressions (Text S3 in Supporting Information S1). Seasonal changes in englacial drainage pathways (Gulley, Benn, Screaton, & Martin, 2009) contribute to the cyclic filling and draining of ponds (E. S. Miles et al., 2017;Watson et al., 2018), and it is estimated that most glacier melt caused by ponds occurs during englacial drainage (E. S. Miles et al., 2016). Englacial drainage helps explain the abundance of depressions without ponds and suggests the possibility of many more transient ponds than can be detected during individual surveys.
If most depressions were not englacially drained, we might expect ponding to be uncorrelated with elevation. However, the median water level elevation in ponded depressions on the Ngozumpa Glacier was ∼20 m below the median minimum elevation in dry depressions ( Figure S2b in Supporting Information S1). Sustained depression growth eventually lowers the depressions to the elevation at which water exits the glacier, known as base level. When englacially drained depressions reach base level, depressions back-fill with water, forming permanent ponds. Our results illustrate how the widespread development of englacially-drained depressions precondition the ablation zone for the proliferation of supraglacial ponds and lakes. During this final stage of deglaciation, ponds and lakes permanently occupy depressions (Thompson et al., 2012;Watanabe et al., 2009).
Although the development of hummocky topography has been attributed to variable debris thickness (Benn & Evans, 2014;Sharp, 1949), landscape evolution models (Fowler & Mayer, 2017;Moore, 2021) only produce moderate, decaying relief from variable thickness debris. We propose an alternative hypothesis: hummocky topography develops from the sustained growth of englacially drained topographic depressions; the surrounding topography becomes elevated relative to the depressions, forming the hummocks ( Figure 4c). As with karst landscapes, isolated drainage basins characterize the hummocky topography on debris-covered glaciers . In the absence of englacial sink points, more ponds and continuous supraglacial drainage networks would develop. Our findings emphasize the role of englacial hydrology, rather than debris thickness, in the development of hummocky topography.

Conclusion
The abundance of topographic depressions is a defining characteristic of hummocky topography on debris-covered glaciers. We examined the size distribution and geometry of depressions on the Ngozumpa Glacier to better understand depression growth. The power-law size distribution, fractal perimeters, and power-law depth-area scaling showed that topographic depressions are scale invariant, and follow power-law perimeter-area and deptharea scaling relationships. Power-law size distributions do not arise in the presence of strong negative feedbacks. With our geometric model, we showed that strong positive feedback depression growth produces power-law size distributions similar to the depressions on the Ngozumpa Glacier. Hillslope debris transport and differential melt combine to produce a negative feedback to depression growth, thus, variable thickness debris alone cannot explain depression development. However, englacial drainage allows for pond drainage and debris removal, increasing depression longevity by reducing the influence of negative feedbacks. Crucially, englacial drainage explains why more than 95% of depressions did not contain ponds. These findings suggest that englacial drainage is a key contributor to the development of hummocky topography. This raises the possibility that predicting melt rates on debris-covered glaciers could hinge on understanding how englacial drainage networks evolve.