Reduced melt on debris-covered glaciers: investigations from Changri Nup Glacier, Nepal

. Approximately 25 % of the glacierized area in the Everest region is covered by debris, yet the surface mass balance of debris-covered portions of these glaciers has not been measured directly. In this study, ground-based measurements of surface elevation and ice depth are combined with terrestrial photogrammetry, unmanned aerial vehicle (UAV) and satellite elevation models to derive the surface mass balance of the debris-covered tongue of Changri Nup Glacier, located in the Everest region.


Introduction
Predicting the future of the Himalayan cryosphere and water resources depends on understanding the impact of climate change on glaciers (Lutz et al., 2014).About 14-18 % of the total glacierized area in the Himalayas is debris covered (Kääb et al., 2012).This ratio increases to between 25 and 36 % in the Everest region of Nepal (Nuimura et al., 2012;Shea et al., 2015;Thakuri et al., 2014).However, the role played by debris on the surface mass balance of glaciers and, in turn, on the glacier response to climate change remains unclear (Kääb et al., 2012).Indeed, this debris layer insulates the glacier surface from the atmosphere when it reaches a sufficient thickness and complicates the response to climate change compared to clean-ice glaciers (Jouvet et al., 2011;Kirkbride and Deline, 2013;Østrem, 1959;Pellicciotti et al., 2015).
In comparison with debris-free (clean) ice, melt is enhanced when the surface is covered by a very thin layer of debris (1-2 cm) as a result of increased absorption of solar radiation and related heat transfer.On the other hand, debris layers thicker than a few centimetres reduce ice melt rates as less surface heat will be conducted through the debris layer and transferred to the ice (Østrem, 1959;Nakawo and Young, 1981;Mattson, 1993;Kayastha et al., 2000;Mihal-Published by Copernicus Publications on behalf of the European Geosciences Union. cea et al., 2006;Nicholson and Benn, 2006;Reid and Brock, 2010;Lambrecht et al., 2011;Lejeune et al., 2013;Brock et al., 2010).However, several studies based on remote sensing data have shown comparable rates of elevation change on debris-covered and clean-ice glaciers at similar altitudes in the Himalayas and Karakoram (Gardelle et al., 2013;Kääb et al., 2012).Some studies hypothesized that this "debriscover anomaly" could be explained by increased ice cliff ablation and englacial melt on debris-covered glaciers (Buri et al., 2015;Immerzeel et al., 2014;Inoue and Yoshida, 1980;Miles et al., 2016).Yet Ragettli et al. (2015) observed different thinning rates on clean and debris-covered glaciers at similar elevations in Langtang Valley (Nepal) using remote sensing techniques.To date, the debris-cover anomaly hypothesis has not been tested with field-based observation.
To add complexity, the surface area of debris-covered tongues has increased in recent decades due to glacier surface lowering and unstable adjacent slopes, processes that are likely associated with climate change (Bhambri et al., 2011;Bolch et al., 2008;Schmidt and Nüsser, 2009;Shukla et al., 2009).Between 1962 and 2011, the proportion of Everest region glaciers covered by rock debris increased by 17.6 ± 3.1 % (Thakuri et al., 2014) and this proportion could further increase in the future (Rowan et al., 2015).
For these reasons, it is urgent to determine the mass balance sensitivity of debris-covered glaciers to climate change.Unfortunately, there are very few surface mass balance measurements which have been carried out on debris-covered glaciers (Mihalcea et al., 2006).First, the surface mass balance field measurements from ablation stakes are sparse.Second, these measurements cannot be expected to be representative given that the ice ablation exhibits a strong spatial variability depending on the debris thickness or type (Azam et al., 2014;Berthier and Vincent, 2012;Hagg et al., 2008;Inoue and Yoshida, 1980;Mihalcea et al., 2006), and measurements can only be made at locations where the ice surface can be reached.Furthermore, geodetic measurements based on the difference between digital elevation models (DEMs) derived from satellite or aerial imagery only determine surface height change and glacier-wide mass balance and are typically unable to resolve the spatial pattern of surface mass balance (Immerzeel et al., 2014).
In this paper, we assess the surface mass balance of the entire debris-covered tongue of a Himalayan glacier (Changri Nup Glacier) using the ice flux method (Berthier and Vincent, 2012;Nuimura et al., 2011;Nuth et al., 2012).DEMs constructed from (i) terrestrial photogrammetry surveys in 2011 and 2014, (ii) an unmanned aerial vehicle (UAV) survey in 2015 and (iii) satellite stereo pair imagery acquired in 2009 and 2014 are used to estimate changes in glacier thickness.The surface mass balance of the debris-covered area is inferred from the difference between the ice flux measured through a cross section at the upper limit of the debriscovered area and the observed elevation changes.Finally, we compare our field-based estimate of the debris-covered glacier mass balance against surface mass balances observed at nearby debris-free glaciers and quantify the overall reduction in ablation due to debris cover.

Site description
Debris-covered Changri Nup Glacier (27.987 • N, 86.785 • E) is located in the Dudh Koshi catchment, 14 km west of Mt.Everest (Fig. 1).The climate in this region is monsoon dominated and 70-80 % of the annual precipitation falls between June and September (Salerno et al., 2015;Wagnon et al., 2013).Changri Nup is a confined valley glacier with a southeastern aspect and it has no ablation zone tributaries.With a total length of ca. 4 km and a total area of ca.2.7 km 2 it presents a reasonable size for field campaigns.The accumulation zone of the glacier is a cirque surrounded by peaks reaching elevations greater than 6500 m a.s.l.Large serac avalanches feed the accumulation zone from steep southfacing slopes.Most of the ablation zone is covered by debris interspersed with supraglacial ponds and ice cliffs.The debris-covered portion of the tongue has a length of 2.3 km and an average width of 0.7 km, with a terminus located at 5240 m a.s.l.
A smaller debris-free glacier known locally as White Changri Nup Glacier (Fig. 1; 27.97 • N, 86.76 • E) lies 200 m south-west of Changri Nup Glacier.White Changri Nup Glacier has a total area of 0.92 km 2 , a north-eastern aspect and ranges in elevation from 5690 to 5330 m a.s.l.
Additional mass balance data used in this study are taken from nearby Pokalde (27.9 • N, 86.8 • E) and Mera (27.7 • N, 86.9 • E) glaciers, located approximately 7 and 30 km southeast of Changri Nup Glacier respectively (Wagnon et al., 2013).Pokalde Glacier is a small (0.1 km 2 ) north-facing glacier that ranges in elevation from 5690 to 5430 m a.s.l.Mera Glacier is larger (5.1 km 2 ), originates at Mera summit (6420 m a.s.l.) and splits into two distinct branches at 5800 m a.s.l.The Mera branch faces north and west and extends down to 4940 m a.s.l., whereas the Naulek branch faces north-east and terminates at an elevation of 5260 m a.s.l.

Data and methods
A suite of field-based and remote sensing methods were used to calculate the mass balance of clean and debriscovered Changri Nup glaciers.These included photogrammetric surveys, field-based differential global position system (DGPS) and ground-penetrating radar (GPR) surveys, unmanned aerial vehicle (UAV) surveys, point surface mass balance (SMB) measurements and satellite-derived height changes.

Differential GPS measurements
For all ground control points (GCPs) and ablation stake measurements, we used a Topcon DGPS unit.Occupation times were typically 1 min with 1 s sampling, and the number of visible satellites (GPS and GLONAS) was greater than 6.GCP and ablation stake locations have an intrinsic accuracy of ±0.01 m.For surface elevation data collected along the transverse profiles (Sect.3.4), a shorter occupation time (∼ 10 s) was used, and these measurements have an accuracy of ±0.05 m relative to a fixed reference point outside the glacier.Maximum uncertainty for the transverse profile measurements is ±0.10 m for horizontal and vertical components, with the horizontal uncertainty being usually lower.

Photogrammetric surveys
Terrestrial photogrammetric surveys were carried out in the last week of October 2011 and in the last week of November 2014.The photographs were made using a Canon EOS5D Mark II digital reflex camera with a Canon 50 mm f/2.8 AF fixed focus lens.The 21.1 million pixel images were captured in a raw uncompressed format.
Oblique terrestrial photographs that covered most of the debris-covered tongue were collected from three bases and under similar conditions in October 2011 and November 2014.Camera positions were between 1100 m and 2000 m from the glacier, which results in a ground-scaled pixel size of 0.14 to 0.25 m.Camera locations were 280, 264 and 253 m apart and the base formed by the camera locations was roughly perpendicular to the sightings.The base-to-distance ratio (∼ 17 %) enabled a good level of stereovision for manual plotting during restitution.
Photogrammetric measurements were used to build DEMs over the glacier tongue downstream of cross section M (Fig. 2).To geometrically correct the images, GCPs that included 28 large (2 × 2 m) white painted crosses and 12 easily identifiable rocks were measured with a DGPS unit (Sect.3.1; Fig. 2).Baseline lengths between the base and the rover used for GCP measurements did not exceed 2500 m.Additional tie points (36 to 60 points depending on the pair of photographs) on the overlapping images were added to improve consistency.Photogrammetric restitutions were obtained using ArcGIS and ERDAS Stereo Analyst software, and the actual geometric correction was performed with Leica LPS software.Control points used for photogrammetric restitutions had a root mean square error of 0.06 m in XY Z.The accuracy of the photogrammetric DEMs was assessed by comparing spot elevations with DGPS measurements at 25 points not used as GCPs (see Sect. 4.3).The photogrammetric restitution was done manually and elevation contours were constructed at 5 m intervals.

Ground-penetrating radar measurements
GPR measurements were performed on 25 October 2011 to measure ice thickness on the transverse cross section M, located upstream of the debris-covered area at ca. 5525 m a.s.l.(Fig. 2).We used a pulse radar system (Icefield Instruments, Canada) based on the Narod transmitter (Narod and Clarke, 1994) with separate transmitter and receiver, a frequency centred near 4.2 MHz and an antenna length of 10 m.Transmitter and receiver were towed in snow sledges along the transverse profile, separated by a fixed distance of 20 m, and measurements were made every 10 m.The positions of the receiver and the transmitter were recorded with static DGPS measurements and have an accuracy of ±0.1 m (combination of the accuracy of the DGPS and the radar antenna locations).
www.the-cryosphere.net/10/1845/2016/The Cryosphere, 10, 1845Cryosphere, 10, -1858Cryosphere, 10, , 2016 To estimate the ice depth, the speed of electromagnetic wave propagation in ice was assumed to be 167 m µs −1 (Hubbard and Glasser, 2005).Field measurements were performed so that reflections from the glacier bed were more or less vertically aligned with measurement points at the glacier surface, which allowed the glacier bed to be determined in two dimensions.Radargrams were constructed by plotting individual radar traces side-by-side using the DGPSmeasured distance along the transverse profiles.Strong reflectors close to the surface were interpreted as either air within the debris layer or the ice surface.At depth, strong reflections were interpreted as the ice-bedrock interface.The bedrock surface was constructed as an envelope of all ellipse functions, which give all the possible reflection positions between sending and receiving antennas.Using this manual migration, estimates of bedrock depths were then interpolated to reconstruct the glacier-bedrock interface in two dimensions to account for the bed slope.See Azam et al. (2012) for details of the methodology and an example of a radargram acquired on Chhota Shigri Glacier (India) using the same device.

Ice flow velocities and elevation changes from DGPS measurements
DGPS measurements were collected on six transverse profiles located on the tongue of the glacier in the last weeks of October 2011, November 2014 and November 2015 (Fig. 2).
To measure ice flow velocities between 2011 and 2015 along profile M, we made repeat DGPS measurements of (i) ablation stakes and (ii) six painted rocks.Seven bamboo stakes were installed up to a depth of 6 m in 2011 along profile M (Fig. 2).Six of these stakes were replaced in 2014 at their original locations, and all were resurveyed in 2015.To delineate the active part of the glacier from stagnant ice, velocities were also obtained with repeat DGPS measurements performed on more than 75 painted or recognizable rocks in 2011, 2012, 2014 and 2015 (Fig. 2).Some measurements performed on painted stones were discarded when the stones slipped on ice or rolled down on steep slopes.

Unmanned aerial vehicle survey
A detailed survey of the glacier surface was conducted on 22-24 November 2015 using the senseFly eBee UAV.Over the course of five survey flights, a total of 582 photos were collected with the onboard Canon Ixus from an average altitude of 325 m above the glacier surface (Fig. 2).Prior to the survey flights, we collected DGPS measurements of 34 ground control points that consisted of (i) red fabrics with painted white squares and (ii) white crosses used also for the photogrammetry (Fig. 2).Twenty-four GCPs were used to process the imagery and create a DEM with Agisoft, and 10 GCPs were reserved as independent checks on the accuracy of the DEM.The original resolutions of the orthomosaic and DEM are 10 and 20 cm respectively.The images from the survey were processed using the Structure for Motion (SfM) algorithm, which is implemented in the software package Agisoft Photoscan Professional version 1.2.0 (Agisoft, 2014).First, a feature recognition and matching algorithm was applied on a set of overlapping pictures resulting in a set of points in 3-D space derived from the matching features and camera positions.This positioning of the sparse point cloud was then corrected using the DGPS control point measurements.Multiview stereo techniques were then used to generate a dense point cloud of the glacier surface.This dense point cloud was used to construct the DEM, which was then used to generate a geometrically corrected mosaic of all input images.A detailed description of the processing steps can be found in Kraaijenbrink et al. (2016).
Based on the 10 independent GCPs, the average error in the UAV-derived DEM is ±0.04 m in the horizontal and ±0.10 m in the vertical.The removal of an outlier GCP with a vertical error of 0.7 m (the GCP is located on the edge of a large boulder) reduces the average vertical error to ±0.08 m.The resulting orthomosaic and DEM derived from UAV imagery are shown in Fig. 3. Before estimating elevation changes the photogrammetric DEM was interpolated to a regular 5 m resolution raster using a Kriging interpolation method, whereas the UAV DEM was aggregated to the same resolution from its original 20 cm resolution.

Surface height changes from satellite images
To calculate surface height changes over a larger area and longer time period, we also used DEMs derived from two satellite stereo acquisitions.The 2014 DEM was derived from two SPOT7 images acquired on 28 October 2014.The ground resolution of each image is 1.5 m and the base to height ratio between the two images is 0.24.The images are slightly covered by snow above approximately 4800 m a.s.l.The 2014 DEM was derived without GCPs using the commercial software PCI Geomatica 2015.The 2009 DEM was derived from two SPOT5 images acquired on 28 October and 4 November 2009.The ground resolution of each image is 2.5 m and the base to height ratio is 0.45.The 2009 DEM was derived using 23 GCPs extracted from the 2014 SPOT7 DEM and the corresponding 1.5 m orthoimage.Output resolution of both DEMs was set to 6 m.The two DEMs were horizontally shifted to minimize the standard deviation of elevation differences on stable terrain (Berthier et al., 2007).Glaciers were masked out using the inventory from Gardelle et al. (2013).We excluded off-glacier pixels for which the elevation difference was larger than 3 times the normalized median absolute deviation.The vertical shift between the two DEMs was calculated as the median elevation difference on flat and stable zones near the glaciers (1.67 km 2 ).The horizontal shifts were −3.0 and 2.3 m in the easting and northing respectively.The vertical shift was +10.0 m.
The uncertainty of the elevation difference between the two DEMs is assessed from the statistical distribution of the elevation differences over stable terrain (Magnússon et al., 2016;Rolstad et al., 2009).The standard deviation of elevation differences on stable ground (σ STABLE ) is 3.6 m.The decorrelation length estimated from the semivariogram is approximately 50 m, which gives 604 independent pixels for the entire debris-covered tongue (n GLA ), 330 independent pixels for the debris-covered tongue common with the photogrammetric survey (n GLA_COM ), and 668 independent pixels on the stable zone (n STABLE ).Conservatively, we also assumed that the error was 5 times higher in the voids of the DEM (Berthier et al., 2014), which represents n VOIDS /n GLA = 6.6 % of the pixels for the entire tongue and less than 4 % of the pixels for the area in common with the photogrammetric survey.Therefore, we assumed that the total uncertainty for the glacier elevation difference could be obtained as the sum of three independent error sources: the uncertainty on the median elevation difference on stable zones, the standard error on the mean elevation change on glacier and an estimate of the error due to voids in the DEM.By summing these three terms quadratically, we obtain the following: We found σ ELEV = 0.21 m for the total debris-covered tongue and 0.25 m for the area overlapping with the photogrammetric survey.

Point surface mass balance (SMB) measurements
Point SMB measurements, with uncertainties of ±0.20 m w.e. were calculated from annual stake emergences recorded between 2011 and 2015 over both Changri Nup glaciers, as well as Pokalde and Mera glaciers (see Wagnon et al. (2013) for details of the methodology).At Changri Nup Glacier, seven stakes were inserted along the debris-free profile M on 25 October 2011 at approximately 5525 m a.s.l.At White Changri Nup Glacier, eight ablation stakes were inserted to a depth of 10 m on 28-29 October 2010 at elevations ranging from 5390 to 5600 m a.s.l.All stakes on both glaciers were measured annually, so annual local surface mass balance measurements are available since October 2011.On 29 November 2014, a "stake farm" was installed over a 2400 m 2 area in the debris-covered tongue at an elevation of 5470 m a.s.l.(Fig. 2).At the "stake farm", 13 bamboo stakes were inserted to a depth of 4 m, with variable artificial debris thicknesses from 0 (bare ice) to 0.41 m.The debris composition ranged from sand to decimetre-sized gravels, and the stake emergences were measured in November 2015.

Calculation of SMB in the debris-covered area
We estimate the ice flux (m 3 a −1 ) through cross section M using the cross-sectional area obtained from both GPR measurements and surface DGPS survey, and the ice velocities measured at ablation stakes and painted rocks along the profile.Average elevation changes ( h) of the tongue over the periods 2011-2014 and 2011-2015 were obtained by differencing the photogrammetric and UAV DEMs.For the portion of Changri Nup Glacier downstream of the flux gate M, the equation of mass conservation (Berthier and Vincent, 2012;Cuffey and Paterson, 2010;Reynaud et al., 1986) states that the change in surface elevation (h) with time (t) between year 1 (yr1) and year 2 (yr2) is the sum of the area-averaged surface mass balance (B) and the flux term (all terms in m ice a −1 ): where ρ is the density of ice (900 kg m −3 ), FG (m 3 ice a −1 ) is the ice flux through cross section M, front is the flux at the glacier front (equal to zero) and A (m 2 ) is the glacier area downstream of the cross section M. FG A yr1−yr2 is the average emergence velocity below cross section M between www.the-cryosphere.net/10/1845/2016/The Cryosphere, 10, 1845-1858, 2016 year 1 and year 2. Note that the emergence velocity refers to the upward or downward flow of ice relative to the glacier surface (Cuffey and Paterson, 2010).Averaged over the entire ablation zone, it would correspond to the average surface mass balance of this zone for a steady state glacier.Taking into account the elevation changes of the tongue, we can calculate the area average surface mass balance between 2011 and 2014, and 2011 and 2015 in this region.

Ice flow velocities and delineation of the tongue
The demarcation between active glacier flow and stagnant glacier ice downstream of cross section M is crucial for our SMB assessment (Eq.2).However, the strongly heterogeneous debris layer covering this tongue may mask the true glacier margin, and buried ice may not be connected to the active part of the glacier.Indeed, this glacier has been in retreat over the last decades and many stagnant ice areas are no longer connected.From remote sensing optical images, it is very challenging to delineate the margins of debris-covered glaciers (Paul et al., 2013).For instance, several previous studies (Quincey et al., 2009;Rowan et al., 2015) have indicated that Changri Nup Glacier was connected to the Khumbu Glacier, a distance of nearly 3.5 km from the terminus delineated in this study.Similarly, the inventories most commonly used in this region connect the debris-covered Changri Nup, the debris-free Changri Nup and the Changri Shar glaciers Bolch et al., 2011;Gardelle et al., 2013;Nuimura et al., 2012Nuimura et al., , 2015)).
For Changri Nup Glacier, zones of active glacier flow were delineated using horizontal velocities derived from repeat DGPS measurements.Velocities derived from freely available optical imagery (e.g.Landsat) cannot resolve velocities less than 5-10 m a −1 (Paul et al., 2015;Quincey et al., 2009;Rowan et al., 2015), and our measured horizontal ice flow velocities range from 12.7 m a −1 in the vicinity of cross section M to zero at the terminus and margins (Fig. 3).
Despite the presence of stagnant ice far downstream of the terminus, the delineation of the terminus is clear in most places (dashed line in Fig. 3).For example, a proglacial stream flows on a thick layer of sand in a flat area immediately below the delineated snout.However, at some locations the boundary between active and stagnant ice is unclear, and here we spatially interpolated measured ice flow velocities using a Kriging interpolation method (Fig. 3).With this approach and obvious features in the field (slope change, visible ice), the debris-covered ablation area was estimated to be 1.494 km 2 with an uncertainty of 0.16 km 2 , taking into account an uncertainty of ±20 m on the delineated glacier outlines.

Ice flux at the upper cross section of the debris-covered area and tongue-averaged emergence velocity
The ice flux at cross section M (Fig. 2) was obtained by multiplying the surface area of this cross section with the mean cross-sectional ice flow velocity.From the GPR measurements (Fig. 4a), the maximum observed ice thickness is 150 m, and the cross-sectional area was estimated to be 79 300 m 2 in 2011.Taking into account the mean thickness decrease of −0.8 m a −1 at cross section M (Table 1) between 2011 and 2015, we calculated a mean cross-sectional area of 78 200 m 2 for 2015.
A mean cross-sectional ice velocity can be calculated from surface velocities and assumptions about the relation between mean surface velocity and depth-averaged velocity.
Here, two approaches are used to estimate the mean surface velocity.The first uses all surface velocities observed along the flux gate between 2011 and 2015, and the mean The Cryosphere, 10, 1845Cryosphere, 10, -1858Cryosphere, 10, , 2016 www.the-cryosphere.net/10/1845/2016/surface velocity is calculated by fitting a second-order polynomial function (Fig. 4b).Unfortunately, surface velocities were not measured near the glacier margins.We thus assume that ice flow velocity decreases linearly to zero at the margin of the glacier (Fig. 4b), and obtain a mean surface velocity of 9.7 m a −1 from an integral calculation.The second approach infers a mean surface velocity from the centre line surface velocity.The ratio between the mean surface velocity and the centre line surface velocity has been estimated to be between 0.7 and 0.8 for other mountain glaciers (Azam et al., 2012;Berthier and Vincent, 2012).Following this approach, and given that the centre line surface velocity is 12.7 m a −1 , the mean surface velocity is assessed to 9.5 ± 0.6 m a −1 , which is in agreement with the first estimate.The next step is the conversion from mean surface velocity to depth-averaged velocity.Without basal sliding, theoretical calculations suggest that the depth-averaged velocity is 80 % of the mean surface velocity (for n = 3 in Glen's law; Cuffey and Paterson, 2010, p. 310).We do not have any information about the thermal regime of the glacier but we assume that basal sliding is negligible.Our assumption is based on the fact that the glacier is probably cold, as ice in the high-elevation accumulation area (> 6200 m a.s.l.) is transported to lower elevations primarily through serac collapses.Taking the mean surface velocity from the polynomial function (9.7 m a −1 ), we therefore assume that the depth-averaged velocity is 7.8 m a −1 .These assumptions and their influence on the resulting uncertainties are discussed in Sect.4.4.Mean cross-sectional velocity and cross-sectional area are then multiplied to compute an average annual ice flux of 609 960 m 3 a −1 at cross section M over the period 2011-2015.This ice flux, distributed over the mean downstream glacier area of 1.494 km 2 , corresponds to an emergence ice velocity of 0.37 m w.e. a −1 .

Surface elevation changes
Elevation changes are directly measured along DGPS profiles and calculated by differencing DEMs from terrestrial photogrammetry, UAV surveys and satellite stereo pair imagery.

Elevation changes over the area between profiles M and N
For the mostly debris-free region between profiles M and N, where photogrammetric measurements are not available, we calculated a mean elevation change from repeat DGPS measurements along profiles M and N.In general, elevation changes in clean-ice areas are expected to have low crossglacier variability (Berthier and Vincent, 2012;Fischer et al., 2005;Vincent et al., 2009).At profile M, this is confirmed by the similarity in elevation profiles between years, and the mean rate of elevation change is −0.8 m a −1 between 2011 and 2015 at this location (Fig. 5).Along the partly debriscovered profile N, elevation change between 2011 and 2015 is not as homogeneous as for profile M, and the mean rate of elevation change is lower (−0.5 m a −1 between 2011 and 2015; Table 1).As the area is relatively small and has a steady slope, we assume that the elevation change of the region between 2011 and 2015 is equal to the mean rate obtained at profiles M and N, i.e. −0.65 m a −1 .The volume change between profiles M and N is thus 87 343 m 3 over the period 2011-2015.

Elevation changes over the debris-covered area
Downstream of profile N, elevation changes were calculated for two periods (2011-2014 and 2011-2015) by differencing DEMs obtained from terrestrial photogrammetric measurements and the UAV survey.Due to terrain obstruction, thickness changes can only be calculated for 60 % of the ablation area downstream of profile N when photogrammetry data is used.Our results show a highly heterogeneous downwasting pattern of the tongue of Changri Nup Glacier (Figs.

Validation of surface elevation changes
Elevation changes obtained from photogrammetry and UAV data have been validated using DGPS measurements and high-resolution satellite stereo pair imagery.First, we directly compare point elevation data from photogrammetric and UAV DEMs with DGPS elevations observed at independent GCPs (i.e.those not used in DEM generation).Differences between DGPS and photogrammetric elevations for 25 independent GCPs near the terminus and profile R have a root mean square error (RMSE) of 0.63 m.A similar comparison between DGPS spot heights and UAV-derived elevations at 10 independent points gives a RMSE of 0.25 m.A direct comparison between photogrammetric and DGPS elevations on cross-glacier profiles (Fig. 5) shows that differences are generally less than 1 m.We also compare surface elevation changes obtained from photogrammetric DEM differencing and repeat DGPS measurements (Table 1).As photogrammetric measurements are incomplete along the transverse profiles due to terrain obstruction, we only consider the sections of the profiles where both DGPS and photogrammetric elevation data are available.At profiles R, P and Z the differences in rates of surface elevation change measured with the two approaches are approximately 0.1 m a −1 .However, repeat DGPS measurements obtained from transverse profiles are not sufficient to obtain a representative mean elevation change of the tongue despite the numerous profiles.This is a direct result of the high spatial variability of elevation changes in the debriscovered area of the glacier.
From repeat DGPS profiles, we found a mean elevation change (2011-2014) of −0.6 m a −1 (Table 1).In comparison, we observed a mean elevation change of −0.95 m a −1 over the debris-covered area from photogrammetry and UAV data.Consequently, there is no agreement between the mean elevation changes obtained from repeated profiles along the debris-covered tongue and the area-averaged elevation change.Given the large spatial variability in surface height changes over the debris-covered tongue (Fig. 6) and the lack of any clear relation between surface height change and elevation (Table 2; Fig. 6), DGPS profiles cannot be used to assess mean elevation change over debris-covered glaciers.
As a final test, elevation changes downvalley of the delineated glacier terminus were calculated from photogrammetry and UAV data.In this small (0.014 km 2 ) region, average thickness changes of −0.07 and −0.18 m were observed over the periods 2011-2014 and 2011-2015 respectively.These are not significantly different from zero, when the margin of error is considered.However, the unconfirmed presence of stagnant ice in the test area may lead to the slightly negative surface height changes (e.g.Fig. 6c).
Finally, photogrammetric and UAV-derived elevation changes (2011)(2012)(2013)(2014)(2015) can be compared to elevation changes measured from the stereo-pair DEMs (2009)(2010)(2011)(2012)(2013)(2014), though the periods of measurement are slightly different.From 2009 to 2014, we find a mean elevation change of −0.88 m a −1 on the debris-covered tongue downstream of profile M (Table 1; Fig. 6c).If we consider only the areas where the UAVphotogrammetric elevation changes are calculated, the mean elevation change is −0.95 m a −1 (Fig. 6c).This compares well with the mean elevation change of −0.96 m a −1 obtained from photogrammetry and UAV data.Given the uncertainty in the ground-based measurements, elevation changes derived from satellite imagery support the assumption that elevation changes measured on 60 % of the tongue are representative of the whole area.
The Cryosphere, 10, 1845Cryosphere, 10, -1858Cryosphere, 10, , 2016 www.the-cryosphere.net/10/1845/2016/The total uncertainty in our estimated SMB is related to uncertainties in (i) the delineation of the surface area of the tongue, (ii) the elevation changes of the tongue, (iii) the thickness of cross section M and (iv) the mean crosssectional velocity at cross section M. Total uncertainty was assessed following the calculation of the area-averaged surface mass balance (B M ): where B M is the mean SMB (m w.e. a −1 ) downstream of cross section M, ρ is the density of ice, A is the glacier area (m 2 ) downstream of cross section M, h 1 is the elevation change (m a −1 ) between the cross sections M and N, A 1 is the surface area (m 2 ) between cross sections M and N, h 2 is the elevation change (m a −1 ) downstream of cross section N, A 2 is the surface area (m 2 ) downstream the cross section N, S M is the cross-sectional area (m 2 ) at M and U is the mean cross section velocity (m a −1 ) through the flux gate M. Using Eq. ( 3), the overall squared error (σ 2 b ) on the calculated SMB is given by the following: (4) Uncertainties relative to the delineation of the surface areas (σ A1 , σ A2 , and σ A for the surface areas A 1 , A 2 , and A respectively) are assigned a value of ±20 m.The uncertainty relative to the elevation changes ( h 1 ) between profiles M and N is estimated to be ±0.20 m a −1 , based on previous DGPS results.Satellite measurements performed between 2009 and 2014 show that the mean elevation change obtained on 60 % of the surface differs by 0.07 m from the mean elevation change calculated on the whole surface area.Consequently, due to this difference, we assumed an uncertainty σ h2 of 0.1 m relative to the average elevation change h 2 obtained from photogrammetry and UAV data.
The uncertainty relative to the cross-sectional area of profile M has been assessed using an ice thickness uncertainty of 10 m (Bauder et al., 2003).Uncertainty relative to the mean cross-sectional velocity is assumed to be 10 % of the calculated velocity (Huss et al., 2007).Following Eq. ( 4), the overall error σ b on the calculated SMB is thus ±0.2 m w.e. a −1 .

Spatial variability of elevation changes over the debris-covered tongue of Changri Nup Glacier
High-resolution surface elevation changes derived in this study from photogrammetry, UAV surveys and satellite stereo pairs highlight the fact that elevation changes over debris-covered glaciers are highly spatially variable (Fig. 6).This is already well known over debris-covered glaciers where elevation changes depend on both the variability in debris thickness and the spatial distribution of ponds or cliffs (Immerzeel et al., 2014;Nuimura et al., 2012).However, www.the-cryosphere.net/10/1845/2016/The Cryosphere, 10, 1845Cryosphere, 10, -1858Cryosphere, 10, , 2016 this study shows that neither repeat DGPS measurements obtained from transverse profiles nor an ablation stake network are sufficient to obtain a representative mean elevation change or surface mass balance over debris-covered glaciers.
The spatial variability in height changes (Fig. 6) also precludes comparisons between direct (glaciological) observations of SMB on clean and debris-covered glaciers.

5.2
The debris cover controversy: SMBs over debris-covered and clean-ice glaciers in the Khumbu area The overall surface lowering rates and mass balances of debris-covered glaciers remains controversial.Several recent studies suggested that elevation changes on debris-covered and debris-free glaciers are similar in the Himalayas and Karakoram (Gardelle et al., 2013;Kääb et al., 2012;Pellicciotti et al., 2015).Conversely, Nuimura et al. (2012) showed that the debris-covered areas are subject to higher rates of lowering than debris-free areas in Khumbu region, though the 400 m difference in mean elevation between the debriscovered and debris-free areas (5102 and 5521 m a.s.l.respectively) may account for this conclusion.
Comparisons between the mass balances of debris-covered and debris-free glaciers (as opposed to comparisons of surface elevation change only) are hindered by methodological deficiencies and uncertainties.First, geodetic studies typically provide only glacier-or region-wide mass balances based on elevation changes (Bolch et al., 2008(Bolch et al., , 2011;;Nuimura et al., 2012).Geodetic methods are unable to determine a separate surface mass balance for debris-covered areas, because they do not account for the emergence velocity.Moreover, the size, altitude and dynamic behaviour of clean and debris-covered glaciers are different and the comparison between glacier-wide mass balances cannot distinguish ablation rates between debris-covered and debris-free areas.In addition, most of these studies in Nepal have been carried out on catchments with a predominance of debriscovered glaciers (Bolch et al., 2011) and cannot be compared with catchments dominated by debris-free glaciers.Second, the uncertainties related to these remote sensing methods (e.g. the delineation of the glaciers, elevation bias due to the radar penetration into the ice, elevation change assessment and snow density) are large (Pellicciotti et al., 2015).Finally, the regional average mass balances obtained from geodetic methods mask strong differences among glaciers and cannot be used to compare ablation rates between debris-covered and debris-free ice.
In contrast with full-glacier geodetic results, our method based on ice flux calculations and surface lowering observations from photogrammetric and UAV DEMs enables the calculation of an average SMB (−1.21 ± 0.2 m w.e. a −1 ) over the whole debris-covered tongue of Changri Nup Glacier.This assessment includes an area of nearly debris-free ice that represents less than 9 % of the total surface area consid- ered and is thus representative of the debris-covered area for the periods 2009-2014, 2011-2014 and 2011-2015.As our estimate of SMB incorporates the spatial variability in surface lowering, we compare the area-averaged SMB obtained for Changri Nup Glacier with direct SMB measurements from debris-free ice and glaciers in the region (Fig. 7).These include point SMB measurements from profile M (Fig. 2), White Changri Nup Glacier (5390 to 5600 m a.s.l.), Pokalde Glacier (5505 to 5636 m a.s.l.), and Mera and Naulek glaciers (5112 to 5415 m a.s.l.).Also displayed on Fig. 7 are the 2014-2015 point SMB measurements from the stake farm located in the debris-covered area of Changri Nup Glacier (Fig. 2).
The average SMB assessed over the debris-covered Changri Nup Glacier tongue (−1.21 ± 0.2 m w.e. a −1 ) is similar to directly observed SMBs at profile M (−1.50 and −0.85 m w.e. a −1 ) and less negative than measurements from the stake farm (−1.35 to −1.98 m w.e. a −1 ).This implies that (i) the average SMB of the tongue would be much more negative if it was debris free and that (ii) the stake farm measurements are not representative of melt rates over the rest of the debris-covered area.
The Cryosphere, 10, 1845Cryosphere, 10, -1858Cryosphere, 10, , 2016 www.the-cryosphere.net/10/1845/2016/To estimate the effect debris cover has on the SMB, we estimate the average SMB for Changri Nup with the vertical gradient of SMB (−1.4 ± 0.5 m w.e.(100 m) −1 a −1 ; Fig. 7) observed at a nearby debris-free glacier (White Changri Nup).Extrapolating from the mean observed SMB at profile M (−1.16 m w.e. a −1 ), we estimate that an area-averaged SMB of −3.0 m w.e. a −1 would be found for the entire debris-covered area if it were debris-free.The difference between the debris-covered and theoretical debris-free SMB estimates (1.8 ± 0.6 m w.e. a −1 ) represents the overall reduction in melt due to debris cover.
Several studies have suggested that supraglacial ponds and ice cliffs considerably enhance glacier ablation for debriscovered glaciers (Benn et al., 2012;Brun et al., 2016;Buri et al., 2015;Miles et al., 2016;Sakai et al., 2000;Zhang et al., 2011;Juen et al., 2014).Although supraglacial ponds and ice cliffs are present on the debris-covered tongue of the Changri Nup Glacier, the overall mass loss is still considerably reduced due to the debris cover and we conclude that the insulating effect dominates at this site.
This conclusion seems to contradict the results of several studies (Gardelle et al., 2013;Kääb et al., 2012) which revealed comparable rates of elevation changes on debriscovered and clean-ice glaciers.However, these previous results came from geodetic measurements and do not account for the effect of ice dynamics (i.e.difference in emergence velocities between debris-covered and clean-ice glaciers).To overcome this issue, Kääb et al., 2012) compared elevation changes between debris-covered and clean ice using neighbouring ICESat footprints (separated by approximately 1 km), in an attempt to minimize differences in emergence velocity.Still, the geodetic method does not permit direct comparisons of ablation rates, and only the ice flux method employed here allows for reliable estimates of average glacier mass balance over the terminus and comparisons with other glaciers.

Conclusions
The calculated surface mass balance of the debris-covered area of Changri Nup Glacier has been obtained from (i) the ice flux at a cross section close to the boundary between debris-free area and debris-covered area and (ii) elevation changes of the tongue.From the calculated ice flux we estimate an average emergence velocity for the debris-covered tongue of +0.37 m w.e. a −1 .The average surface elevation change between 2011 and 2015, derived from photogrammetric and UAV DEMs, is equal to −0.84 m w.e. a −1 .Consequently, the average emergence velocity does not compensate the surface mass balance, and we infer an average SMB of −1.21 ± 0.20 m w.e. a −1 over the debris-covered area of Changri Nup Glacier (5240-5525 m a.s.l.).
A vertical mass balance gradient derived from nearby debris-free glaciers in the studied region suggests that the av-erage SMB would be −3.0 m w.e. a −1 if the glacier was debris free.This net mass loss reduction of 1.8 ± 0.6 m w.e. a −1 indicates that the surface mass balance is strongly influenced by the debris cover, and we infer that the insulation effect of debris cover largely dominates the enhanced ice ablation due to supraglacial ponds and exposed ice cliffs at this site.
Our method to obtain the surface mass balance of the debris-covered area is reliable.However, the application of the method requires accurate and extensive field data and is hard to transpose to numerous or larger glaciers.A precise delineation of the debris-covered glacier tongue is required.For this purpose, ice flow velocities derived from DGPS field measurements are needed given that ice flow velocities are very low in the debris-covered areas in the vicinity of the margins.In addition, GPR measurements performed on a transverse cross section are also mandatory.
Our results have important implications for studies modelling the future evolution of debris-covered glaciers (Rowan et al., 2015;Shea et al., 2015).An empirical model of debriscovered glacier melt that takes into consideration the relevant processes (surface melt, englacial/subglacial melt and ice cliff and surface pond migration and density) will be an important development.

Data availability
The DEM data can be accessed upon request by contacting Christian Vincent (christian.vincent@univ-grenoblealpes.fr).
and Pakistan.The views expressed are those of the authors and do not necessarily reflect their organizations or funding institutions.

Figure 1 .
Figure 1.Study area overview showing the general location (inset map) and delineation of debris-free and debris-covered Changri Nup glaciers.Background is from ESRI basemap imagery.

Figure 2 .
Figure 2. Map of debris-covered Changri Nup Glacier showing the glacierized area (light blue), DGPS cross sections (blue), delineated debris-covered tongue (dashed black line) and UAV imagery extent (black line).TP is terrestrial photogrammetry and background is from ESRI basemap imagery.TP control markers are painted crosses, and TP control features are characteristic boulders.

Figure 3 .
Figure 3. Map of measured glacier surface velocities (m a −1 ) and location of the glacier margins (dashed line).

Figure 4 .
Figure 4. (a) Cross section of glacier thickness derived from GPR measurements at profile M on 25 October 2011.(b) Measured surface velocities across section M over the period 2011-2015.The dashed line corresponds to a second order polynomial function using all the measurements and forced linearly to zero at the right and left margins.

Figure 5 .
Figure 5. Surface elevation profiles (m a.s.l.) for 2011 (black), 2014 (blue) and 2015 (red) from DGPS measurements (dots), terrestrial photogrammetry (black and blue lines) and UAV survey (red lines).Note that the right (left) bank is on the left (right) of each profile.

Figure 7 .
Figure 7. Surface mass balance as a function of elevation for Changri Nup, Mera and Pokalde glaciers over the period 2011-2015.The upper cross corresponds to the mean surface mass balance obtained on the debris-covered tongue.The blue dashed line represents the mean vertical gradient of mass balance observed at White Changri Nup glaciers and is extrapolated from the mean of SMB measurements at profile M. The lower large cross corresponds to the surface mass balance of a hypothetical clean-ice glacier.Note that surface mass balances of the stake farm on Changri Nup Glacier were measured in 2014-2015 only.The heights of each cross correspond to the uncertainty on inferred SMB.

Table 1 .
Mean elevation changes (m a −1 ) estimated from repeat DGPS measurements and DEM differencing (photogrammetry, UAV and satellite) on cross sections and over the debris-covered tongue (entire and common areas).The letters M to Z refer to cross sections as in