Variations in near‐surface debris temperature through the summer monsoon on Khumbu Glacier, Nepal Himalaya

Debris surface temperature is a function of debris characteristics and energy fluxes at the debris surface. However, spatial and temporal variability in debris surface temperature, and the debris properties that control it, are poorly constrained. Here, near‐surface debris temperature (Ts) is reported for 16 sites across the lower elevations of Khumbu Glacier, Nepal Himalaya, for the 2014 monsoon season. The debris layer at all sites was ≥1 m thick. We confirm the occurrence of temporal and spatial variability in Ts over a 67‐day period and investigate its controls. Ts was found to exhibit marked temporal fluctuations on diurnal, short‐term (1–8 days) and seasonal timescales. Over the study period, two distinct diurnal patterns in Ts were identified that varied in timing, daily amplitude and maximum temperature; days in the latter half of the study period (after Day of Year 176) exhibited a lower diurnal amplitude (mean = 23°C) and reduced maximum temperatures. Days with lower amplitude and minimum Ts were concurrent with periods of increased seasonal variability in on‐glacier air temperature and incoming shortwave radiation, with the increased frequency of these periods attributed to increasing cloud cover as the monsoon progressed. Spatial variability in Ts was manifested in variability of diurnal amplitude and maximum Ts of 7°C to 47°C between sites. Local slope, debris clast size and lithology were identified as the most important drivers of spatial variability in Ts, with inclusion of these three variables in the stepwise general linear models resulting in R2 ≥0.89 for six out of the seven sites. The complexity of surface energy fluxes and their influence on Ts highlight that assuming a simplified relationship between air temperature and debris surface temperature in glacier melt models, and a direct relationship between debris surface temperature and debris thickness for calculating supraglacial debris thickness, should be undertaken with caution. © 2018 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd.


Introduction
Debris-covered glaciers exhibit a continuous mantle of rock debris over the full width of at least some of their ablation zone (Kirkbride, 2011). These glaciers are common in mountainous regions across the world, including in the European Alps (Mihalcea et al., 2006), Andes (Glasser et al., 2016), Southern Alps of New Zealand (Kirkbride, 2000) and the Himalaya (Scherler et al., 2011). The presence of a supraglacial debris layer influences glacier ablation, acting as a thermal buffer between the atmosphere and glacier ice surface, and modifying the energy available for melt (Jansson and Fredin, 2002;Kirkbride, 2000). The extent to which a supraglacial debris layer controls ablation is primarily dependent on the thickness of the debris layer (Clark et al., 1994;Mattson, 2000;Østrem, 1959). While a thin layer of debris below a critical thickness causes an increase in ablation due to a reduction of the surface albedo (Nakawo and Rana 1999), ablation exponentially decreases with increasing debris thickness above a critical thickness, as the debris layer inhibits glacier melting by attenuating and reducing thermal energy transfer to the underlying ice surface (Brock et al., 2010;Mihalcea et al., 2008a;Nicholson and Benn, 2006;Reid et al., 2012).
Supraglacial debris surface temperature is a function of the surface energy balance and modulates heat transfer through the debris layer (Nakawo and Young, 1981). Therefore, debris surface temperature can provide useful insight into the extent to which debris properties affect energy transfer at the surface of and through a debris layer. To date, little focus has been given to the influence of spatial and temporal variability in surface temperature across supraglacial debris layers, which can be affected by incoming energy fluxes and debris properties including albedo, surface roughness, sediment porosity, and moisture content (Evatt et al., 2015;Reznichenko et al., 2010;Rounce et al., 2015). Nicholson and Benn (2013) highlighted the occurrence of spatial and temporal variability in supraglacial debris properties and their influence on surface temperature and temperature gradients through the debris layer, and therefore glacier mass balance. However, many of the previous studies concerned with the measurement of debris surface temperature on glaciers have had limited spatial or temporal extent. For example, Nakawo and Young (1982) measured debris surface temperature at six plots over a 48 h period, while Nicholson and Benn (2006) measured debris surface temperature at a maximum of 11 plots on one glacier, but only for a maximum period of 11 days. Steiner and Pellicciotti (2015) presented one of the most extensive debris surface temperature data sets to date, from 13 locations over three ablation seasons on Lirung Glacier, Nepal. However, the study focused on describing the relationship between air temperature (T a ) and debris surface temperature rather than exploring spatial variability in debris surface temperature. Moreover, Steiner and Pellicciotti (2015) did not state the thickness of the debris layer underlying each of the sensors measuring debris surface temperature, an important factor in the consideration of spatiotemporal variability in debris surface temperature and the influence of underlying ice (cf. Nicholson and Benn, 2006). Consequently, the nature of and controls on debris surface temperature variability remain poorly constrained in glacial environments.
Conversely, ground surface temperature variability has been relatively well studied in other cold region environments (Gubler et al., 2011;Guglielmin, 2006;Romanovsky and Osterkamp, 2000) where significant spatial variation arises from localised changes in surface properties and environmental conditions. These studies have concluded that such variability influences the accuracy of surface energy balance modelling in these environments. We therefore contend that such variability may also be applicable to numerical modelling of debriscovered ice ablation and the response of these glaciers to climate change.
The importance of studies of debris surface temperature on debris-covered glaciers is manifested in the recent application of temperature-index models to debris-covered glaciers, which determine debris surface temperature from T a (Carenzo et al., 2016). Furthermore, debris surface temperature has previously been used to determine debris layer thickness through two approaches: the use of an empirical relationship between debris surface temperature and debris layer thickness, based on field data (Mihalcea et al., 2008a(Mihalcea et al., , 2008bMinora et al., 2015); and a surface energy balance approach also using debris surface temperature (Foster et al., 2012;Rounce and McKinney, 2014). Currently, neither approach has been considered robust, as the empirical approach is only applicable for debris layers thinner than 0.5 m (Mihalcea et al., 2008a) and the energy balance approaches exclude consideration of spatially variable debris properties such as albedo, surface roughness or moisture content that will affect energy exchange and therefore surface temperature at the debris surface (Collier et al., 2014;Evatt et al., 2015;Rounce et al., 2015). To understand the validity of these methods, and discern how to develop them further, confirmation of both the spatiotemporal regime of debris surface temperature and its controls is needed.
Considering these shortcomings, here we aimed to characterise the spatial and temporal variability in debris surface temperature on a debris-covered glacier using data collected from temperature sensors located in the debris near-surface and distributed over the lower ablation area of Khumbu Glacier, Nepal, in areas of thick (≥1 m) debris cover. The primary objectives of the study were to: (i) examine the temporal and spatial variation of debris surface temperature during an ablation season; and (ii) determine the controlling factors underlying variations in debris surface temperature.
The ablation area is almost entirely debris covered below 5400 m a.s.l., with the debris layer >2 m thick in places (Gades et al., 2000). The debris-covered ablation area displays a wide range of clast sizes comprising of granitic and schistose lithologies derived from the surrounding hillslopes (Iwata et al., 1980;Nuimura et al., 2011). The debris-covered area is topographically complex and dynamic being characterised by an undulant surface punctuated by numerous supraglacial ponds and associated ice cliffs, which changes over seasonal and interannual timescales (Nuimura et al., 2011;Watson et al., 2016). The more stable, lowermost region of the ablation area shows the early stages of soil formation and is partially vegetated .

Central Himalayan climate
The South Asian Summer Monsoon (hereafter, 'the monsoon') dominates the climate of the Khumbu Glacier catchment, and the Central Himalaya. The highest annual air temperatures occur between May and October (Ageta, 1976;Nayava, 1974) and~80% of precipitation falls between June and September (Bookhagen and Burbank, 2010). During the onset and progression of the monsoon season, high pressure over the Tibetan Plateau results in an increased temperature and pressure gradient southward towards the Indian subcontinent (Yasunari, 1976). This pressure gradient produces seasonally variable wind patterns in the Central Himalaya region and localised synoptic weather systems are dominated by mountain and valley winds, which vary on sub-diurnal timescales (Bollasina et al., 2002). As the monsoon season progresses, increases in regional precipitation frequency, air temperature, relative humidity and incoming longwave radiation occur, and are coupled with a decrease in shortwave radiation attributed to increasing cloud cover (Salerno et al., 2015;Shea et al., 2015).

Data acquisition
Near-surface debris temperature Temperature sensors Near-surface debris temperature (T s ) was measured as a robust proxy for true debris surface temperature using Maxim iButton ™ Thermochron temperature sensors (model number DS1921G: http://datasheets.maximintegrated.com/en/ds/ DS1921G.pdf), which record instantaneous temperature from -30 to +70°C with a manufacturer-stated accuracy of ±1.0°C. iButton sensors were chosen due to their low cost, reliability (Hubbart et al., 2005) and previous successful applications in a number of environmental settings including permafrost landscapes (Gubler et al., 2011). Gemini Tiny Tag ™ Plus2 data loggers (model number TGP-4520) with encapsulated thermistor probes were used for sensor calibration prior to fieldwork and have a manufacturer-stated accuracy of ±0.4°C. The iButtons were placed in waterproof polycarbonate plastic containers to protect from water damage following the method of Gubler et al. (2011). The effect of polycarbonate plastic waterproof casing on temperatures recorded was tested in laboratory conditions prior to fieldwork. In laboratory conditions, temperatures recorded by contained and uncontained iButtons in the same environments varied by <2°C, and more typically by ≤0.5°C, which is within the manufacturer's stated accuracy (see Supplementary information; Figure S1).
Field experiment design Near-surface debris temperature (T s ) was measured at hourly intervals at 16 sites between the 21 May and 29 July 2014 (Day of Year (DOY) 141 and 210). The first 48 h of each T s timeseries were discarded to allow the sensors to equilibrate with local conditions. For all sites, iButtons were placed in the immediate near-surface of the debris layer, typically between 0.01 and 0.05 m below the surface, using a single layer of clasts of representative size for each site from the immediate surrounding area as a shield from direct solar radiation as is common practice in ground surface temperature studies (Apaloo et al., 2012;Gisnås et al., 2014). Using a handheld Garmin 64 GPS, the iButton temperature sensors were distributed across the lowermost 2 km 2 of Khumbu Glacier's ablation area in a gridded pattern (Figure 1(c)). The elevation of sensor sites varied across the study area by 49 m between 4903 m a.s.l. and 4952 m a.s.l. (±3 m due to vertical accuracy of the handheld GPS) and each site had a unique combination of site characteristics, varying in slope, aspect, elevation, clast size, sorting, roundness, and clast lithology (Table I; see also the section 'Ancillary data').
To allow examination of the influence of additional debris layer properties and incoming energy fluxes on T s other than debris layer thickness, all iButton temperature sensors were installed in locations where the debris layer had a thickness of ≥1 m where the effect of cold propagation from underlying ice on T s is insignificant (Foster et al., 2012;Nicholson and Benn, 2006). Debris thickness was established by excavating the debris layer adjacent to the iButton location to a depth of 1 m; if no ice was present, debris thickness was reported as >1 m. At each site, a textural description of the debris was made, and digital photographs were taken before and after emplacement of the sensors (Figure 2). The iButton temperature sensors at Sites 7 to 13 were placed within a 90 × 90 m area to investigate variability in T s across an area typical of the resolution of remotely sensed thermal satellite data (e.g. ASTER) often used for supraglacial debris thickness mapping.
On retrieval of the iButton temperature sensors at the end of the monsoon season, comparison with the initial site photographs was used to evaluate any surface change at each site. For all 16 sites reported, the debris showed little or no disruption after sensor installation, and none of the temperature sensors were exposed at the time of collection. A further 42 iButton sensors were installed on the glacier surface but, due to topographic change during the monsoon season, they could neither be located or retrieved. Despite following standard methods for measuring ground surface temperature (Apaloo et al., 2012;Gisnås et al., 2014), placing clasts on the contained iButtons to shield them from direct incoming shortwave radiation created an additional source of uncertainty in the 16 retrieved T s data. Consequently, our measurements of T s do not necessarily reflect absolute debris surface temperature (Conway and Rasmussen, 2000) as the emplacement of sensors beneath clasts may mean that the sensors record temperature below rather than at the debris surface. Without detailed knowledge of the specific thermal properties of the debris at each site, more accurate assessment of the uncertainty between near-surface and true surface temperature is challenging. However, here we assumed our T s data were sound proxies for absolute T s . To identify any data which were likely to be less representative of true surface temperature, uncertainty at each site was estimated using the diurnallyaveraged temperature gradient calculated through a debris layer by Nicholson and Benn (2006) from data collected on nearby Ngozumpa Glacier of -10.5°C m -1 , and mean clast size for each site. These uncertainties ranged from 0.03°C to 4.39°C (Table I). Temperature metrics (mean T s , maximum T s , minimum T s and T s amplitude) were also regressed against estimated sensor depth. No significant relationship was identified meaning T s variability between sites cannot be attributed directly to sensor depth. Consequently, sites at which the calculated near-surface to surface temperature difference was greater than 0.5°C (the assessed uncertainty in our iButton sensor data) were considered to be less reliable in reflecting absolute surface temperature (Sites 1, 2, 9, 11 and 13), and were therefore either noted or omitted from subsequent analyses to avoid potential influence of misrepresentative data.
Mean clast size was considered a proxy for sensor burial depth, although it is probable that clasts covering the sensors were smaller than the mean clast size as a bias towards the smaller clasts would have occurred on emplacement. It is therefore expected the uncertainty calculated using mean clast size overestimates burial depth, and consequently the uncertainty in temperature with depth is less than estimated. However, this method of uncertainty calculation does not include consideration of diurnal variability in temperature gradient through the debris layers, which may cause mean temperature differences calculated here to be larger at certain times of day (as observed by Nicholson and Benn, 2006). The influence of this diurnal variability on results is discussed in the section 'Spatial variability in near-surface debris temperature'.

Ancillary data
Clast size and lithology Clast size at each site was estimated from 18.0 Mpix digital site photographs acquired using a Canon 550D camera and processed in ImageJ, v. 1.49 (Rasband, 2008), following the method outlined by Igathinathane et al. (2009). At all sites, images covered approximately 1 m 2 and a known scale in each photograph was used to define the metre: pixel ratio. Clasts were selected using a random sampling method. For each site photo, every clast identified was assigned a number, and a random number generator was used to subsample 25 clasts for measurement within ImageJ. Assuming from the 2D imagery that the long and intermediate clast axes were visible, the intermediate axis length was retrieved and a mean representative clast size for each site calculated (Table I). Where the intermediate axis of a clast was larger than the photo (Sites 9 and 13) the maximum length measurable from the scaled image was used.
Clast lithology was determined in the field using clast size, colour and mineral composition. Two major lithologies were identified; granite and schist. The dominant lithology at each site (Table I) was determined by manually classifying the lithology of all clasts in each of the site photographs in ImageJ and then calculating the percentage of granite for each site (Solano et al., 2016).
Local meteorological data Meteorological data were collected at four locations: on the debris-covered glacier surface of Khumbu Glacier at an elevation of 4950 m a.s.l. (Figure 1 Off-glacier air temperature (T aP ) was recorded at hourly intervals 2 m above the ground surface, using an artificially ventilated LSI-Lastem DMA 570 sensor (accuracy ±0.2°C) at the Pyramid Observatory. On-glacier air temperature (T aG ) was recorded at 30 min intervals in a location with schistose debris lithology (Figure 1(c)) using a Gemini Tiny Tag ™ Plus2 data logger (model number TGP-4520) and thermistor probe with a stated accuracy of ±0.2°C. The on-glacier thermistor probe was placed in a naturally aspirated radiation shield mounted on a tripod 1 m above the debris surface. T aG was resampled to give hourly values corresponding to the resolution of the T s data. Incoming shortwave (SW in ) and longwave (LW in ) radiation (Kipp&Zonen CNR4 sensor, 1.0 m above debris surface, stated accuracy ±3%) and relative humidity data (RH: Vaisala HMP45C sensor, 2.15 m above debris surface, stated accuracy ±2%) were recorded at an automatic weather station at the Changri Nup Glacier. Meteorological data from the Changri Nup station were collected at 30 min intervals and resampled to 1 h resolution using an hourly mean algorithm. Precipitation (P) was measured using a Geonor T-200 all-weather rain gauge at the Pheriche site where summer precipitation predominantly occurs as rainfall; these data were corrected for undercatch of solid precipitation using air temperature and wind speed (Sherpa et al., 2017) and the resultant corrected data have an estimated accuracy of ±15%.
Local topography The digital elevation model (DEM) from which slope and aspect were extracted for each sensor site was derived from a series of Surface Extraction from Triangulated Irregular Network Searchspace Minimization (SETSM) DEMs sourced from the Polar Geospatial Centre (University of Minnesota) at 8 m resolution, collected between 8 February and 4 May 2015 (Noh and Howat, 2015). The DEM correction method is detailed in King et al. (2017). Due to the complex and dynamic nature of the glacier surface, topographic parameters at each iButton site were estimated a posteriori from the DEM and are presented here as a generalised local proxies rather than absolute, sitespecific values (Table I). Slope (in degrees) and terrain curvature were extracted for the pixels corresponding to the sensor locations using ESRI's ArcMap v10.1 Spatial Analyst toolbox. Relative terrain roughness was derived using the 'vector ruggedness measurement toolbox', which considers slope and aspect variability for the nine pixels on and around each site location (Sappington et al., 2007). Curvature and roughness metrics both ranged between -1 and +1. In situ observations of the local aspect of each iButton site, measured relative to north, were collected in the field using a magnetic compass with an uncertainty of ± 2°.

Near-surface debris temperature
Daily mean near-surface debris temperature (T s ) for all 16 sites typically exceeded air temperatures (T aP and T aG ) throughout the monsoon period ( Figure 3(a)). Mean T s for the period of observations at the 16 sites ranged from 7.0 to 11.6°C. T s remained close to 0°C between DOY 146 and 152, which was coincident with heavy snowfall in Khumbu valley and the ensuing persistence of a~0.4 m snow layer on the glacier surface. Following DOY 152, the snow cover melted, with the rate and timing of the return to T s >5°C at each site highly varied. Subsequently, from DOY 156 onwards, all T s timeseries exhibited a broadly similar quasi-parallel pattern of change until the end of the observation period. T s appeared to follow a generally rising trend from DOY 156-166, and then a seasonal decrease of approximately -0.1°C d -1 until DOY 210. However, these seasonal rising and falling trends were superimposed with 5 to 8 day cycles in T s , potentially reflecting synoptic variations, and intermittent, shorter (1-3 day) periods with lowered T s . At all 16 sites, T s exhibited marked diurnal variation (Figure 3(b)). Zero amplitudes persisted during the brief period of snow cover (DOY 147-151), the highest daily amplitudes of up to 47°C were found before DOY 170, and progressively declining amplitudes (reducing to a mean of 15°C) characterised the period following DOY 170. Over the monsoon season, the contrasts in T s between the sites were greatest at the start of our observations and between DOY 153 and 170, and declined thereafter, with the least difference between sites seen during the short periods of reduced T s .

Meteorology
Mean daily on-and off-glacier air temperature (T aG and T aP ) followed a similar, but subdued, pattern to the T s data (Figure 3(a)). Air temperature increases of the order of 3°C occurred over the entire study period in both T aP and T aG . The seasonal pattern in T aG and T aP was overlain by a subtle synoptic periodicity with a 5-8 day recurrence. The diurnal amplitudes seen in the T a series were less than those observed for T s . Daily variation in amplitude ranged from 2.1 to 10.4°C for T aP , and from 5.4 to 20.2°C for T aG . In both T a records, diurnal amplitude was greatest during the period of snow cover, and showed a general reduction over the course of the observation period albeit punctuated by short (1-3 day) variability. Off-glacier T aP was consistently lower than on-glacier T aG by a mean difference of 5°C between DOY 145 and 190, and 3°C from DOY 190 onwards. Mean daily SW in displayed an overall seasonal decrease from 405 W m 2 to~217 W m 2 over the observation period, with short-term (<5 days) variability of the order of 200 W m 2 over the study period (Figure 3(c)). Between DOY 148 and 149, SW in was lowest at 123 W m 2 , which corresponded to snowfall and a coincident decrease in T s to 0°C. In contrast, mean daily LW in increased from 253 W m 2 to 320 W m 2 from DOY 143 to 210. Total net incoming radiation (NR in ) was primarily influenced by the pattern of SW in . All three series of radiative energy displayed synoptic (3-8 days) and short-term (1-3 day) variability. Relative humidity displayed a seasonally increasing trend from around 60% on DOY 143 to around 95% by the end of the observation period; this seasonal change was superimposed with shorter-term variability including a brief increase in relative humidity (to >80%) between DOY 146 and 150, aligned with the snowfall and snow cover event (Figure 3(c)). During the snowfall event, total daily precipitation peaked on DOY 150 at 34 mm, but subsequently remained low until DOY 170 and then, as the monsoon progressed further, the magnitude and frequency of precipitation events increased (Figure 3(d)). Increases in total daily precipitation were typically concurrent with decreased SW in and increased LW in and relative humidity.

Timeseries Analyses
A Kolmgorov-Smirnov normality test showed that none of the temperature timeseries (T s or T a ) were normally distributed at 95% confidence level. Therefore, non-parametric analyses were required to interrogate these data further.

Comparison of timeseries
The overall average of mean and standard deviation of T s for all timeseries was 9.2 ±1.3°C, or 9.6 ±1.2°C if the data considered less representative of T s were excluded. Analytical tests indicated that the mean T s timeseries was highly correlated with both T aP (Spearman's r = 0.85, P < 0.05) and T aG (r = 0.78, P < 0.05) but was significantly higher than both the two T a timeseries.
The broad similarity in the individual T s timeseries (Figure 3 (Table II). The generally high correlation (r ≥ 0.88) between timeseries indicated that all sites exhibited a broadly similar general pattern in both periodicity and seasonal trend. However, further comparison using a Kruskal-Wallis test (which tests whether samples originate from the same distribution) showed the T s populations were significantly different (χ 2 = 308.9, or χ 2 = 201.1 excluding the timeseries that were less (c) mean daily incoming shortwave, longwave and total radiation (SW in , LW in and NR in , respectively); (d) total daily precipitation and mean daily relative humidity across the study period.

NEAR SURFACE DEBRIS TEMPERATURE VARIABILITY ON KHUMBU GLACIER, NEPAL
representative of T s , both P << 0.05). To explore the underlying nature and causes for these differences, we: (i) examined the temporal variability in the T s series; (ii) conducted a more detailed assessment of the spatial differences between timeseries; and (iii) explored any associations between T s and the local meteorological and site-specific data. Each of these three sets of analyses are detailed in the following sections.

Temporal variability in near-surface debris temperature
The similarity in the daily T s means and their seasonal pattern, with the exception of the period of snowfall (DOY 146-152), was underlain by a marked reduction in the daily amplitude of variability in T s at all sites over the study period (Figure 3(b)). To test this observation further, regression analysis was employed, with omission of data from the snowfall period. Sites 1, 4, 7, 10, 12 and 16 showed a significant (P < 0.05) decrease in daily mean T s over the observation period, while all other sites showed no such temporal trend (Table III). However, all sites showed a statistically significant increase in daily minimum T s during the monsoon season, averaging 0.08°C d -1 ; and with the exception of Site 13, all sites also showed a significant decrease in daily maximum temperature (mean -0.19°C d -1 ). The concomitant increase in minimum and decreasing maximum T s between timeseries was reinforced by the significant decreasing trend in daily amplitude by a mean of -0.26°C d -1 over the monsoon period at all 16 sites (Table III). These changes were in contrast to air temperature trends, where daily minimum and mean T aG increased by 0.1°C d -1 and 0.04°C d -1 . No significant trend in mean daily maximum T aG, was present, although daily amplitudes decreased by -0.1°C d -1 .
To further examine these seasonal trends in T s amplitude, and to ascertain if there was systematic change in the diurnal pattern of T s fluctuation, we adopted the approach commonly used to analyse synoptic climatology (Brazel et al., 1992;Davis and Kalkstein, 1990), hydrological timeseries (Hannah et al., 2000;Irvine-Fynn et al., 2005;Swift et al., 2005) and ground surface temperature (Lundquist and Cayan, 2007). These previous published analyses used principal components analysis (PCA) to classify patterns of change or modes of variation in diurnally fluctuating timeseries. Here, rather than analyse all 16 T s timeseries, and given the high correlation between all sites (excluding timeseries less representative of T s ) (Table II), a 'representative' timeseries from the data set was used. The most representative T s timeseries was identified using a Nash-Sutcliffe efficiency coefficient (E) typically used to determine the fit of modelled to observed data (Krause et al., 2005;Legates and McCabe, 1999). E was calculated for each T s pair and then summed and averaged for each individual site (Table II). The timeseries with the highest similarity to all other T s series was from Site 14 (∑E = 12.4, mean E = 0.83), and was therefore considered representative.  Table II. A matrix of Spearman rank correlation coefficient (r) and Nash-Sutcliffe efficiency coefficient (E) for each pair of raw (hourly) T s timeseries. All correlations displayed P < 0.05. The greyed rows (Sites 1, 2, 9, 11 and 13) are those identified as being less representative of debris surface temperature due to site clast size. Correlation between each raw T s series and the mean T s is shown, along with the sum and average E for each Spearman's correlation coefficient (r) Debris temperature data from Site 14 were divided into individual diurnal periods of 24 measurements commencing at midnight (00:00). Diurnal periods in which T s was consistently 0°C (DOY 146 to 152) due to lying snow cover were omitted from the analysis. The resultant 61 diurnal data series were reduced and simplified into a number of 'modes' of variation, or principal components (PCs), using PCA without rotation. The first two PCs provided the primary modes of diurnal variation in T s (Figure 5(a)). PC1 accounted for 81.3% of the variance and PC2 for 8.8%. The remaining PCs were discounted as 'noise' because they represented less than 10% of the total variance in the data set (Hannah et al., 2000;Irvine-Fynn et al., 2005). Although absolute loadings were relatively weak (<0.5) for both PCs, a total of 30 days were described best by PC1 and 19 days associated with PC2. A total of 11 days were very weakly related to either PC1 or PC2 (absolute loadings of < 0.09), and were considered to have an undefined diurnal T s cycle ( Figure 5(b), (c)). Of note were the 11 days described by negative loadings on PC2, which contrasted with the consistently positive loadings for PC1, and were suggestive of lagged relationships between the mode of variation and true diurnal T s pattern. Days associated with PC1 predominantly occurred during the former half of the observation period (76% before DOY 176), while 78% of days associated with PC2 and 90% of days with an undefined cycle both occurred following DOY 176 ( Figure 5(c)).
The contrast between the days assigned to the two main PCs and the undefined diurnal cycles were illustrated through a comparison of descriptive statistics (Table IV). The mean diurnal T s was greatest for those days defined by PC1 at 10.9°C, while the mean maximum temperature and diurnal amplitude Not included Figure 5. (a) The two modes of variability in T s for Site 14, described by PC1 and PC2; (b) plot to identify days described by PCs 1 or 2, filled circles identify days with a negative or lagged relationship to PC2 and greyed circles mark days not described by either dominant PC; (c) T s timeseries for Site 14 highlighting each day's mode of variation. Table III. Results of regression analyses to identify seasonal trends in minimum, mean, maximum T s and the associated daily amplitude. Seasonal trend slope (b, in°C d -1 ) is given with the associated P-value, and statistically significant slopes are indicated in italic. was highest compared with days with an undefined T s variation and those associated with PC2. Days that were best described by PC2 exhibited relatively low mean daily amplitude, and mean and maximum diurnal temperatures. The 11 days that were less well defined by PCs had lowest mean, maximum and amplitude in T s . Days described by PC1 were characterised by a lower mean minimum T s (0.9°C) while all other days experienced similar minimum values of T s . The mean time at which T s peaked for each group of days associated with the PCs varied by less than one hour (Table IV).
Subtle variation in diurnal patterns was present in the T s timeseries. There was a clear progressive shift during the monsoon season towards T s exhibiting a lower daily mean, maximum and amplitude, but with a seasonal increase in the minimum T s . The combination of E and PCA analyses explored this further, showing that all sites displayed a regular diurnal pattern of T s during the former part of the monsoon, while there was a systematic shift to more variable and delayed diurnal cycles in the latter half of the observation period. These shifts in magnitude of T s were aligned with the observed seasonal changes in meteorological conditions, specifically with increased precipitation, relative humidity and LW in from around DOY 180.

Spatial variability in debris surface temperature
With evidence of spatial variability between sites most clearly evidenced by the differences in diurnal amplitude between the T s timeseries, further exploration of the spatial contrasts was undertaken. Following the identification of significant difference by a Kruskal-Wallace test, a signed rank pairwise Wilcoxon test provided further detail on spatial variations by comparing pairs of timeseries populations. The representative series from Site 14 was the most similar to all other timeseries, being statistically dissimilar to only Sites 1, 3, 4 and 16 (Table II). Removal of the timeseries considered as less representative of T s made relatively minimal difference to the analysis, suggesting that even the outlying data (Sites 2,9,11,13) were broadly similar to the remaining T s despite the uncertainty arising from varying depth of sensors. A further set of Wilcoxon tests were undertaken on the positively skewed distribution series of maximum, minimum and mean diurnal amplitude of T s . The results of the site comparison data showed 86% and 89% of site pairs had significantly different diurnal amplitudes and maximum T s from one another (P < 0.05), while 39% of the site pairs displayed significantly different minimum T s (P < 0.05).
Daily mean minimum T s for all timeseries varied by -1°C to -4°C between sites, while daily mean maximum T s varied between 10°C and 17°C. While non-parametric correlation coefficients (r) suggested minimal variability between sites with 86% of correlations displaying r ≥0.90 (Table II), such correlations only reveal similarity in timeseries patterns rather than magnitude (Borradaile, 2013). Consequently, notwithstanding the sensitivity of the efficiency criterion (Krause et al., 2005), E was used to compare the strength of each relationship with regards to similarity in both value and pattern for the T s timeseries (Table II). The E values displayed high variability and ranged from -0.42 (Sites 5 and 9) to 0.96 (Sites 7 and 12). The timeseries less representative of T s displayed predominantly lower E values, particularly in their relationships with each other. Spatial variability between the sites appeared relatively small with 84% of E values ≥0.75, suggesting a good similarity in pattern and magnitude between pairs of T s timeseries. For sites located in close proximity to one another (Sites 7-13, omitting those that were less representative of T s ) all the site pairs displayed r ≥0.87 and 80% of these site pairs displayed an E value ≥0.81. However, the contrast in E between timeseries suggests subtle spatial variability in T s did exist between study sites. The correlations between T s remained high (>0.87) even when they were detrended to remove diurnal cycles (following Kristoufek, 2014). This further shows that T s exhibited similar short-term and seasonal variations despite varying sensor locations.
Cross-correlation between the detrended timeseries was used to identify any lag between T s (Table V). Lag times were present for Sites 1 and 2 and a number of other different sites, and with both Sites 8 and 15 for a number of sites. All sites lagged the timeseries at Site 8 by 1 or 2 h, while Site 15 displayed a 1 h lag with 7 sites. Site 8 and 15 were located under 0.010 m and 0.042 m of debris, neither of which are sites where mean clast size, and therefore burial depth, were greatest, and neither sites had been identified as less representative of T s or statistically dissimilar. With regards to the site characteristics, Site 8 was placed in the most northerly aspect and lowest elevation of all iButton locations, while Site 15 had one of the highest elevations and roughness metrics (Table II). Despite a broad statistical similarity in the T s data, there were a number of contrasts in the magnitude, distribution Table V. Correlation coefficient and lag time for pairs of detrended T s timeseries for which the persistent 24-h diurnal cycles have been removed. The grey rows are those identified as being less representative of debris surface temperature due to site clast size Correlation coefficient (r)

Ts1
Ts2 Ts3  Ts4  Ts5  Ts6  Ts7  Ts8  Ts9  Ts10  Ts11  Ts12  Ts13  Ts14  Ts15  Ts16   Ts1   and timing between timeseries. The analysis of the T s data suggested subtle spatial variability in T s was primarily manifested in variability in diurnal T s amplitude, which was principally controlled by variability in maximum T s between sites.

Controls on temporal and spatial variability in near-surface debris temperature
To investigate whether meteorological conditions and site characteristics were associated with controlling T s , and particularly maximum T s , assessment of the influence of meteorological drives and site-specific traits was undertaken using multivariate analysis techniques.
Controls on temporal variability in near-surface debris temperature Controls on temporal variability in T s over the monsoon season were investigated for all hourly timeseries, omitting the period of sustained 0°C in T s in which the debris surface was snow covered. Analysis was undertaken using stepwise multilinear regression (SMR), with meteorological timeseries as predictor variables, to determine the control and combined control of meteorological variables on T s . SMR iteratively adds and removes variables included in the model based on their statistical significance in regression (Draper and Smith, 1998), therefore enabling the relative importance of meteorological variables to be identified. This method is superior to simply regressing individual variables against T s as it can give insight into the extent to which different combinations of meteorological variables control T s . Assessment of the meteorological data demonstrated that none of the timeseries were normally distributed, as for all T s and T a data. Consequently, to transform the T s and meteorological variables to more approximately normal distributions, a simple natural logarithmic conversion was applied. The multivariate models described *T s (where * reflects a log-transform) as a function of *SWin, *LWin, *T aG , *RH (relative humidity) and *P (precipitation). The output from the primary SMR is detailed in Table VI highlighting the relative strength of the relationships between T s and each of the meteorological variables between sites. *T aG was ranked as the most influential predictor of *T s for all sites, with coefficients of determination between R 2 = 0.44 and R 2 = 0.67. The addition of *SW in , *LW in , *RH and *P resulted in only minimal incremental increases in the strength of the correlation between predictor variables and *T s , in all cases resulting in an increase in R 2 of ≤ 0.1. In all cases, *RH was only the third or fourth most significant predictor variable. *P was not significant in terms of contributing to improving prediction of *T s for any site, and was therefore omitted from the model and not included in the first set of results (SMR1) in Table VI. Typically, the sites with the weakest SMR model were those classed as less representative of T s , although Site 16 had similarly low results relative to all sites. One of the potential weaknesses in the first pass SMR models is the collinearity between variables, particularly SW in and T a , for which r = 0.84 (P << 0.05). There is typically a positive relationship between incident solar radiation and T a , due to the direct influence SW in has on T s (Hock, 2003), and the strong covariant relationship present between T s and T a (Foster et al., 2012;Shaw et al., 2016). Consequently, the SMR analyses were re-run with *T aG removed from the model to explore whether additional variables influence T s independent of T aG (Table VI: SMR 2). Results highlighted that, in the absence of T aG , all models exhibited *SW in as the dominant predictor for T s , but with coefficients of determination much reduced (0.17 ≤ R 2 ≤ 0.40). Inclusion of the other meteorological variables, while increasing the models' performance (with R 2 increasing to ≤0.49) Table VI. Results of SMR models describing natural logarithm transformed T s timeseries (*T s ) from meteorological variables and additional predictors derived from the meteorological timeseries (see text for full details).
Predictive variable importance (e.g. 1, 2, etc.) or sequence (e.g. variables 1+2, or all indicated by +3+) is shown, with coefficients of determination and root mean squared error for each model given in parentheses (R 2 , RMSE). The grey rows are those identified as being less representative of debris surface temperature due to site clast size SMR 1: raw transformed meteorological variables  (Table VI). Collinearity between P and RH, or between LW in and RH may also be present but due to the minimal influence of these predictor variables on the SMR results identifying whether such collinearity existed here would be challenging, and so has not been considered further. Conflating the radiation terms (SW in and LW in ) into 'net incident radiation' (NR in ) and continuing the omission of T aG in a third set of SMR analyses (SMR 3) yielded similar results to SMR 2, with *NR in being the dominant predictor variable; moreover, opting for inclusion of 'rate of change in T aG ' (dT a ) for the preceding hour, and cumulative radiation variables (∑SW in and ∑LW in ) and 'time since precipitation' (tP) as potential drivers for T s in SMR 3 showed similarly incremental improvements but only to R 2 = 0.51. In all cases in SMR 3, dT a was the second most significant predictor variable. A final SMR model (SMR 4) excluded all radiation terms and utilised *RH, *P and tP. Despite the close association between incident radiation and T a , the multivariate models using SW in , LW in and NR in were less effective in describing T s change over the monsoon season.
To gain a deeper understanding of the extent to which T s and T aG were related, and whether the two parameters have a varying temporal relationship, T s and T aG were also investigated for daytime (06:00-17:00) and night-time (18:00-05:00) periods separately. A number of previous studies have investigated the seasonal and diurnal variability of T aG (Brock et al., 2010;Steiner and Pellicciotti, 2015) , and in some cases its relationship to T s (Fujita and Sakai, 2000). As elsewhere, days when T s was consistently 0°C (DOY 145-153) were excluded from the correlation analysis. The relationship between T s and T aG varied across the study period for both day and night (Figure 6). The relationship between T s and T aG was predominantly stronger at night (r = 0.86) than in the day (r = 0.75). Daytime T s -T aG correlations varied between r = -0.01 (DOY 190) and r = 0.97, while night-time correlations varied between r = 0.48 (DOY 188) and r = 0.99 (DOY 199). The seasonal and diurnal variation in the relationship between T s and T aG therefore suggests that T aG was the dominant driver of T s but that the strength of this relationship varied across a diurnal period and seasonally, due to diurnal and seasonal variation in additional incident or outgoing energy fluxes that also influence T s . Controls on spatial variability in near-surface debris temperature To determine whether statistically significant relationships between site characteristics and between timeseries existed, as suggested by contrasting diurnal amplitudes and the lags between T s timeseries, a two-step process of analysis was undertaken. Initially, stepwise generalised linear models (SGLMs) were explored to investigate possible controls on variability in T s . SGLMs were undertaken rather than SMR due to the small sample size, and therefore the need to relax the assumptions of normal distribution of each timeseries. The SGLMs examined debris temperature metrics that included means for daily mean T s , maximum T s , minimum T s and the daily mean amplitude of T s for each site as the dependent variables. Site characteristics were used as predictor variables, including elevation, slope, aspect, mean clast size, lithology, terrain curvature and terrain roughness. A simple linear model was used, and potential interactions between site characteristics were not included. The less-representative timeseries (1,2,9,11,13) were omitted from the SGLMs, and 5% significance levels were used to eliminate weaker predictors. Second, following identification of the possible important predictor variables influencing T s identified by the SGLM, linear bivariate regression (LBR) analysis was undertaken between T s variables and the debris variables identified as important in the SGLMs. While the SGLM results give an insight into the combinations of debris characteristics that control the temperature variables, the LBR analysis enable the relationship between the predictor and T s variables to be analysed in isolation.
Results of the SGLMs are given in Table VII, which includes variables that were identified as statistically significant in prediction of T s . None of the models were improved through inclusion of site curvature or roughness, which may be due to the resolution of the DEM causing specific site metrics to be less than exact. The combination of clast size, lithology and slope played significant roles in the SGLMs, with coefficients of determination of around 0.9 for mean, maximum and amplitude T s . Aspect was only considered important for predictions of minimum T s , in which elevation was also critical. The LBR analysis results (Table VIII) show that the relationship between T s variables and debris characteristics identified as influential in the SGLMs were not statistically significant in isolation. The exception was minimum T s and elevation, which had an R 2 of 0.44 (P = 0.02). Consequently, although clast size, lithology and slope influence T s metrics in conjunction with one another, they have little influence on T s independently. Specifically, debris size and lithology are considered to have an impact on the absorption and transfer of solar radiation through a debris layer through their influence on albedo, porosity and moisture content, while slope is a critical factor influencing solar radiation receipt. The southerly aspect of the majority of the sites reported here may undermine identification of the merit in describing T s metrics using aspect. In addition, the lack of prediction of minimum T s by the debris variables except for elevation suggests that minimum T s may be independent of the majority of variables considered, but may be most appropriate for identification using a lapse rate. While the sample set was relatively small, the SGLMs illustrated the potential for physical site characteristics to modulate T s , the importance of considering a suite of debris characteristics and their combined influence in control of T s .

Discussion
The timeseries analyses detailed above identified a number of key aspects in the variability in T s for thick (>1 m) debris on the debris-covered ablation area of Khumbu Glacier. A seasonal trend of decreasing maximum and mean T s was identified at the majority of sites, while an increase in minimum T s was in contrast to seasonal changes in T a . A systematic shift from a dominant smooth diurnal cycle in T s early in the monsoon season to a lagged cycle as the monsoon progressed occurred, alongside which meteorological conditions became more varied. In terms of spatial contrasts, there was evidence of subtle differences between sites, illustrated by disparities in how closely the T s timeseries paralleled each other, and short term (≤2 h) lags in T s between sites. Exploring these differences through consideration of meteorological drivers and potential site characteristic controls enabled identification of a dominant association between T a and T s and the influential role of clast size, lithology and slope on T s metrics at each site. Here, we discuss the processes that may underlie the observed variability in T s on a debris-covered glacier.

Temporal variability in near-surface debris temperature
The near-surface debris temperature (T s ) timeseries were notably perturbed between DOY 145 and 153, during which a period of sustained 0°C occurred following an observed major snowfall event. Following the period of 0°C, short-term variability on timescales of around 3-8 days and a seasonal trend in decreasing maximum T s were observed in all T s timeseries. The timing of short-term variability in T s and SW in , LW in , RH and precipitation was simultaneous, while the seasonal decrease in maximum T s occurred alongside a trend of decreasing SW in , increasing T a , LW in and RH, and increased frequency of precipitation ( Figure 3). The coincidence of the seasonal trends in meteorological variables provide a strong indication of increased cloudiness over the study period (Mölg et al., 2009;Sicart et al., 2006;Van Den Broeke et al., 2006).
Increasing cloud cover results in a decreasing amount of SW in reaching the debris surface, causing maximum T s to decrease, which occurs in all timeseries presented here, and a delay in the time at which maximum T s is achieved as the incoming energy flux to the debris surface is reduced and the debris therefore takes longer to heat up. Consequently, such an increase in cloudiness over the study period would have resulted in the decrease in the diurnal amplitude of T s , and a delay in the timing of peak diurnal T s , both of which are observed in changing modes of variation in T s identified in the PCA (Figure 4). An additional control on decreasing SW in would be that following midsummer (DOY 172) regional SW in and solar angle would decrease, reducing the intensity and duration of SW in a debris surface would receive. However, the decrease in SW in was initiated before DOY 172, suggesting this trend was primarily dependent on increasing cloud cover.
A seasonal increase in cloud cover, relative humidity and the frequency of precipitation would also increase the moisture content of the debris layer. Moisture content of the debris layer has the potential to affect T s considerably (Collier et al., 2014), but is challenging to quantify and not reported here. The presence of moisture in a debris layer affects its effective thermal conductivity and therefore the energy needed to increase bulk temperature. An increased amount of energy would therefore be needed to heat water-filled pores to the same temperature as air-filled pores within the debris layer (Collier et al., 2014;Evatt et al., 2015). Consequently, as incoming energy to the Table VII. Stepwise generalised linear models (SGLMs) for describing debris temperature metrics based on environmental variables for the iButton sensor sites. Models detail the coefficients for each significant (P < 0.05) predictor variable, and summarise the model performance using the coefficient of determination and root mean square error (R 2 , RMSE) debris surface decreased during the monsoon season, and the amount of energy needed to maintain debris layer temperature would increase due to presence of moisture-filled rather than air-filled pores, and mean T s would decrease. In addition, an increasingly moist debris layer would have decreased T s due to enhanced latent heat exchange and subsequent loss of heat through evaporation in the debris surface layer (Cuffey and Paterson, 2010;Takeuchi et al., 2000). These trends in T s are observed in the timeseries presented here, and alongside the precipitation timeseries, suggest debris moisture content may have been a factor in controlling T s . However, direct collection of data for moisture content is needed to confirm the link between T s and debris moisture content.
While the 1-3 day cycles are considered to be the passing of localised weather systems in the Khumbu valley, the 5-8 day cyclic perturbations of T s were synchronous with periods of markedly lower SW in , higher LW in and relative humidity, and higher P. These perturbations suggest the intensity of cloud cover was also temporally variable, resulting in periods of T s with decreased diurnal amplitude and lower maximum T s . The perturbations of T s were increasingly frequent in the latter half of the study period, evidenced by the majority of days loaded to PC2 present in this period. These perturbations suggest that alongside seasonal increase in cloud cover due to progression of the monsoon, more localised weather patterns still contribute to variability in meteorological parameters that also affect T s .

Spatial variability in near-surface debris temperature
Despite the period of asynchronous snow melt and subsequent spatial variation in T s between sites for the period DOY 145-153, all T s data displayed strong similarity for the majority of the study period, evidenced in the r and E values for the raw data and the r values for the detrended timeseries. E values suggested subtle variability did exist between sites, which was primarily manifested in the amplitude and magnitude of temperature recorded at each site rather than the pattern of T s .
Variability in sensor depth may have caused some variability in E between site pairs. Although sensor depth variability was accounted for using the temperature gradient through a debris layer, which was calculated by Nicholson and Benn (2006), their gradients were means of a day (24 h) period. As mentioned previously, applying a daily gradient to determine uncertainty in T s due to depth does not reflect the diurnal variability of temperature with depth, which would affect the magnitude and pattern of T s recorded between sites (Nicholson and Benn 2006). However, after the sites identified as less representative of T s were omitted, sensor depth varied by <0.03 m, which would have produced a maximum uncertainty of 0.44°C between sites (excluding less representative sites) even for the steepest gradients previously identified (at 13:00 by Nicholson and Benn 2006). Variability of T s between sites reached 10°C throughout the study period, which exceeds discrepancies exclusively due to sensor depth and so instead suggests other drivers of spatial variability in T s between sites.
Controls on variability in near-surface debris temperature Coincident trends in T s and meteorological variables suggest a high level of interconnection between meteorological variables and T s . T aG explained the majority of the relationship identified between meteorological variables and T s through SMR for all sites (Petersen et al., 2013), while the other meteorological variables identified to be statistically significant in the SMR1 model (SW in , LW in and RH) were less effective as predictors (Table VI). Omission of T aG in SMR models identified SW in , LW in and RH as contributory drivers of T s , and thus reiterates the complexity of the energy balance at a debris-covered surface where all of the meteorological parameters play some role in controlling T s . However, within the SMR models, the strongest relationship between T aG and T s was R 2 = 0.67, and inclusion of additional variables only improved model performance to a maximum R 2 of 0.68 (Table VI), suggesting T aG is the most important driver of T s , and that temperature-index melt models that calculate T s from T aG will account for at least twothirds of temporal variability in energy input to the debris surface. Consequently, these results suggest that when debris surface temperature is being modelled at a daily time step, a temperature-index model using only the relationship between T aG and T s is considered appropriate, as the incorporation of additional parameters such as SW in and NR in would provide minimal improvement in model performance. However, variation in the relationship between T aG and T s over a diurnal cycle (e.g. Figure 6) suggests the strength of this relationship would be less applicable when modelling T s on sub-daily time steps.
Unravelling the relationship between T aG and T s is complex, as the two variables are interdependent on one another (Shaw et al., 2016), particularly when T a is collected below the standard height of 2 m above the glacier surface in the surface boundary layer (Reid et al., 2012;Wagnon et al., 1999). Critically, here, T aP and T aG were highly correlated (r = 0.72, P < 0.05), but accounting for the elevation difference using a lapse rate of -0.0046°C m -1 appropriate for the monsoon season on Khumbu Glacier (Shea et al., 2015) and a standard lapse rate of -0.0065°C m -1 , exhibited mean residuals between T aP and T aG of -1.9°C and -1.3°C, evidencing the observation that T aG was consistently significantly higher than T aP . This on-/offglacier contrast is due to heat loss from the thick supraglacial debris layer to the near-surface atmosphere through turbulent heat exchange (Takeuchi et al., 2000). Our results mirror those of Steiner and Pellicciotti (2015) where T aP from equivalent elevations was consistently lower than T aG over a debris-covered surface, highlighting the need to use off-glacier temperature records with caution when driving numerical models of glacier ablation, and wherever possible use on-glacier measurements.
The influence of specific meteorological controls of T s was also spatially variable (Table VI). Although a difference in elevation between the T s sensors and the T a sensor existed, variability in the relationship between T aG and T s is predominately attributed to spatial variability between the sites at which T s was recorded. The maximum elevation variation between T s and T aG sensors was 47 m, which, using the range of lapse rates described above, would result in variations in T aG of up to 0.3°C across the study site, which is below the T aG sensor uncertainty. Differences between T a and T s were greater than 0.3°C for all sites. The spatial variability in T s is therefore attributed to variation in a combination of slope, lithology and clast size between sites, variables found to be important for variability in maximum T s between sites, which would result in varying effective thermal conductivity between sites.
The results of the SGLM analysis support previous work on debris-free and debris-covered glaciers, and in permafrost environments, where topographic controls including aspect, slope (Gao et al., 2017;Gubler et al., 2011;Guglielmin et al., 2012;Hock and Holmgren, 1996;Strasser et al., 2004), albedo and surface roughness (considered a factor due to the importance of clast size; Brock et al., 2000;Mölg and Hardy, 2004) were found to influence spatial variability in the incoming energy flux to the ground surface, and would therefore be anticipated to control T s . The most dominant variables describing T s metrics from each site on Khumbu Glacier were slope, clast size and lithology. These variables would be expected to control incident radiation receipt through solar geometry and albedo, moisture content and evaporation, and affect local thermal conductance. However, these debris properties were only found to influence T s metrics in conjunction with one another and were not found to independently control T s . Without further data such as site-specific moisture content and SW in values for each site, the exact controls on such variability cannot be identified. In addition, elevation and aspect were only found to influence minimum T s . The majority of sites reported here were south facing and therefore provide a systematic bias, hindering ultimate identification of the influence of this variable. However, the relatively strong, and statistically significant, relationship between the elevation and minimum T s suggests estimation of minimum T s using lapse rates, and potentially night time temperatures when T s is at its minimum, to estimate spatial variability in T s would be appropriate. The diurnal and seasonal variability in the relationship between T aG and T s identified here builds on the conclusions of Steiner and Pellicciotti (2015), who identified a variation in relationship between the two parameters between night and day and with differing climatic conditions. The occurrence of a seasonal influence in this variable relationship is attributed to variability in meteorological parameters, with decreased strength of the relationship between T aG and T s occurring concurrently with perturbations in SW in , and peaks in LW in and RH (e.g. around DOY 173). Such variability is attributed to differences in the capacity of air and debris to hold thermal energy, and the addition of moisture in either or both environments, causing the relationship to vary between T aG and T s seasonally as well as diurnally. Understanding the importance of the high RH values and precipitation is also important for understanding the effect of turbulent heat flux on glacier ablation for these monsoon-influenced debris-covered glaciers (Suzuki et al., 2007). The correlation coefficients for the T s -T aG relationship presented here also reinforce the findings of Steiner and Pellicciotti (2015), displaying stronger relationships at night due to T s increasing at a greater rate and magnitude than T aG . Consequently, temperature-index melt models with a sub-daily time, which rely on the relationship between T aG and T s , need to consider additional controls on T s such as diurnal and seasonal fluctuations in incoming radiative fluxes, particularly for monsoon-influenced debris-covered glaciers that experience large variability in seasonal weather patterns. However, further investigation into the relationship between meteorological variables, T aG and T s over diurnal cycles is needed to quantify the relative influence of radiative fluxes on T s .

Implications of variability in near-surface debris temperature
While the results of this study provide an interesting insight into the extent of temporal and spatial variability in T s for thick (>1 m) supraglacial debris layers, there is a need to carry out a similar study on thinner debris layers as debris-covered glaciers exist in a range of climatic conditions. Following such studies, a development of surface energy balance models to incorporate spatiotemporal variations in debris properties would be appropriate for modelling mass balance, and also for constraining surface energy balance models used for estimating debris thickness (Foster et al., 2012;Rounce and McKinney, 2014). Our findings advocate the use of a surface energy balance approach for calculating debris layer thickness rather than a direct empirical relationship between T s and debris layer thickness as used by Mihalcea et al. (2008aMihalcea et al. ( , 2008b and Minora et al. (2015). The latter of these approaches oversimplifies the relationship between T s and debris thickness, and omits additional factors such as the influential relationship between SW in and T s , and spatial variability of T s due to varying slope, lithology and clast size of the debris layer. This study suggests that inclusion of site characteristics such as slope and aspect and debris characteristics such as moisture content, porosity, lithology and thermal conductivity would increase the accuracy of results using the surface energy balance approach. Further investigation into the extent of spatial variations in site and debris properties on a glacier scale and the influence of these characteristics on debris surface temperature is therefore needed to constrain such model development. In the meantime, application of either method for estimation of debris thickness (empirical and energy-balanced methods) should acknowledge the possible uncertainty involved in disregarding spatial variability in debris properties and compare their debris thickness estimates with direct field measurements of debris thickness. For energy balance models that calculate glacier mass balance, temporal variation in debris properties, specifically moisture content, have the potential to influence energy exchange at a supraglacial debris surface and through a debris layer, and therefore ablation that occurs under a supraglacial debris layer. However, as highlighted in this study, little is known about such temporal variation in debris properties and so constraining this variability should be the focus of future investigations into supraglacial debris properties.

Conclusions
This study presents the most comprehensive analysis of nearsurface debris temperature (T s ) data for a Himalayan debriscovered glacier to date. The timeseries presented extend beyond describing the influence of debris layer thickness on near-surface debris temperature, and confirm both temporal and spatial variability in T s on Khumbu Glacier. Sixteen sites across Khumbu Glacier's debris-covered ablation area displayed a marked daily cycle in T s , overlying seasonal, short-term and spatial variation in maximum T s and diurnal amplitude. A clear transition in the mode of diurnal variation was associated with increasing cloud cover and precipitation; the latter considered to control debris moisture content. Differences in the magnitude and range of variation in T s were apparent between sites, and were indicative of contrasts in response of T s to meteorological or environmental variables. A close association between on-glacier air temperature (T aG ) and T s was evident while radiative energy had a lesser influence on T s . Analyses of these timeseries also demonstrated the role that the site characteristics slope, lithology and clast size hold in controlling spatial variability in T s when in conjunction with one another, but have little controlling influence on spatial variability of maximum T s in isolation, and that minimum T s is influenced by elevation and aspect. Consequently, this study specifically identified the variables controlling temporal and spatial variability in T s for debris-covered glacier surface with a debris layer thickness of over 1 m.
Our results reinforce the complexity and interconnected nature of the surface energy balance at a supraglacial debris surface, identifying that energy fluxes such as ambient air temperature and incoming radiative flux at the debris surface, as well as debris characteristics such as lithology and clast size to a degree, regulate debris surface temperature but are not independent of one another. Hence, these results suggest that, although temperature-index melt models can be useful for 2711 NEAR SURFACE DEBRIS TEMPERATURE VARIABILITY ON KHUMBU GLACIER, NEPAL estimating supraglacial debris thickness or ablation for daily time steps, these models should follow an enhanced approach in which additional aspects of energy exchange such as incoming solar radiation are included when modelling at a sub-daily (e.g. hourly) resolution (Carenzo et al., 2016). These models also need to consider spatial and temporal variation in the controlling variables used (e.g. air temperature and incoming solar radiation), and use on-glacier air temperature to reduce uncertainties in estimates of ablation. Studies that simulate ablation or derive debris thickness should consider including spatial variability in T s and debris thickness in model calibrations, and consider the influence of variability in site characteristics on these results, in particular with regards to their influence on bulk effective thermal conductivity of the debris layer. Finally, the data presented here were limited to debris layers >1 m thick, and future studies should assess the role of debris characteristics and local topography in defining the energy exchange and T s across thinner debris layers to enable the variability of and controls on surface temperature to be understood across an entire debris-covered glacier surface.