The Role of Differential Ablation and Dynamic Detachment in Driving Accelerating Mass Loss From a Debris‐Covered Himalayan Glacier

Sustained mass loss from Himalayan glaciers is causing supraglacial debris to expand and thicken, with the expectation that thicker debris will suppress ablation and extend glacier longevity. However, debris‐covered glaciers are losing mass at similar rates to clean‐ice glaciers in High Mountain Asia. This rapid mass loss is attributed to the combined effects of; (a) low or reversed mass balance gradients across debris‐covered glacier tongues, (b) differential ablation processes that locally enhance ablation within the debris‐covered section of the glacier, for example, at ice cliffs and supraglacial ponds, and (c) a decrease in ice flux from the accumulation area in response to climatic warming. Adding meter‐scale spatial variations in supraglacial debris thickness to an ice‐flow model of Khumbu Glacier, Nepal, increased mass loss by 47% relative to simulations assuming a continuous debris layer over a 31‐year period (1984–2015 CE) but overestimated the reduction in ice flux. Therefore, we investigated if simulating the effects of dynamic detachment of the upper active glacier from the debris‐covered tongue would give a better representation of glacier behavior, as suggested by observations of change in glacier dynamics and structure indicating that this process occurred during the last 100 years. Observed glacier change was reproduced more reliably in simulations of the active, rather than entire, glacier extent, indicating that Khumbu Glacier has passed a dynamic tipping point by dynamically detaching from the heavily debris‐covered tongue that contains 20% of the former ice volume.

. Low climatic sensitivity promoted the development of long, low-angle debris-covered tongues that have a greater longevity than climatically equivalent clean-ice glaciers (Rowan et al., 2015). However, debris-covered glaciers that have historically persisted at lower elevations than climatically equivalent clean-ice glaciers due to their reduced climatic sensitivity now show a differing response to recent climate forcing; as air temperatures rise, debris-covered glacier termini stabilize while rates of mass loss accelerate further upglacier (Nuimura et al., 2011). Regional mass loss from debris-covered glaciers in High Mountain Asia is occurring at similar rates to mass loss from clean-ice glaciers (e.g. Brun et al., 2019;Gardelle et al., 2013) an observation called the "debris-cover anomaly" (Pellicciotti et al., 2015). A paradox exists in the hypothesis that such high rates of mass loss from debris-covered glaciers is the result of enhanced climatic forcing at the glacier surface, as greater surface melting across the ablation area of a debris-covered glacier will increase the rate of englacial debris emergence, resulting in thickening of supraglacial debris and reduced ablation over decadal timescales. The following mechanisms have been suggested as possible causes of the debris-cover anomaly (Pellicciotti et al., 2015): 1. Increasing climatic sensitivity of low-angle debris-covered tongues due to suppressed, or reversed, ablation gradients as the climate warms (Benn et al., 2012). 2. The development of ice cliffs and supraglacial ponds that disrupt the insulating debris layer and promote differential ablation Buri et al., 2021). 3. Decreasing mass flux from glacier accumulation areas to debris-covered tongues (Banerjee, 2017) which is most pronounced when the active glacier detaches from the debris-covered tongue (Rowan et al., 2015).
Differential ablation occurs because debris thickness varies laterally on a meter-scale due to undulating glacier surface topography, from close to zero at ice cliffs to several meters thick where debris accumulates in topographic depressions (Benn et al., 2012;Mölg et al., 2020). The term differential ablation is used herein to refer to the combined effect of debris thickness variations, the presence of ice cliffs and supraglacial ponds, and other processes occurring within supraglacial debris that affect mass balance (e.g., seasonal supraglacial streams). Given this spatial heterogeneity, because insulation is a non-linear function of debris thickness, a positively skewed distribution of debris thickness will generally increase net ablation relative to that calculated using a mean debris thickness value . Ignoring these local variations in debris thickness can lead to underestimates of sub-debris ablation rates of 11%-30% (Nicholson et al., 2018). Ice cliffs form ablation "hotspots" on the surface of debris-covered glaciers (Watson et al., 2017) and although ice cliffs occupy a relatively small part of the debris-covered area, they contribute disproportionately to total ablation as the dominant mechanism of mass loss where debris is thick . On Ngozumpa Glacier, ice cliffs occupy 5% of the area and account for 40% of ablation from the debris-covered tongue (Thompson et al., 2016). On Changri Nup Glacier, ice cliffs occupy 8% of the debris-covered tongue and account for 24% of total ablation, with a net ablation rate three times higher than the glacier-wide mean . Most ice cliffs are associated with supraglacial ponds that also enhance ablation (Thompson et al., 2016;Watson et al., 2017). In the Langtang valley, Nepal, ponds cover 0.8% of the debris-covered glacier area and enhance melt by a factor of 14, accounting for 13% of total ablation (Miles, Willis, et al., 2018). The area occupied by ice cliffs and supraglacial ponds has expanded on Khumbu Glacier by at least 20% across the debris-covered area between 1984 and 2016, suggesting that the impact of these features on ablation is increasing (King, Turner, et al., 2020).
The rate of mass loss from glaciers in the Himalaya has accelerated over the last 40 years, with the rate of ice mass loss doubling to −0.43 ± 0.14 m water equivalent (w.e.) a −1 during 2000-2016 compared to 1975-2000 in response to climate change (Maurer et al., 2019). The acceleration in mass loss is expressed by the rate of surface elevation change, which has increased dramatically at Khumbu Glacier since 1969 (King, Bhattacharya, et al., 2020). Glacier mass loss is occurring due to atmospheric warming and is regionally variable, due to topographic feedbacks between climate systems and ice flow (Dehecq et al., 2018). Debris-covered glaciers in the Everest region have slowed during recent decades, resulting in a typical pattern of fast ice flow (>50 m a −1 ) through steep icefalls and slow-flowing or stagnant ice (<10 m a −1 ) in their debris-covered tongues (Quincey et al., 2009). Spatial patterns in velocity therefore differ from those expected for a debris-covered glacier closer to equilibrium with climate or for a clean-ice glacier indicating an ongoing and pronounced reduction in mass flux from the accumulation area to the debris-covered ablation area termed "dynamic detachment" (Quincey et al., 2009). Dynamic detachment occurs as a precursor to physical detachment of the upper active section of a glacier from the debris-covered tongue, after which point, no further mass is added to the tongue by ice flow and the former accumulation area becomes the entire (flowing) glacier (cf. Rippin et al., 2020).
Numerical models can simulate the dynamic response of glaciers to climate change (e.g. Egholm et al., 2012;Jouvet et al., 2009); however, modeling debris-covered glaciers is challenging because important feedbacks between debris transport, ice flow and mass balance are often not represented, leading to contradictory predictions of mass change, that, in the Everest region, range from loss of 10% to more than 80% of glacier ice over the same period (e.g. Rowan et al., 2015;Shea et al., 2015;Soncini et al., 2016). A small number of recent modeling studies have included these feedbacks to explore the processes affecting the behavior of debris-covered glaciers (Anderson & Anderson, 2018;Wirbel et al., 2017). Here, we use an ice flow model that includes the feedbacks between debris transport, ice flow and mass balance with a novel parameterization of differential ablation to test the hypothesis that the mechanisms suggested to cause the debris-cover anomaly caused the acceleration in mass loss from Khumbu Glacier observed over a 31-year period from 1984 to 2015.

Glacier Morphometry
Khumbu Glacier is 16 km long, 19 km 2 in area and has a 12 km ablation area extending from the base of the icefall at 5,400 m above sea level (a.s.l.) to the Little Ice Age (LIA) terminal moraine at 4,900 m a.s.l.
(RGI Consortium, 2017) ( Figure 1a). The LIA moraines formed during the last advance 500 years before the present day when mean annual air temperatures were about 1.5°C colder (Rowan, 2017). Ice thickness measured along seven transects using radio-echo sounding was about 440 m in the upper ablation area decreasing to 110 m in the lower ablation area and less than 20 m close to the terminus (Gades et al., 2000;Moribayashi, 1978). Ice thickness from a global consensus estimate gives a maximum of 341 m and a mean of 88 ± 20 m (Farinotti et al., 2019). The glacier has thinned and stagnated as the debris layer thickened and the debris surface has remained relatively stable over the last 40 years (Iwata et al., 2000;Nakawo et al., 1986). The ablation area lost up to 65 m of ice between 1984 and 2016 ( Figure 1b) (King, Turner, et al., 2020). The area between the LIA terminal moraine and the debris-covered terminus (about 1.5 km distant) is occupied by thin, decaying ice covered with debris several meters in thickness.
The geodetic mass balance was −0.26 ± 0.13 m w.e. a −1 between 1984 and 2001, and -0.44 ± 0.12 m w.e. a −1 between 2001 and 2018 (King, Bhattacharya, et al., 2020) in agreement with geodetic values from other studies of −0.27 ± 0.08 m w.e. a −1 between 1970 and 2007 (Bolch et al., 2011) and −0.45 ± 0.52 m w.e. a −1 between 1992 and 2008 (Nuimura et al., 2012). Ice cliffs are found in close proximity to areas occupied by thick debris as the geometry of supraglacial hummocks promotes both the formation of ice cliffs where hummocks collapse and the accumulation of thick debris in the troughs between hummocks, with a mean hummock spacing of 162 m across glacier and 196 m along glacier (Bartlett et al., 2021). Ice cliffs and supraglacial ponds occupy 5% of the debris-covered area, reaching a maximum density 5 km upglacier from the terminus and becoming less extensive in the upper ablation area where surface lowering reaches a maximum (Watson et al., 2017).

Supraglacial Debris Thickness
Supraglacial debris thickness varies dramatically over short length scales on the surface of Khumbu Glacier, and significant surface roughness reflecting differential ablation occurs on a decimeter scale (Bartlett et al., 2021). Glacier-wide debris thickness, excluding areas occupied by ice cliffs, estimated remotely using inversions of satellite measurements of surface temperature and elevation change gave a median debris thickness of 1.9 m and a 75% centile of 2.7 m in the lowermost 1.5 km of Khumbu Glacier (Rounce et al., 2018). Including supraglacial ponds and ice cliffs in this estimate reduced the median debris thickness to 1.1 m with debris thickness decreasing to less than 0.2 m above 3.4 km upglacier of the terminal moraine (Rounce et al., 2018). We measured debris thickness across the lower ablation area in 2014 and 2015 in areas with continuous debris cover, excluding ice cliffs, supraglacial ponds and large boulders often several meters in diameter. From the terminus to 3.5 km upglacier, 80% of points measured (n = 143) had debris greater than 1 m thick, with thinner debris (0.04-1.0 m) found around supraglacial ponds or overlying ice cliffs (Gibson et al., 2018). Other studies have measured debris thicknesses of 0.5-2.0 m overlying ice cliffs on Khumbu Glacier (n = 50) (Nakawo et al., 1986) and a mean debris thickness of 0.35 m with a maximum of 3.0 m (n = 64) across a similar area (Soncini et al., 2016). We note that there are limitations in estimating glacier-wide supraglacial debris thickness. Remote sensing methods typically underestimate debris thickness, as they rely on the assumption that the thermal properties of supraglacial debris are stable over representative timescales for calculating thermal flux through the debris layer or ablation beneath debris (e.g. Rounce et al., 2018). Field measurements typically overestimate debris thickness, as they are usually made where debris is thickest in the lower ablation area of the glacier and typically do not sample the unstable sections of the glacier surface where thin debris occurs (e.g. Gibson et al., 2018). Both methods are unlikely to capture the thickest debris, which is poorly resolved in remote sensing methods and challenging to measure directly. However, the two methods can be combined to estimate glacier-wide debris thickness as they indicate the shape of the distribution Nicholson et al., 2018;Stewart et al., 2021) and the likely range of values (Rounce et al., 2018).

Glacier Dynamics
Velocities calculated from surface displacements obtained by feature tracking indicate that ice flow through Khumbu Glacier declines rapidly below the upper ablation area where the glacier surface becomes debris covered (Quincey et al., 2009;Rowan et al., 2015). Velocities in the icefall (10.8-11.5 km upglacier from the terminal moraine) reached a maximum of about 100 m a −1 between 2016 and 2017 (Watson & King, 2018) and over one meter per day in May 2018 (Altena & Kääb, 2020). Velocities for the ablation area generated using auto-RIFT  provided by the NASA MEaSUREs ITS_LIVE project indicate ice flow of 30-50 m a −1 in the upper ablation area (6.0-10.8 km from the terminal moraine) declining rapidly to less than 10 m a −1 the lower ablation area (Gardner et al., 2019). ITS_LIVE data from above the base of the icefall (10.8 km upglacier from the terminal moraine) were not used as they suffer from "surface skipping or locking" such that robust displacements are not recorded (Dehecq et al., 2018). Mean emergence velocities across the ablation area estimated from ITS_LIVE data are 0.28 ± 0.32 m w.e. a −1 , although a measurement gave a higher value of 5.0 m w.e. a −1 (Nuimura et al., 2011). The pronounced decrease in velocity between the upper and lower ablation areas around 6.0 km upglacier of the terminal moraine correlates with the downglacier margin of the region where surface lowering is greatest between the Changri Nup palaeoconfluence and the base of the icefall (Figure 1b). Observations of skewed internal debris layers, indicating englacial thrusting, and the presence of basal ice at the glacier surface in a similar location  suggest that the upper active glacier may have dynamically detached from the heavily debris-covered tongue.

Model Description and Experimental Design
Experiments were performed using the ice-flow model iSOSIA (Egholm et al., 2012) that couples the transport of debris with mass balance as presented in Rowan et al. (2015). Debris accumulates at the surface of Khumbu Glacier from high-energy avalanches that can spread rock and snow across large sections of the accumulation area (Takeuchi et al., 2020). To avoid assumptions about the distribution and frequency of avalanches in our model, debris was added to the glacier surface in the accumulation area at a constant and spatially uniform rate equivalent to headwall erosion of 1 mm a −1 . Debris was advected through the ice in a 3-D finite volume scheme until it emerged at the glacier surface in the ablation area. Our 3-D model accounts for across-glacier variations in ice thickness and stress, and therefore, as all debris is sourced from ROWAN ET AL.
10.1029/2020JF005761 5 of 20 the accumulation area, the distribution of debris in the ablation area is a function of mass conservation such that a thinner debris layer is found, for example, where ice flow lines are not parallel but diverge, as they often do in the ablation area. In a 1-D model (e.g. Anderson & Anderson, 2016 debris would consistently thicken downglacier, however, our 3-D model has more realistic across-and along-glacier controls on debris thickness. We allowed the evacuation of debris from the glacier by erosion of the outer slopes of the ice-marginal moraines when they reached a critical slope of 40°. The critical slope was quantified from observations of the maximum slope of the LIA lateral and terminal moraines from an 8-m Digital Elevation Model (DEM) obtained via the National Snow and Ice Data Center (Shean, 2017). Ice-marginal moraine erosion was implemented using the hillslope sediment transport functions described by Egholm et al. (2012).
Ice thickness for Khumbu, Changri Nup and Changri Shar Glaciers was estimated in Rowan et al. (2015) using a perfect-plasticity relationship defined by Nye (1952) assuming that ice thickness is determined by surface slope and basal shear stress, using a mean value for basal shear stress of 150 kPa and constrained by geophysical measurements of ice thickness (Gades et al., 2000;Moribayashi, 1978). Maximum estimated ice thickness was 391 m with a mean of 123 ± 29 m. Estimated ice thickness was subtracted from the 30-m ASTER DEM and resampled to 100-m grid spacing to give subglacial topography for the model domain. Surface mass balance was calculated for an initial clean-ice surface then the simulated debris thickness in each grid cell was used to estimate sub-debris melt. Mass balance was calculated using values from Benn and Lehmkuhl (2000) to give an accumulation gradient of 0.2 m w.e. m −1 and an ablation gradient of 0.005 m w.e. m −1 limited to a maximum of ±2.0 m w.e. a −1 . Snow was removed by avalanching from slopes steeper than 28° and redistributed following a path of steepest descent. Mean annual air temperature (MAAT) was specified to define the equilibrium line altitude (ELA), and change in MAAT (hereafter ΔT) used to impose climatic forcing of mass balance. The ELA was prescribed using a lapse rate of −0.004°C m −1 to give 5,625 m a.s.l. during the LIA and 6,000 m a.s.l. during the present day if the glacier was in mass-balance equilibrium with regional climate.
The ice-flow model was spun up to simulate the Late Holocene advance of Khumbu Glacier (∼1,000 CE) and recession to the LIA maximum using the same approach as Rowan et al. (2015). We perturbed the mass balance function to investigate the impact of differential ablation on glacier evolution. The impact of differential ablation depends on the development of supraglacial debris, which was minimal during the LIA when debris was efficiently exported from the glacier to form ice-marginal moraines. Debris began to accumulate across the glacier surface about 100 years after the onset of climate warming following the LIA maximum (Fushimi, 1977;Nakawo et al., 1986;Owen et al., 2009;Rowan et al., 2015). We made transient simulations of glacier evolution from the LIA maximum to the present day (2015 CE) using a range of parameterizations of sub-debris melt to represent the differences between a continuous and a discontinuous debris layer (i.e., one that is punctuated by ice cliffs and supraglacial ponds). In addition, we simulated the behavior of the upper active glacier separately from the heavily debris-covered tongue to explore the possibility of dynamic detachment.

Ablation Beneath a Continuous Debris Layer
The reduction in ablation from clean-ice values beneath continuous supraglacial debris can be represented as an exponential function (Rowan et al., 2015) or a reciprocal function (Anderson & Anderson, 2016). A reciprocal function represents a satisfactory fit between average cell-values of debris thickness (h) and sub-debris melt rate (b debris ): where h 0 is a constant that represents the characteristic debris thickness at which the reduction in ablation due to insulation by supraglacial debris is 50% of the value for an equivalent clean-ice surface (b clean ). The h 0 value for sub-debris melt was calculated by fitting Equation 1 to measurements of sub-debris ablation normalized to measurements of clean-ice ablation made using ablation stakes at Khumbu Glacier (Inoue, 1977;Inoue & Yoshida, 1980;Kayastha et al., 2000) and Changri Nup Glacier (Vincent et al., 2016). These measurements indicate a reduction in ablation up to 40% of clean-ice values beneath 0.5 m of debris (Inoue & Yoshida, 1980) and a mean value for sub-debris ablation in the lower ablation area where debris was around 0.8 m thick of 25% ablation from clean ice in the upper ablation area (Inoue, 1977) (Figure 2b).
We used 12 ablation stake measurements for debris thicknesses ranging from 0.06 to 2.15 m to calculate h 0 . Where sub-debris ablation was given as ice thickness, we converted these values to water equivalent using a density of ice of 900 kg m −3 . Where ablation was given as daily values, we converted to annual values assuming the length of the sub-debris ablation season at 5,300 m a.s.l. is 160 days . We normalized the measurements of sub-debris ablation to the clean-ice ablation rate at the same site and fitted Equation 1 to the data. The best fit of Equation 1 to these data was given where h 0 = 0.23 m (r 2 = 0.72; Figure 2b) similar to the h 0 values of 0.25 m calculated from inversion of surface elevation change (Rounce et al., 2018) and 0.17 m used in a regional model that included Khumbu Glacier (Rounce et al., 2020). The ablation stake measurements indicated a smaller reduction in sub-debris melt with debris thickness at Khumbu Glacier relative to debris-covered glaciers outside High Mountain Asia (e.g. Östrem, 1959;Anderson et al., 2021). We assumed no enhancement of ablation where debris was thin, that is, below the critical thickness where sub-debris melt is equal to clean-ice melt for a climatically equivalent surface (Anderson & Anderson, 2016) estimated as 0.02 m in this region (Rounce et al., 2018). Although some studies report an increase in ablation on debris-covered glaciers in the Everest region relative to that expected for a climatically equivalent clean-ice surface (e.g. Kayastha et al., 2000), there is only a small area where a sufficiently thin continuous debris cover is found on Khumbu Glacier close to the base of the icefall where no such increase in melt was observed, which was attributed to high daytime air temperatures and high humidity during the ROWAN ET AL.
10.1029/2020JF005761 7 of 20 to reduce sub-debris melt from clean-ice values, compared with the exponential function used in Rowan et al. (2015). The clean-ice melt reduction factors for sub-debris melt calculated from 12 ablation stake measurements for Khumbu and Changri Nup Glaciers are presented in (b) as dots (Kayastha et al., 2000), crosses (Inoue, 1977;Inoue & Yoshida, 1980) and a triangle (Vincent et al., 2016). Results for Isfallsglacier in arctic Sweden from Östrem (1959) are shown for comparison.
ablation season (Inoue & Yoshida, 1980). Thin discontinuous debris layers are also found on the surface of ice cliffs and their impact was incorporated into our calculations of differential ablation.

Differential Ablation Beneath a Discontinuous Debris Layer
We tested a novel parameterization of sub-debris melt that represents the impact of differential ablation on mass balance. We  (Table 1). This approach allowed us to calculate sub-debris melt for each distribution of h from the value in each cell without calculating the distribution of h within the cell. The distributions of h were chosen to give realistic values from debris thicknesses observed on glaciers in the Everest region. The mean and standard deviation of each distribution were chosen to ensure that the total debris thickness was the same and only the spatial distribution of debris changed. The skewness and kurtosis gave the same form for each distribution and were chosen based on field measurements of debris thickness around ice cliffs from Ngozumpa Glacier, located 10 km west of Khumbu Glacier ) and a similar study using ground-penetrating radar (Nicholson et al., 2018) that observed positively skewed distributions of h (Table 1). The simulated positively skewed distributions were symmetrical to the negatively skewed distribution, with the maximum value of h ranging from 1.0 to 4.0 m to account for uncertainty in defining the maximum debris thickness found on Khumbu Glacier. To account for the presence of very thin debris or clean-ice surfaces on ice cliffs where h is at or close to zero, we set the minimum h value to zero for all distributions. The positively skewed distributions gave the highest values for h 0 as they represent the ablation "hotspots" that are found where debris is thin around ice cliffs and supraglacial ponds. Our values for h 0 are higher than the value of 0.066 ± 0.029 m used for a continuous debris layer to represent the propagation of heat through a granular medium by Anderson and Anderson (2016) because, in addition to accounting for the reduction of ablation with increasing debris thickness beneath a continuous debris layer, we explicitly include the impact of differential ablation in our h 0 values as thin or debris-free points (representing ice cliffs) contribute to the effective mean ablation in each grid cell. For example, for a 100 m by 100 m cell with a maximum debris thickness of 2.0 m, observations of the glacier surface indicate that it will contain ice cliffs where melt is locally 3-14 times higher . The mean sub-debris melt for that cell calculated using h 0 = 0.94 m would therefore be double that calculated assuming a continuous debris layer but 68% lower than that for a cell with an entirely clean-ice surface (Figure 2b).
ROWAN ET AL.  (2017) and using ground-penetrating radar (GPR) by Nicholson et al. (2018). The form of the distributions is illustrated in Figure 2a.

Dynamic Detachment of the Active Glacier From the Debris-Covered Tongue
We tested the hypothesis suggested by Rowan et al. (2015) that Khumbu Glacier can maintain higher velocities in the upper ablation area compared to the lower ablation area because the upper active glacier is dynamically detached from the debris-covered tongue close to the Changri Nup palaeoconfluence. The timing of the onset of dynamic detachment is unknown but can be estimated from observations of the surface expression of the englacial structure. In a structural map of Khumbu Glacier from the 1970s, the active terminus was identified approximately 1 km downglacier of the Changri Nup palaeoconfluence with a 1-km long section above this point identified as containing a high concentration of debris and evidence of intense shear bringing material from the glacier bed to the surface (Fushimi, 1977). From Fushimi's observations we can infer that dynamic detachment occurred during or prior to the 1970s. Fushimi (1977) notes that shear occurs between the lowermost 30% of the ice thickness and the overlying ice along the active part of Khumbu Glacier. We simulated dynamic detachment of Khumbu Glacier from the debris-covered tongue by subtracting the simulated active ice volume from the simulated LIA ice volume and incorporating the remaining ice volume into the subglacial topography of the model domain. The new model domain therefore contained a stagnant ice volume that occupied the LIA lateral and terminal moraines downglacier of the Changri Nup palaeoconfluence. The glacier bed above this point dipped upglacier and had a higher elevation than the bed topography used in the previous simulations, representing the ice thickness underlying the shear plane that separates the active glacier from the underlying basal ice (Fushimi, 1977). The simulation was spun up to represent the active glacier with detachment occurring after 1900 CE, forced with the same values for ΔT and h 0 as previous simulations of the entire glacier to reach the present day 200 years after the LIA maximum.

Model Evaluation
Results were evaluated against measurements of ice thickness, supraglacial debris thickness, net mass balance, and horizontal (surface) and vertical (emergence) velocities described in Section 2.1 to evaluate how well each simulation represented the present-day glacier. Simulations of the active glacier were evaluated against data with the same extent. As velocities calculated by the glacier model are depth-integrated, the simulated deformation rate was multiplied by 1.2 then added to the sliding rate to give surface velocities equivalent to those measured using feature tracking. Ice flux was calculated from surface velocity and simulated ice thickness for each cross-section of the 500-m wide swath profile (Figure 1b, inset) and used to calculate emergence velocity by dividing ice flux at each cross-section by the downglacier area. To evaluate how well the model reproduced the recent acceleration in mass loss, each simulation was evaluated by minimizing the root mean squared error (RMSE) and bias between simulated with observed surface elevation change ( The 2015 DEM of Khumbu Glacier (hereafter WV-DEM) was generated from several pairs of WorldView-1 and WorldView-2 images using the Surface Extraction with TIN-based Search-Space Minimization (SETSM) algorithm described by Noh and Howat (2015). The WV-DEM covering the Khumbu valley is a composite of two DEMs, one derived from imagery acquired on January 31, 2015 and one acquired from imagery acquired on February 2, 2015. We mosaicked the two DEMs assuming no significant surface elevation change occurred during this two-day gap. The 1984 DEM of Khumbu Glacier (hereafter AP-DEM) was derived from fine resolution (0.5 m ground resolution) aerial photographs acquired in 1984 using a Wild RC-10 camera (Washburn, 1989). We used seven images from the Bradford Washburn aerial photo set, purchased from Swissphoto and scanned at 1693 dpi from the original diapositives. DEM processing was carried out using PCI Geomatica software (2018 version). The geolocation of our AP-DEM was fixed using 14 Ground Control Points collected in the Khumbu valley in October 2015 which had a mean 3-D positional uncertainty of 5.1 mm.
We removed any remnant geolocation error between the location of the WV-DEM, SRTM DEM and the AP-DEM using the methods of Nuth and Kääb (2011). The AP-DEM was used as the master DEM during the co-registration process because of the higher accuracy of GCPs used in its generation. After co-registration, we differenced each pair of DEMs and used a similar approach to the two-step filtering process described by Gardelle et al. (2013) to remove outliers resulting from errors produced during the DEM extraction process. We applied an initial threshold of >±100 m to remove obvious errors in DEM difference data and then a secondary threshold of ±3 times the standard deviation of remaining difference data where surface slope was less than 30°. The mean and standard deviation of off-glacier AP-DEM-WV-DEM elevation differences were 1.6 and 7.9 m following co-registration and outlier removal. The mean and standard deviation of off-glacier SRTM-DEM-WV-DEM elevation differences were −0.90 and 2.25 m following co-registration and outlier removal.

Simulations Assuming a Continuous Debris Layer
Simulations assuming a continuous debris layer using a h 0 value of 0.23 m were forced by ΔT of 1.5°C, 2.0°C, 2.5°C or 3.5°C between the LIA maximum and the present day (Table 2). In each simulation, ice thickness decreased along the entire length of the glacier (Figure 3a). The simulation where ΔT = 2.5°C produced greater ice loss from the ablation area compared to those simulations that used either a higher or lower value of ΔT. Simulated present-day ice thickness was greater than observed values as a thick debris layer formed across the entire ablation area (Figure 3c). The greatest amount of mass loss occurred at the base of the icefall above the debris-covered section with relatively little mass loss below this point compared to the ROWAN ET AL.   For the same simulations, mean simulated debris thickness in 2015 CE for (c) simulations assuming a continuous debris layer, (d) simulations assuming a discontinuous debris layer. Note that the decrease in debris thicknesses between the terminal moraine and up to 2 km upglacier occurs because supraglacial debris is exported when it reaches the ice margin and no longer contributes to the results, therefore the supraglacial debris thickness appears less extensive than that presented in Rowan et al. (2015). Simulated and observed velocities from the NASA MEaSUREs ITS_LIVE project (Gardner et al., 2019)  icefall and lower ablation area but lower than ITS_LIVE values through the upper ablation area (Figure 3e). The mean simulated horizontal velocity in the upper ablation area (6.0-10.8 km upglacier from the terminal moraine) was 10 m a −1 and the mean simulated emergence velocity through the ablation area was 0.11 m w.e. a −1 .

Simulations Assuming a Discontinuous Debris Layer
Simulations assuming a discontinuous debris layer accounted for the impact of differential ablation where the distribution of debris thickness (h) within each cell had a uniform, skewed or normal distribution. For example, sub-debris melt calculated using a positively skewed distribution of h with a maximum value of 2.0 m was 29% greater where debris was 0.2 m thick and 34% greater where debris was 0.5 m thick compared to a continuous debris layer of the same thickness (Figure 2b). Debris thickness increased with increasing h 0 , from a mean of 0.56 ± 0.11 m where h 0 = 0.13 m to a mean of 0.80 ± 0.25 m where h 0 = 1.55 m (Figure 3d). A value for ΔT higher than 1.5°C did not improve the fit between simulated and observed δz, instead most ice was lost at the base of the icefall rather than further downglacier resulting in physical detachment of the active glacier before 2015 CE. The simulations representing uniform, normal or negatively skewed distributions of h (h 0 < 0.5 m) reproduced no more than 43% of observed δz ( Table 2). The simulations with positively skewed distributions of h (h 0 > 0.5 m) gave better results, reproducing 60%-80% of observed δz between 1984-2015 CE and 50%-60% between 2000-2015 CE (Figures 3h and 3j).
The RMSE and bias decreased and proportion of observed δz simulated increased with increasing h 0 (Table 2) and all h 0 values representing a positively skewed distribution of debris thickness gave a reasonable match to observations (Figures 3b and 3d). We consider that the simulation using a value of 0.94 m for h 0 that assumes a maximum debris thickness of 2.0 m gave the most convincing result without overestimating δz (Figures 3h and 3j). For this simulation, the present-day ELA was 6,431 m a.s.l., the mean ice thickness was 110 ± 26 m, the mean debris thickness was 0.76 ± 0.23 m, and the net mass balance for a ten-year period up to the present day was −0.52 ± 0.40 m w.e. a −1 . Simulated velocities were similar to observed values in the icefall and lower ablation area but lower than ITS_LIVE values through the upper ablation area (Figure 3f). The mean simulated surface velocity in the upper ablation area was 13 m a −1 and the mean simulated emergence velocity through the ablation area was 0.13 m w.e. a −1 .

Simulations Assuming Dynamic Detachment of the Active Glacier From the Debris-Covered Tongue
The active glacier terminus remained in contact with the stagnant debris-covered tongue during all simulations, initially expanding slightly for 100 years as the model was spun up to the LIA maximum before the onset of sustained mass loss through the present day. After a further 100 years, the simulation of the active glacier lost mass reproduced observations of the present-day glacier (Figure 4) suggesting that dynamic detachment has occurred within the last 100 years and that the glacier has changed more rapidly since the end of the LIA (∼1900 CE) than can be represented using  (Table 2; Figures 4d and 4e).
In the simulation of the active glacier with a value for h 0 of 0.94 m, the present-day ELA was 6,239 m a.s.l., the mean ice thickness was 68 ± 29 m, the mean debris thickness was 0.51 ± 0.26 m and the net mass balance for a ten-year period up to the present day was −0.13 ± 0.11 m w.e. a −1 . The mass balance of the active glacier was less negative than that for the entire glacier because the debris-covered tongue is excluded from this simulation (Figure 5g). Simulated surface velocities were similar to observed values in the icefall and lower ablation area, and more similar to ITS_LIVE values through the upper ablation area than for the simulations of the entire glacier (Figure 4c). The mean simulated surface velocity in the upper ablation area was 15.6 m a −1 with a maximum of 72 m a −1 , and the mean simulated emergence velocity through the ablation area was 0.25 m w.e. a −1 ; similar to the value of 0.28 m w.e. a −1 estimated from the ITS_LIVE data. The present-day ice volume for the active glacier simulation is 80% of the volume of the present-day ice volume  for the entire glacier with the same parameterization. Allowing the active glacier simulation to reach equilibrium with present-day climate indicated a committed volume loss with no further change in climate of 77% of the active glacier volume (81% of the entire glacier volume) in the next 200 years (Figure 4i).

Discussion
In this section, we discuss the role of each of the three factors suggested to contribute to the debris-cover anomaly (a low or reversed mass balance gradient, differential ablation, and reduced mass flux leading to the extreme case of dynamic detachment) in driving mass loss from Khumbu Glacier between 1984 and 2015.

Uncertainties Associated With Glacier Dynamics
The simulations presented for the active and entire glacier assuming differential ablation beneath a discontinuous debris layer, that is, where debris thickness varies on a meter-scale and each model cell includes ablation "hotspots" such as ice cliffs and supraglacial ponds, resulted in values that better reproduced the present-day ice thickness and extent, debris thickness and net mass balance observed for the present-day glacier, and observations of recent glacier surface change (δz), than those simulations assuming a continuous debris layer without ablation "hotspots" when a positively skewed distribution of supraglacial debris including ice cliffs and supraglacial ponds was used (h 0 > 0.5 m) (Figures 3 and 4). Ice flow was similar to observations in the lower ablation area in all simulations but underestimated observations in the upper ablation area in simulations of the entire glacier. Simulations of the active glacier extent resulted in faster ice flow through the upper ablation area that was more similar to the observed pattern of surface velocity ( Figure 4c).
However, although our higher-order ice-flow model captures the important stresses controlling the flow of ice though rugged mountainous terrain, it was not possible to precisely represent the transient dynamics of a debris-covered glacier in an advanced state of decay. Reducing the temperature-dependent flow constant (A) in Glen's law in simulations of the entire glacier to invoke the relatively warm ice found within Khumbu Glacier (Miles, Hubbard, et al., 2018) increased velocity through the icefall and above the Changri Nup palaeoconfluence but did not reproduce the high rates of ice flow through the section between these points where the ice thickness was greatest. Therefore, the high velocities observed in feature-tracking measurements (e.g. ITS_LIVE; Quincey et al., 2009) cannot be attributed to the plastic flow of ice through the entire glacier. We chose not to impose extreme parameterizations of basal sliding and hydrology that could potentially reproduce the observed patterns of ice flow in the upper ablation area for two reasons; (a) because data describing these properties with which to constrain our parameterization of basal conditions are not available, and (b) instead of setting up an unusual representation of glacier dynamics to reproduce short-term variability we chose to work with a system that we know is representative of longer-term glacier behavior.
There is much to be learned about the behavior of debris-covered glaciers in general from a more theoretical approach (e.g. Anderson & Anderson, 2018) and some of these limitations could be addressed using a model containing the full Stokes equations that govern ice flow (e.g. Wirbel et al., 2017).
The LIA advance of Khumbu Glacier occurred about 500 years before the present day when the glacier surface was relatively clean as debris was efficiently exported to the glacier margins to build substantial ice-marginal moraines (Owen et al., 2009;Rowan, 2017). When climate warmed after the LIA, net glacier mass balance became negative but the resulting increase in ablation led to increased exhumation of englacial debris and expansion of supraglacial debris that forced a short-lived expansion of Khumbu Glacier (Rowan et al., 2015). After the LIA maximum, the glacier entered a period of sustained mass loss as the climate continued to warm to the present day. The end of the LIA, indicated by smaller ice-marginal moraines inset within those formed during the LIA maximum, is suggested to have occurred around 1900 CE (Rowan, 2017). We therefore considered two different timescales to simulate the behavior of the entire glacier and the active glacier. For the entire glacier, these simulations represent the response of Khumbu Glacier to climate change from the LIA maximum when moraine building finished and the supraglacial debris layer established, which required 400 years of glacier evolution. For the active glacier, these simulations are assumed to start later in the glacier's history, after the LIA maximum and when the debris layer has expanded across the lower ablation area sufficiently to reduce ice flow through the tongue. The evolution of the glacier from this point required 200 years to simulate, an initial 100 years to allow the active terminus to stabilize followed by 100 years of mass loss, suggesting that dynamic detachment of the active glacier from the debris-covered tongue occurred within the last century. Therefore, ΔT of 0.75°C over 200 years was required to force the active glacier from the onset of detachment to the present day, which is the same rate of change as ΔT of 1.5°C over 400 years required to force the entire glacier from the LIA maximum to the present day. If detachment of the active glacier did occur after 1900 CE, then we would expect mass loss from the debris-covered tongue to exceed that produced in the simulations of the entire glacier, as after detachment occurred there would be no further input of ice from the accumulation area. The debris-covered tongue would then be subject to mass loss as a function of environmental conditions with limited or no accumulation by direct snowfall. Ablation would result in thickening of the debris layer that may allow the tongue to persist as an inactive ice mass.

Low or Reversed Mass Balance Gradient
Thick supraglacial debris insulates the underlying ice surface from atmospheric warming. However, the low-angle or reversed mass balance gradient that occurs due to the presence of such debris layers make a glacier extremely sensitive to small changes in air temperature once the climate has warmed sufficiently to give positive MAAT at the terminus. Compared to a climatically equivalent clean-ice glacier, which would have a steeper mass balance gradient, the low mass balance gradient of a debris-covered glacier allows relatively small increases in MAAT to affect a large part of the ablation area (Pellicciotti et al., 2015) although thick supraglacial debris may keep ablation rates relatively low unless differential ablation occurs. Additional experiments were made (results not presented) that varied the debris supply to the entire glacier and the shape of the initial mass balance gradient before the impact of sub-debris melt was imposed. These experiments did not improve the model-data fit for the surface and emergence velocities and gave a weaker fit to the other data used to evaluate the model.
Mass balance is tuned to the extreme elevations of the accumulation area in the Western Cwm (>6,000 m a.s.l.) and therefore the model does not reliably simulate the geometries of Changri Nup and Changri Shar Glaciers as these accumulation areas have lower elevations (Figures 5e and 5f). As Changri Nup and Changri Shar Glaciers have not contributed mass to Khumbu Glacier since the LIA maximum the representation of these glaciers does not affect our results for evolution of Khumbu Glacier over the last 100 years. Our simulations of the evolution of debris-covered Khumbu Glacier captured the low-angle, reversed mass balance gradient that occurred across the ablation area because of the overall increase in debris thickness towards the terminus (Figure 5g). In 2015 CE, the mass balance of the entire glacier reached a minimum just below the base of the icefall with the most negative mass balance for the active glacier occurring 2-3 km further downglacier (Figures 5e and 5f).

Differential Ablation
Parameterization of differential ablation improved the ability of the model to reproduce the magnitude and distribution of ice volume loss more accurately than driving the model solely by forcing ΔT (Figure 3; Table 2). Positively skewed distributions of debris thickness with a form similar to that measured for Ngozumpa Glacier Nicholson et al., 2018) that accounted for the presence of ice cliffs and supraglacial ponds where debris thickness was at or close to zero, gave the closest model-data fit compared to negatively skewed, uniform and normal distributions without requiring additional ΔT forcing beyond that expected between the LIA and the present day (Rowan et al., 2015;Rowan, 2017). Indeed, additional ΔT forcing resulted in a weaker model-data fit (Table 2). These results illustrate the strong feedback whereby enhanced ablation forced by increased ΔT results in a reduction in sub-debris melt over decadal to centennial timescales as former englacial debris melts out to reduce ablation at the glacier surface. represents a larger set of processes that influence mass change than only a relative increase in melting at these hotspots, which include: enhancement of ablation above clean-ice values at ice cliffs, thermal erosion by supraglacial ponds (Miles, Willis, et al., 2018;Thompson et al., 2016;Watson et al., 2017), englacial and subglacial mass loss (Benn et al., 2017), warm englacial ice temperatures (Miles, Hubbard, et al., 2018), reduced snowfall, and a shorter or more variable monsoon season .
Ice cliffs and supraglacial ponds are discrete features on the surface of debris-covered tongues and their number and dimensions fluctuate both seasonally and spatially (Watson et al., 2017;Miles, Willis, et al., 2018). We did not include any possible feedbacks between the evolution of the glacier and changes in the number or distribution of such ablation hotspots as there are few observations of changes in ice cliff and pond magnitude and distribution (Iwata et al., 2000;Watson et al., 2017) although the extent and density of these features are expanding both up-and down-glacier (King, Turner, et al., 2020). More importantly, the processes controlling the evolution of these features are not yet understood, underlining the need for robust, long-term assessments and model development .

Reduced Mass Flux
Maximum surface lowering on Khumbu Glacier occurs above the zone where ice cliff densities are highest (Watson et al., 2017) and differential ablation is therefore highly unlikely to be the sole driver of the highest rates of ice loss. Instead, additional mass loss could result from decreased mass flux from the accumulation area to the debris-covered tongue, reflecting multidecadal changes in accumulation rates and dynamic detachment as a precursor to physical detachment of the active glacier from the debris-covered tongue. The disintegration of the upper ablation area could also be enhanced by the increasing input of meltwater from adjacent glaciers to the glacier bed, a process that has been identified for glaciers in the Indian Himalaya that display similarly anomalous dynamics (Singh et al., 2020). We observed a large volume of water in streams that flow from the Changri Nup palaeoconfluence into Khumbu Glacier (Miles et al., 2019). We suggest that dynamic detachment, driven by local glacier thinning has occurred close to the Changri Nup palaeoconfluence as suggested by shearing of the active glacier and the presence of basal ice at the glacier surface (Fushimi, 1977;Miles et al., 2021). Moreover, the decrease in ice flow through the upper ablation area corresponds to the downglacier limit of the area of enhanced glacier thinning, which is considered to be the approximate location of the termination of active ice flow (Nakawo et al., 1986;Rowan et al., 2015).
Dynamic detachment represents a precursor to physical detachment of the debris-covered tongue from the active glacier, which will result in a complete loss of ice from the middle section of the glacier as the active terminus recedes (Rowan et al., 2015). Simulations of Khumbu Glacier beyond the present day indicate that physical detachment will take place at the base of the icefall where mass balance is most negative and the glacier surface is not protected by supraglacial debris (Figure 4i). After the active glacier physically detaches from the debris-covered tongue, supraglacial debris will accumulate on the newly formed ablation area of the active glacier. Physical detachment is evident for other glaciers in Nepal, such as neighboring Lobuche Glacier in the Khumbu valley, Shalbachum Glacier in the Langtang valley, Nepal (Pellicciotti et al., 2015) and in the Canadian Rockies (Rippin et al., 2020) where the former glacier beds are exposed between the upper active sections and the former ablation areas. This is not yet the case for Khumbu Glacier, but the increase in observed δz between 1984-2000 and 2000-2015 suggests that mass loss here will rapidly lead to physical detachment.

Conclusions
Recent observations have demonstrated that debris-covered Himalayan glaciers are losing mass at similar rates to clean-ice glaciers, but the mechanisms driving this anomalous mass loss remained elusive. A recently identified and newly simulated set of processes that enhance ablation across the debris-covered area of such glaciers is shown to increase net mass loss by 29%-47% of the amount expected with a continuous debris layer. However, a low mass balance gradient and differential ablation across a discontinuous debris layer punctuated by ice cliffs and supraglacial ponds are not sufficient to explain the total observed loss of ice mass indicated by surface elevation change for Khumbu Glacier over a 31-year period . Instead, using a higher-order ice flow model with a novel representation of differential ablation we demonstrated that dynamic detachment of the upper active glacier from the heavily debris-covered tongue is also required to reproduce the recent evolution of Khumbu Glacier. Dynamic detachment occurred within the last 100 years along a shear plane underlying the upper ablation area, and the active glacier terminus is now located in the middle of the former ablation area. At the active terminus, velocities rapidly decline over a short length scale and the upper glacier overrides the thick slow-flowing ice contained within the debris-covered tongue. Dynamic detachment represents a tipping point in the evolution of Khumbu Glacier whereby the glacier lost 25% of its length (half of the ablation area) and 20% ice volume in a few decades. Furthermore, dynamic detachment is a precursor to a second tipping point-physical detachment of the accumulation area from the ablation area will occur close to the base of the Khumbu icefall and has already occurred for debris-covered glaciers elsewhere in the Himalaya. The accumulation area of Khumbu Glacier is the highest of any glacier on Earth, and we expect that these tipping points have already or will soon be crossed by many other debris-covered glaciers in the Himalaya and more widely.