Surface and subsurface hydrology of debris-covered Khumbu Glacier, Nepal, revealed by dye tracing

Surface and subsurface hydrology of debris-covered Khumbu Glacier, Nepal, dye tracing. While the supraglacial hydrology of debris-covered glaciers is relatively well studied, almost nothing is known about how water is transported beneath the glacier surface. Here, we report the results of sixteen ﬂuorescent dye tracing experiments conducted in April–May 2018 over the lowermost 7 km of the high-elevation, debris-covered Khumbu Glacier, Nepal, to characterise the glacier’s surface and subsurface drainage system. Dye breakthroughs indicated a likely highly sinuous and channelised subsurface hydrological system draining water from the upper part of the ablation area. This ﬂowpath was distinct from the linked chain of supraglacial ponds present along much of the glacier’s lower ablation area, through which water ﬂow was extremely slow ( ∼ 0.003 ms − 1 ), likely reﬂecting the study’s timing during the pre-monsoon period. Subsurface drainage pathways emerged at the glacier surface close to the terminus, and ﬂowed into small near-surface englacial reservoirs that typically delayed meltwater transit by several hours. We observed rapid pathway changes resulting from surface collapse, indicating a further distinctive aspect of the drainage of debris-covered glaciers. We conclude that the surface and subsurface drainage of Khumbu Glacier is both distinctive and dynamic, and argue that further investigation is needed to reﬁne the characterisation and test its regional applicability to better understand future Himalayan debris-covered glacier meltwater delivery to downstream areas.


Introduction
Meltwater from Himalayan glaciers and snow feeds some of Earth's largest river systems, influencing the supply of water to ∼1.4 billion people (Barnett et al., 2005;Bolch, 2017;Immerzeel et al., 2010). Approximately 30% of Himalayan glaciers have a supraglacial debris cover (Thakuri et al., 2014) that influences mass loss processes, the loci of maximum melt, and the volume of meltwater produced (Luckman et al., 2007;Østrem, 1959;Thompson et al., 2016). Where the debris cover exceeds several centimetres in depth and covers a considerable portion of the glacier's ablation area, it influences the hydrological system both on the surface and below it (Fyffe et al., 2019). These extensive debris covers produce a range of surface features not commonly found on clean-ice glaciers, such as supraglacial ponds and ice cliffs. * Corresponding author.
The supraglacial hydrology of debris-covered glaciers has received increasing attention in recent years. For example, it is welldocumented that the formation of proglacial moraine-dammed lakes, and the consequent presence of a local base-level, can be facilitated by the development, growth, and coalescence of supraglacial ponds into a chain of linked ponds Mertes et al., 2016;Sakai, 2012;Thompson et al., 2012). Melt rates are disproportionately high at pond margins due to the continued horizontal and vertical incision (Miles et al., 2016;Sakai et al., 2000), and recent studies have found that supraglacial ponds are expanding to cover an increasing proportion of the surface of debris-covered glaciers (Gardelle et al., 2012;Watson et al., 2016). This has implications for greater meltwater production and water storage, since ponds moderate diurnal glacier runoff (Irvine-Fynn et al., 2017).
On debris-covered glaciers, supraglacial streams do not tend to persist for long distances, instead incising to become englacial features Iwata et al., 1980). Supraglacial and shallow englacial conduits located towards the centre of debriscovered glaciers have been suggested to be discontinuous due to https://doi.org/10.1016/j.epsl.2019.02.020 0012-821X/© 2019 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). the variable surface topography, and are therefore most likely to transport meltwater between supraglacial ponds (Benn et al., 2017;Miles et al., 2017a;Narama et al., 2017;Thakuri et al., 2016). Longer-distance transport has been observed through perennial sub-marginal channels (Benn et al., 2017;Thompson et al., 2016), while more conduits open and transport water more effectively through the monsoonal melt season than through the dry season Hewitt et al., 1989;Miles et al., 2017aMiles et al., , 2017bSakai et al., 2000).
Aside from a limited number of speleological investigations (Benn et al., 2017;Gulley et al., 2009;Gulley and Benn, 2007;Narama et al., 2017), most work on Himalayan debris-covered glacier englacial drainage has been inferred from proxies. Even less is known about possible subglacial drainage networks; such systems have been deduced on the basis of methods such as proglacial sediment and water analyses (Haritashya et al., 2010;Hasnain and Thayyen, 1994) and remote-sensing observations of seasonal changes in glacier surface velocity (Benn et al., 2017;Copland et al., 2009;Kraaijenbrink et al., 2016;Quincey et al., 2009). Seasonal system evolution has, however, been determined directly from a few dye-tracing studies (Hasnain et al., 2001;Liu et al., 2018;Pottakkal et al., 2014). Fyffe et al. (2019) also conducted a dye tracing study on a debris-covered glacier in the Italian Alps (Miage Glacier), reporting that the continuous debris cover in the glacier's lower ablation area produced an inefficient subsurface drainage network, which joined a more efficient, channelised network draining water from the cleaner ice farther upglacier.
Consequently, while the intricacies of the englacial and subglacial drainage systems of Alpine and Arctic clean-ice glaciers have been well-studied and documented (e.g. reviews by Hubbard and Nienow, 1997;Fountain and Walder, 1998;Irvine-Fynn et al., 2011), there is limited knowledge of how water flows through and beneath debris-covered glaciers. Our aim herein is to use fluorescent dye tracing experiments to investigate the hydrological system of the debris-covered Khumbu Glacier, Nepal. Our specific research objectives are to: i) determine whether a subsurface hydrological system exists within and/or beneath Khumbu Glacier; ii) if such a system exists, determine its hydraulic characteristics, including likely flowpaths; and iii) further elucidate the nature of meltwater transport through the glacier's linked chain of supraglacial ponds.

Methods
Dye tracing experiments were carried out on Khumbu Glacier ( Fig. 1) during the 2018 pre-monsoon season. The glacier has a well-documented, expanding area of supraglacial ponds that in 2015 covered 3.2% of the glacier's 7.1 km 2 debris-covered area (Watson et al., 2016). Ponds are particularly prevalent along the glacier's eastern margin, where they have been coalescing and connecting hydrologically over recent years to form a linked chain (Irvine-Fynn et al., 2017;Watson et al., 2016). A large, perennial supraglacial channel originating in the upper clean-ice region of the glacier's ablation area (>9 km upglacier from the terminus) has been present since at least 2005 . Field observations in 2017 found that this channel progressively incises downglacier, becoming englacial just above the confluence with Changri Shar Glacier (Fig. 1). Khumbu Glacier's hydrological system receives an additional input of meltwater from Changri Shar and Changri Nup Glaciers (hereafter the Changri catchment), the proglacial streams of which coalesce and cascade down a steep gorge into, and likely beneath, Khumbu Glacier (Benn, pers. comm., 2018). An outburst flood from Changri Shar Glacier in 2017 entered Khumbu Glacier at this point and appeared to transit much of the glacier's length below the surface (E.S. . Little else is known about how water flows through Khumbu Glacier, but currently, only one dominant active supraglacial channel drains the glacier, forming the proglacial stream on the eastern margin of the terminus (Fig. 1B). An analysis of available satellite imagery confirms that this configuration has prevailed since at least the early 1990s.
Sixteen dye tracing experiments were undertaken between 27th April and 14th May 2018 across the lower ∼7 km of Khumbu's 9 km-long ablation area. Fluorescein dye was used due to its photo-degradation, resulting in minimal downstream impacts. Dye was injected into supraglacial streams or channels at selected locations (dye injection points (DIP * ), Fig. 1), in volumes of 1-150 ml according to the tracing distance (Table 1). More dye would have allowed clearer breakthrough curves, however injection volumes were restricted to minimise dye visibility beyond Khumbu Glacier due to the glacier and its proglacial stream being located in a National Park near popular trekking routes. In all cases, dye was detected using Turner Designs Cyclops-7 fluorometers (F * , Fig. 1) logging at one-minute intervals, located at strategic junctions in the supraglacial hydrological network. All fluorometers were shielded from direct sunlight using an inbuilt shade cap, and were fully submerged beneath the water. Images and a video of a dye injection are presented in the Supplementary Material ( Fig. S1 and Video S1, respectively). Dye tracing experiments were carried out at four different spatial scales over the glacier: i) terminus (300 m from the terminus), ii) short-range (within ∼600 m of the terminus), iii) long-range (∼7 km from the terminus), and iv) pond-based experiments (conducted in various supraglacial ponds near the terminus; Table 1). One terminus trace was carried out (dye injected at DIP1), with dye detected by a single fluorometer on the terminus stream (F0; Fig. 1C). Three pond-based traces were carried out just upglacier of this stream (DIP2, DIP3 and DIP4), also detected by the same fluorometer which was subsequently removed. Six short-range traces were conducted a short distance upglacier again (DIP3, DIP5 and DIP6), detected by two fluorometers located on two distinct inlets to the final supraglacial pond (F1 and F2). These fluorometers also detected dye returns from five long-range traces (DIP7 and DIP8). An additional fluorometer (F3) was positioned near the middle of the linked supraglacial pond chain for the later three long-range traces; one final pond-based trace was carried out above this fluorometer within the linked pond chain (DIP9).
The fluorometer used for the first six traces (one terminus trace and three pond-based traces at location F0, and the first two shortrange traces at location F1) could not be calibrated to record absolute concentration. This fluorometer was therefore replaced and not used again; these first traces are hereafter referred to as 'pretests' (PT * ) due to the extra data correction that was required to set the background to ∼0 parts per billion (ppb; an offset of +1,250 ppb was applied). This offset did not influence any subsequent analysis. All tests with the remaining three fluorometers, which all functioned and were calibrated correctly, are referred to as 'dye traces' (DT * ).
Measured fluorescein concentrations were corrected for water temperature, which was also recorded by the fluorometer loggers, as follows: where F r is the calculated fluorescence at the reference temperature, T r ; F s is the observed fluorescence at the time of reading the sample temperature, T s ; and n is the temperature coefficient for fluorescein (0.0036 • C −1 ) (Turner Designs, 2018). A small number of extreme data points, including negative values or readings that exceeded the maximum detection of the fluorometer, were removed and the remaining data were interpolated linearly. In all cases, such outliers comprised individual, isolated points, which we assume resulted from rare electronic disturbance. water seepage from/into supraglacial ponds beneath the debris layer, a shallow moulin and englacial reservoirs. The pond below DIP6 only appeared following fieldwork, and is inferred to have been located englacially during the field season (see Section 4.2). Supraglacial ponds and streams were mapped manually from the background PlanetScope Ortho Scene, captured during the field season on 24.04.2018 (Planet Team, 2017).

Table 1
Key data from the breakthrough curves and subsequent analysis of the successful pre-and dye traces (PT * and DT * , respectively). The long-range tests returned two distinct dye breakthrough curves (BTC 1 and 2).

Dye trace type
Fluorometer 0 for all pre-traces (PT * ) and Fluorometer 3 for DT7 only. ** Dye was a 41% solution. Each dye breakthrough curve was identified and dye transit time (t m , s) calculated from the time of injection to the maximum peak. Minimum transit velocity estimates (u m , m s −1 ) were calculated using the transit distance (x, m): for the terminus, pondbased, and short-range traces, the straight-line distance was used (Hubbard and Glasser, 2005;Seaberg et al., 1988). For the longrange traces, the straight-line distance was adapted slightly to follow the glacier's centreline. The dispersion coefficient (D, m 2 s −1 ) was calculated: where t j (s) represents t 1 and t 2 , the time of half the peak dye concentration on the rising and falling limbs of the breakthrough curve, respectively; t m (s) is the time to peak concentration, obtained by defining the above equation for j = 1 and 2 and solving iteratively for t m until a common value of D is obtained (Seaberg et al., 1988;Willis et al., 1990). Dispersivity (d, m) was then calculated: The fluorometer data were filtered using a 15 min window fast Fourier transform to remove daytime noise, which was likely a consequence of an increased and more variable suspended sediment concentration with a similar fluorescence wavelength to fluorescein during hours of greater discharge (Smart and Laidlaw, 1977). No concurrent discharge data were collected to corroborate this, which also precludes the calculation of dye recovery. The raw fluorometer time series are provided in the Supplementary Material (Fig. S2). Due to the limited duration of some dye breakthrough curves, the filtered data were not suitable for transit velocity or any subsequent calculations. The filtered curves were, however, used to verify the values of t j selected to calculate the dispersion coefficient, which were often difficult to discern due to the background noise. The difference in dispersivity between the filtered and unfiltered data was negligible (<1 m).
Discharge data are available from salt dilution tests during the 2017 pre-and post-monsoon seasons (Table S1), both from the upglacier supraglacial stream (DIP8) and at the terminus (by F0; Fig. 1). Whilst acknowledging that this dataset was collected the year before this study, exploration of the glacier's supraglacial drainage system in both years did not show any significant reconfiguration between 2017-2018. These discharge measurements are therefore presented as an approximation of the relative discharges of the glacier's supraglacial stream and outflow at the terminus, both before and after the monsoon.

Results
The full time series of the fluorescein concentrations recorded by three fluorometers, F0, F1 and F2, are shown in Fig. 2. Thirteen of the sixteen dye experiments gave successful returns, while the three that did not were not unexpected (Table 1): PT3 and PT4 likely merged with PT2, and the missing breakthrough of DT10 was probably caused by the fluorometers being removed before the dye emerged. The peak maxima of the long-range traces are substantially smaller (>6× smaller) than those of the short-range and pond-based traces. This is unsurprising given the considerable distances travelled for the long-range traces (∼70× greater than the short-range traces) while dye volumes were only scaled by a factor of ∼10.
A single terminus trace was conducted to obtain a measurement of the outlet stream velocity and dispersivity (Fig. 3). The breakthrough curve shows a short transit time with a narrow breakthrough curve width, implying low dispersivity. This is supported by the test data (Table 1), indicating a velocity of 0.24 m s −1 and a dispersivity of 3.0 m. The velocities of the short-range experiments are an order of magnitude lower than the terminus trace (Table 1; Fig. 4), but with similarly low dispersivity values. The breakthroughs recorded at both fluorometers from DIP3 and DIP5 tend to be faster and less dispersed than those from DIP6. These traces also show a difference in peak timing between F1 and F2, with breakthrough maxima at F2 occurring 11-35 min after F1. The final short-range traces (Fig. 4C) show a lower velocity from DIP3, and a much higher dispersivity for both injections (∼8.0 m increase for DIP3, and ∼6.0 m increase for DIP6 relative to the previous trace; Table 1).
The fluorometer returns from the pond-based traces are shown in Fig. 5. All four breakthroughs are uncertain and PT3 and PT4 are indeterminable having merged with PT2 (Fig. 5A). Where peaks seem apparent (particularly for DT7; Fig. 5B), the very small concentration range should be considered alongside the short transport distances for these traces (and the background stream fluorescence, also shown in Fig. S2). Velocities were calculated to the first major peak of PT2 (0.02 m s −1 ) and DT7 (0.003 m s −1 ; Table 1), but dispersivities could not be calculated due to the indistinct end points of all breakthrough curves. Fluorometer 3 was left running until the end of the measurement period but showed only background fluorescence values (see Fig. S2).
The breakthrough curves for the long-range experiments are presented in Fig. 6. All traces from DIP7 returned two distinct breakthroughs recorded at both fluorometers, separated by 10-20 h (Table 1; Fig. 6A, C, D). Two breakthroughs may also have been returned from the DIP8 trace ( Fig. 6B; clearer from the raw data) but the period between the breakthroughs was much shorter (<1 h). The first breakthrough consistently produced a peak of shorter duration, while the second breakthrough was slower (by 0.001-0.003 m s −1 over 5-7 km) and more dispersed. Average trace velocities (of both breakthroughs) are similar to the short-range traces (∼0.03 m s −1 ). The difference in breakthroughs between the fluorometers was present but less distinct compared to the short-range experiments: the F2 breakthrough maxima were, on average, 5.5 min later than for F1. The final DIP7 trace (Fig. 6D) occurred in approximately half the time of the first two traces. Summary velocity and dispersivity data are presented in Fig. 7. The differences between the extremes of transit distance are highlighted at either end of the plot: the shortest distance trace (DIP1) has the fastest velocity and a low dispersivity (156 m, 0.24 m s −1 and 3.0 m, respectively; Table 1); the longest transit distance (DIP8) shows low velocities and dispersivities (7,084 m; <0.018 m s −1 and <1.25 m, respectively). The repeat short-range experiments show a moderate decrease in velocity and significant increase in dispersivity over time (e.g. from 0.072 to 0.036 m s −1 and from 1.3 to 9.1 m for DIP3, respectively). The repeat long-range traces from DIP7 show a three-fold increase in velocity (e.g. from 0.017 to 0.052 m s −1 for the first breakthrough of DT3 and DT9, respectively).

Terminus trace
The single terminus experiment was conducted to aid the interpretation of all subsequent dye traces. Velocities ≥0.2 m s −1 and dispersivities ≤10 m have traditionally been interpreted to show drainage through an efficient, fast-flowing, integrated and channelised drainage system (Burkimsher, 1983; Seaberg et al., 1988; Willis et al., 2012Willis et al., , 1990. Both the moderately high velocity (0.24 m s −1 ) and low dispersivity (3.0 m) of this test confirm our observations that this large, single supraglacial channel evacuates meltwater rapidly from a small supraglacial pond into the proglacial stream (approximate discharge 1-2 m 3 s −1 , Table  S1; image of injection and channel shown in Fig. S1A). This dye trace therefore provides a reference for the maximum efficiency we might expect within the system (Fig. 7).

Short-range traces
As the only visible exit for meltwater on Khumbu Glacier is from the supraglacial system at the terminus, the short-range experiments were carried out to characterise this sector, which is also the end of the long-range drainage system. The velocities of 0.02-0.07 m s −1 indicate slow transit, with consistently faster and less dispersed transport from DIP3 than from DIP6. Direct field observations indicated that DIP3 was located at the head of a short length of supraglacial channel draining a supraglacial pond, which flowed into a shallow moulin feeding a larger, shallow englacial reservoir ( Fig. 1C; Fig. S1C-D). This reservoir was visible through a surface fracture, flowing into the small supraglacial pond above F1. After injection, dye could be followed down the stream and viewed in the englacial reservoir through a surface fracture, only slowing in the supraglacial pond immediately above F1. To produce slow velocities but reasonably low dispersivities, we infer that the reservoir and supraglacial pond provided temporary storage but transit was relatively undisturbed by flow complexities, such as turbulence and/or eddies, limiting dispersion. In addition to the main outflow past F1, dye was observed to leave this pond by seepage under the debris along the eastern pond margin towards the next downstream pond (Fig. 1C), a process that has previously been suggested but not observed (Irvine-Fynn et al., 2017). The dye was   Table 1). Each dye injection point (DIP * ) has a separate symbol: for the long-range tests, the larger symbol indicates the first dye breakthrough, and the smaller symbol the second dye breakthrough. Where the dye breakthrough was recorded by two fluorometers, values were averaged to show the difference between the dye breakthroughs more clearly. The dotted lines and shaded areas indicate values that determine when a drainage system is traditionally considered to be channelised rather than distributed (velocity ≥0.2 m s −1 and dispersivity ≤10 m, Burkimsher, 1983;Seaberg et al., 1988;Willis et al., 1990). well diffused by the slow, constant flow through the supraglacial pond; it is assumed that our breakthroughs from the main pond outlet, on which F1 was located, are representative of all pond outflows. The multiple visible and concealed outflows from this one small pond highlight the distinctive near-surface complexity of debris-covered glacier hydrology.
The slower velocities and greater dispersivities from DIP6 compared to DIP3 imply greater storage in the 100 m upstream of DIP3 (transit time of hours rather than minutes). DIP6 was located at the downglacier end of a large supraglacial pond, at the end of the linked pond chain. After injection, dye moved slowly towards the down-flow tip of the pond, exiting the pond beneath the debris (image in Fig. S1B). How water was transported between DIP6 and DIP5 is not known, but we infer that this flowpath included an additional, larger, englacial reservoir(s) due to: i) an absence of surface water in the vicinity; ii) the far slower speed and greater dispersion of the DIP6 breakthrough curve relative to DIP5 and DIP3; and iii) the only visible input to the pond at DIP5 being a small inlet, also beneath the debris layer. Satellite imagery shows the formation of a supraglacial pond between DIP6 and DIP5 shortly following the field experiments (Fig. 1C), which may have resulted from the flooding and/or surface collapse of an englacial reservoir.
The difference in the time of breakthrough at the two fluorometers was consistent for all the short-range experiments, revealing a slight variation in drainage between F1 and F2 downstream of DIP3. The drainage network leading to F1 was noted above, but the system feeding F2 is unknown, despite exploration. We infer that a portion of the flow is diverted from the F1 path within the englacial reservoir, and follows a less direct route before emerging as the F2 stream. Due to the dispersivity variations between the short-range tests, we are unable to confirm whether the network comprises an additional englacial reservoir, a sinuous conduit or a small distributed channel network.
The repeat short-range traces reveal a decrease in system efficiency over the field season, with the final experiment recording 50% of the velocity measured in previous tests, and dispersivity values 140-700% greater (Table 1; Fig. 7). Although the limited number of traces should be acknowledged, these results do align with changes observed in the field. Before the final repeat trace on 8th May, we found that the moulin below DIP3 had collapsed and re-routed flow into a shallower gradient moulin ∼5 m upstream.
In this section, the stream had a greater discharge and faster flow than had been previously observed. As the dye still reached F2, the system changes did not cut this drainage route off, confirm-ing that the flow divergence was downstream of these moulins. The breakthrough curves, particularly from DIP3, display multiple peaks (Fig. 4), suggesting that this new pathway routed water less efficiently through multiple conduits into the englacial reservoir, despite the greater discharge into the moulin. Debris in the conduit may have contributed to the higher dispersivity (Gulley et al., 2014), but the newly-formed pathway may simply have been more convoluted, possibly due to it exploiting voids or pre-existing weaknesses in the ice. Dye may also have been delayed and dispersed in the englacial reservoir due to its higher entry point. Our observations of these pathways and their rapid changes in hydraulic transfer displays a further hydrologic feature that may be distinctive to debris-covered glaciers: re-routing occurring due to surface collapse, which is extremely prevalent on debris-covered glaciers due to their highly spatially variable rates of surface lowering (Benn et al., 2017;E.S. Miles et al., 2017a.

Pond-based traces
The dye for the pond-based experiments at the terminus (PT2 to PT4; Fig. 5A) was injected over 70 min between the three DIPs, but there is only one clear breakthrough maximum (from PT2), yielding a low transit velocity (0.02 m s −1 ). The half-life of fluorescein depends on a large number of factors, but is short when exposed directly to bright sunlight (Smart and Laidlaw, 1977): small sample traces suggest complete photo-degradation within a few hours (Turner Designs, 2018). When dispersed in more turbid water, through which light is attenuated, the half-life will be correspondingly longer. We suggest that dye transit was so slow through the ponds that the separate dye injections merged and much of the dye photo-degraded in the large ice-free pond immediately beyond F1/F2. The same conclusion can be drawn from the other pond-based test (DT7 from DIP9; Fig. 5B), the indeterminate breakthrough of which produced a velocity of 0.003 m s −1 to the largest first peak, ∼40 h after injection (cf. Fig. S2). This is an order of magnitude slower than terminus pond traces and two orders of magnitude slower than the terminus trace (Fig. 3). The ratio of throughflowing discharge to pond volume is therefore inferred to be exceptionally slow prior to the monsoon melt season, supporting previous findings (Irvine-Fynn et al., 2017). No dye was detected at F1 or F2; it may have become very well dispersed within the pond and/or been largely destroyed by photo-degradation between F3 and DIP4.

Long-range traces
Although of low concentration, the long-range tests yielded breakthrough curves (Fig. 6), confirming that meltwater is transported from the upper ablation area to the terminus of Khumbu Glacier, exiting the glacier at the surface rather than (entirely) being lost to groundwater. The slow velocities for the long-range tests (mean of ∼0.03 m s −1 ) and unexpectedly low associated dispersivities (all <1.5 m; Table 1, Fig. 7) indicate minimal long-term storage and/or eddying within the transport path, such as might be found in supraglacial or englacial ponds. The absence of long-range trace breakthroughs at F3 further indicates that this subsurface drainage system does not link into the supraglacial pond chain.
Compared to previous dye tracing studies on other debriscovered and clean-ice glaciers that inferred flow through an efficient, channelised drainage network partly from fast throughflow velocities (Burkimsher, 1983;Fountain, 1993;Hasnain et al., 2001;Nienow et al., 1998;Pottakkal et al., 2014;Schuler et al., 2004;Seaberg et al., 1988), our trace velocities are one, and in some cases two, orders of magnitude lower. Our results are more akin, for example, to those from the similar-sized but debrisfree Midtdalsbreen, Norway, which were interpreted in terms of drainage through a linked-cavity system (Willis et al., 1990). However, these networks produced markedly higher dispersivity values (10.0-71.1 m) than we record at Khumbu Glacier. It therefore appears that at least part of the subsurface drainage network at Khumbu Glacier differs from each of the two forms previously reported from clean-ice glaciers: in showing a relatively low velocity and low dispersivity it is neither completely channelised/efficient (high velocity and low dispersivity) nor completely distributed/inefficient (low velocity and high dispersivity).
We propose that Khumbu Glacier's drainage system across the ablation area initiates supraglacially in the clean-ice region beneath the ice fall above DIP8. Gulley et al. (2009) mapped several supraglacial channels incising into 'cut-and-closure' englacial conduits in this region in 2005, the pathways of which are very similar to parts of the sinuous supraglacial channel network we observed in 2017-18 ( Fig. 1B; image in Fig. S1F). We suggest that rapid surface melt has caused these englacial conduits to become re-exposed at the surface, as has been observed on other debriscovered glaciers (Miles et al., 2017a). Indeed, this allowed us to hike along the progressive incision of the supraglacial stream in 2017 to where the stream now disappears to become englacial, ∼300 m upglacier of the large meltwater input from the Changri catchment. At this point, multiple relict conduits were visible above the active channel, supporting the interpretation of an incising cut-and-closure conduit. The location of the stream submergence to become englacial, and the similarities of the trace breakthroughs from DIP7 and DIP8, lead us to infer that this channel likely continues to incise downwards to join the stream input from the Changri catchment. The Changri input was observed to reach and follow the bed of Khumbu Glacier in 2006 (Benn, pers. comm., 2018), and we expect that this is still the case due to the substantial discharge of this torrent during our field observations both before and after the monsoon season ( Fig. S1E shows an image of dye injection at DIP7). This has allowed the subglacial stream to adopt a stable position at the bed in this section of the glacier (Benn, pers. comm., 2018).
Khumbu Glacier, like most glaciers in this region, has a largely impermeable terminal moraine that has previously been noted to provide a high local base-level for the glacier's englacial drainage . On the basis of the dye breakthroughs at the surface near the terminus and only one supraglacial channel draining the glacier, we suggest that the high hydrological baselevel for the lower glacier has prevented the subglacial channel continuing along the bed further downglacier, resulting in an 'uprouting' of the system back to the surface (of no more than a few tens of metres, Gades et al., 2000). This would produce an englacial depth limit to the cut-and-closure mechanism in the lower part of the ablation area -possibly following the coldtemperate ice boundary inferred by K.E. . Between the subglacial-englacial transition and the reappearance at the surface, the drainage continues beneath the surface, bypassing the supraglacial and near-surface hydrological network in the lower ablation area. The subsurface drainage system downglacier of the Changri Shar confluence may therefore be a continuation of Khumbu's sinuous, upglacier cut-and-closure channel. Speleological investigations on the neighbouring Ngozumpa Glacier have shown that cutand-closure englacial conduits can produce highly sinuous channels (Benn et al., 2017;Gulley et al., 2009). Indeed, this may also be prevalent within Khumbu Glacier, due to the glacier's debriscovered tongue having a very low surface gradient (Fig. 1B). These may combine to produce a low hydraulic gradient and thus encourage meandering. A sinuous system could potentially result in the actual transit velocities being up to an order of magnitude higher (using the actual, rather than straight-line distance): such velocities may even be similar to those used to infer channelised, rather than distributed, drainage on clean-ice glaciers (e.g. Nienow et al., 1998). They would also correspond better to the associated trace dispersivities indicating undispersed dye transit through the drainage system (Fig. 7). However, such velocities would still be lower than those observed on the debris-covered Miage Glacier, which has a much steeper ablation area, potentially encouraging faster flow (Fyffe et al., 2019), and the overall transit through Khumbu Glacier is still slow and inefficient. We suggest there is a short section beneath the Changri Shar confluence that is not sinuous (Benn, pers. comm., 2018), having adapted to the Changri catchment input. A notable, straight, surface depression follows the western glacier margin for a kilometre or so, beginning immediately below the confluence, which may be similar to the submarginal channels observed on Ngozumpa Glacier (Benn et al., 2017;Thompson et al., 2016).
The dual breakthroughs for our long-range traces are most likely caused by a division of the drainage system (Burkimsher, 1983;Willis et al., 1990). Given that we interpret this flow as predominantly englacial, these dual pathways could be similar to the englacial side passageways reported on Ngozumpa Glacier by Benn et al. (2017).
Further complexity in Khumbu's subsurface drainage system is suggested by the 2017 pre-monsoon salt dilution discharge measurements (Methods ; Table S1). At 13:00 on 24th May 2017, the discharge of the supraglacial stream near DIP8 was ∼1.5 m 3 s −1 ; two days later, also at 13:00, the discharge of the proglacial stream near F0 was ∼1.3 m 3 s −1 . Comparable measurements were found in the post-monsoon season: on 22 nd October 2017, the supraglacial stream discharge was ∼2.4 m 3 s −1 ; two days later at 13:00, the proglacial stream discharge was ∼2.6 m 3 s −1 . The similarity of these values indicates that the principal supraglacial stream alone could contribute almost all of the glacier's proglacial discharge, both pre-and post-monsoon. Khumbu's full subsurface meltwater discharge, including inputs from the Changri catchment, may not therefore emerge from the glacier's terminal portal; yet, no notable upwelling emerges from the glacier's terminal moraine. It is possible that the 2017 discharge measurements do not reflect longer-term conditions. Alternatively, if these discharge measurements are representative more generally, then a sizeable proportion of Khumbu's subsurface meltwater must emerge elsewhere from the system. We note that this cannot hold for all the glacier's subsurface flow because the long-range dye traces gave successful breakthroughs at the surface near the terminus. The lost component could be stored within the glacier, flow as groundwater that emerges farther down-valley than the terminal moraine, or emerge diffusely across the outer slope of the terminal moraine at a rate that is insufficient to overcome local evaporation. All of these processes deserve further investigation. The discharge measurements indicate further that the linked supraglacial pond chain provides a relatively small proportion of the glacier's discharge.
The bypassing of the supraglacial hydrological system by the subsurface system for much of the lower ablation area is similar to the pathway observed for perennial sub-marginal conduits on Ngozumpa Glacier, which also route back to the surface very close to the terminus into the proglacial Spillway Lake (Benn et al., 2017). Khumbu Glacier's subsurface drainage is routed to the surface between the end of the linked pond chain (F3) and F1/2. The time difference between the fluorometer breakthroughs is similar to that of the short-range traces (F2 being consistently later), confirming that the flow into F1/2 only diverges very close to the fluorometers and that the subsurface drainage system joins the short-range drainage system we observe, perhaps near/into the englacial reservoir visible from the surface. However, the upglacier subsurface drainage system may not be perpetually separate from the supraglacial system: evidence for the periodic drainage of perched ponds into the englacial network has been observed on both Khumbu and Ngozumpa Glaciers, as well as at other debriscovered glaciers (Benn et al., 2017Miles et al., 2017a). Further, E.S.  observed a lake outburst event that propagated through Khumbu Glacier in 2017, suggesting an overflow pathway created by the flood waters emerged from the subsurface drainage system into the lower supraglacial pond chain.
The breakthrough for the final trace from DIP7 (DT9) occurred notably sooner than the first two traces (DT3 and DT8), with the first dye breakthrough showing a velocity over three times greater ( Fig. 6D; Fig. 7). On the nights of the 1st and 8th May (5 and 6 days before DT4 and DT9, respectively), there were large snowfall events. We therefore interpret that the faster velocities for DT4 and DT9 may have been a short-term system response to greater meltwater inputs. Snow melt from the second event began around midday on the 9th May on the glacier itself; the snow cover may have reduced surface ablation, resulting in the slower velocity of DT8 than DT4. The influence on the subsurface drainage system would have been greater several days later due to the delayed melting of the snow at higher elevations, contributing to the greater velocity of DT9.

Summary and conclusions
To our knowledge, this study reports the first successful dye tracing experiments at a debris-covered glacier in the Nepal Himalaya and the first dye tracing-based investigation exploring the intricacies of debris-covered glacier hydrology, both at and beneath the surface. We conducted sixteen dye tracing experiments on Khumbu Glacier, Nepal, which reveal previously unknown features and complexities in the surface and subsurface hydrology of a debris-covered glacier, many of which differ from current notions of clean-ice glacier drainage. We highlight the following conclusions: • A likely highly sinuous and channelised subsurface drainage system exists at Khumbu Glacier, flowing for some distance along the glacier's bed below the confluence with Changri Shar Glacier. The system does not appear to involve long-term storage, and may have the potential to transport water rapidly, particularly after heavy precipitation events.
• Flow through the linked chain of supraglacial ponds along the eastern margin of Khumbu Glacier is extremely slow (velocity ∼0.003 m s −1 ) and comprises only a small proportion of total flow from the glacier, during the pre-monsoon season.
• Subsurface flow is a parallel system that bypasses the linked supraglacial pond chain. It is ultimately re-routed back to the surface close to the terminus, where it joins the supraglacial system feeding the proglacial stream. A proportion of the meltwater inputs to the subsurface flow do not appear to reemerge by this route, but the ultimate destination remains unclear.
• We observe pathway changes in the short-range linked supraglacial and shallow englacial drainage routes near the terminus, triggered by channel collapse induced by continued topographic evolution of the surface.
While noting that the number and length of the traces was limited due to the difficulties of working in such an environment and the complexity of the system, we suggest there is great scope for future investigations of both the surface and subsurface hydrological systems of Khumbu Glacier and other debris-covered systems in the greater Himalaya and beyond. Our results are influenced by the timing of our experiments early in the melt season: much more could be learned by repeating these tests during or after the monsoon season when the subsurface hydrology may become more developed after sustained large inputs to the system. Our dye tracing experiments have revealed previously unknown features of the subsurface hydrology of a Himalayan debris-covered glacier, which appears to provide the principal drainage network for much of the glacier's ablation area. Combined with the primarily pond-based supraglacial drainage in the lower ablation area, Khumbu Glacier has a heterogeneous and dynamic hydrological system. Further investigation is thus required to improve understanding and prediction of future meltwater production from, and transit through, Himalayan debris-covered glaciers. Such knowledge is essential given the substantial population relying upon Himalayan snow and ice melt.

Author contributions
KM and BH designed the study and carried out the fieldwork with assistance from DQ and EM. KM processed and analysed the data, and wrote the manuscript. All authors contributed to interpretations and manuscript editing.