A distributed energy-balance melt model of an alpine debris-covered glacier

Abstract

Distributed energy-balance melt models have rarely been applied to glaciers with extensive supraglacial debris cover. This paper describes the development of a distributed melt model and its application to the debris-covered Miage glacier, western Italian Alps, over two summer seasons. Sub-debris melt rates are calculated using an existing debris energy-balance model (DEB-Model), and melt rates for clean ice, snow and partially debris-covered ice are calculated using standard energy-balance equations. Simulated sub-debris melt rates compare well to ablation stake observations. Melt rates are highest, and most sensitive to air temperature, on areas of dirty, crevassed ice on the middle glacier. Here melt rates are highly spatially variable because the debris thickness and surface type varies markedly. Melt rates are lowest, and least sensitive to air temperature, beneath the thickest debris on the lower glacier. Debris delays and attenuates the melt signal compared to clean ice, with peak melt occurring later in the day with increasing debris thickness. The continuously debris-covered zone consistently provides 30% of total melt throughout the ablation season, with the proportion increasing during cold weather. Sensitivity experiments show that an increase in debris thickness of 0.035 m would offset 18C of atmospheric warming.

Introduction

The ablation rate of ice beneath a layer of debris has a nonlinear dependence on the debris’s thickness, with the relationship determined by the debris’s thermal and radiative properties. At the plot scale, these relationships have been derived for several different glaciers and surface covers (e.g. Østrem, 1959; Mattson and others, 1993; Kirkbride and Dugmore, 2003). At a glacier scale, the significant role played by supraglacial debris in determining the magnitude and rate of ablation, and its influence on volume and thickness changes, has been demonstrated (e.g. Smiraglia and others, 2000; Thomson and others, 2000; Diolaiuti and others 2006, 2009).

Debris-covered glaciers are particularly extensive in high-relief mountain regions such as the Hindu Kush–Himalaya, Caucasus and central Andes (Bown and others, 2008; Lambrecht and others, 2011; Scherler and others, 2011; Bolch and others, 2012), but are also found in many other mountain ranges such as the Western Alps (Deline, 2009), with their melt providing runoff to downstream areas (Xu and others, 2009; Mayer and others, 2010).

Initial distributed melt models developed for debris-covered glaciers applied degree-day methods (e.g. Mihalcea and others, 2006; Singh and others, 2006), or else relied on simplification of the surface energy balance, or measurements of surface temperature, in order to estimate sub-debris ablation rates (e.g. Nakawo and others 1993, 1999; Nicholson and Benn, 2006; Mihalcea and others, 2008a). To date, DEB-Model (Reid and Brock, 2010) and Crocus-DEB (Lejeune and others, 2013) are the only melt models that attempt a process-resolving simulation of energy fluxes at a debris-covered ice surface, forced solely by hourly meteorological observations. These point models also require knowledge of the key debris properties of thickness and conductivity. DEB-Model was incorporated in a distributed melt model for three patches of debris on the mainly debris-free Haut Glacier d’Arolla, Switzerland, by Reid and others (2012). However, no study has yet applied a distributed physically based energy-balance melt model to a glacier with an extensive and variable thickness debris cover over its ablation zone.

Constructing a distributed energy-balance model for a debris-covered glacier is complicated because it requires knowledge of the distribution of the debris thickness across the glacier; the distribution of different surface cover types (e.g. clean ice, dirty ice or debris-covered ice); and the lapse rates of meteorological variables, particularly air temperature, both across debris-covered regions and over the transition between debris-covered and debris-free ice.

We present a distributed energy-balance melt model of Miage glacier, Italy, which, once constructed and evaluated, is used to address the following aims:

1. 1. to quantify the spatial and temporal variations in ablation over different surface cover types (clean ice, snow, dirty ice and debris-covered ice) and their varying contributions to total runoff;

2. 2. to assess the effect of the debris on the amplitude and magnitude of the diurnal ablation signal, which may have important implications for the glacier’s hydrology, dynamics and proglacial runoff hydrograph; and

3. 3. to assess how changes in air temperature and debris thickness affect ablation over the entire glacier.

Study Site

Miage glacier (45º47^ʹ N, 6º52ʹ E) is a 10.5 km² glacier in the Mont Blanc Group of the western Italian Alps, which has a continuous mantle of mainly crystalline supraglacial rock debris (gneisses, mica-schists and granite (Franzetti and others, 2013)) over most of the lower 5 km of the ablation zone (Fig. 1). It has four main tributaries, Mont Blanc, Dome, Bionnassay and Tête Carrée glaciers, which form crevassed icefalls prior to joining the main glacier. Above the confluence of Mont Blanc glacier with the main tongue at ˜2500 m a.s.l., the debris is confined to the central and lateral moraines, with the intervening valleys having only discontinuous and patchy debris (hereafter ‘dirty ice’) ranging in size from silt to boulder. Approximately 42% of the glacier has a continuous debris mantle, with 7% of the glacier surface classified as ‘dirty ice’. Down-glacier of this confluence the mainly continuous debris cover averages 0.20–0.25 m in thickness (Mihalcea and others, 2008b; Foster and others, 2012). The debris originates from rockfalls and mixed snow and rock avalanches from the steep valley sides (Deline, 2009). As the main tongue enters Val Veny it bends 90º east, before splitting into northern and southern lobes and a small central lobe. The glacier’s maximum elevation is ˜4640 m a.s.l., with its lowest elevation ˜1740 m a.s.l.

Fig. 1. Location map of Miage glacier, showing the positions of the meteorological (LWS, UWS and IWS) and gauging stations (GS) and the ablation stakes. Inset gives the position of the glacier in the Alps.

Methods

Meteorological measurements

Data from three weather stations, named LWS (2066 m a.s.l.), UWS (2357 m a.s.l.) and IWS (2411 m a.s.l.) (Fig. 1), were used as model inputs. LWS and UWS are fully automatic weather stations, both measuring incoming shortwave radiation, air temperature, humidity, wind speed and direction, with additional measurements of incoming and outgoing longwave and outgoing shortwave radiation being measured at LWS only. The variables measured by LWS and UWS were sampled at 1 s intervals and recorded as an hourly average. Full details of the instrument specifications at LWS and UWS are given by Brock and others (2010, table 2) and will not be repeated here. Table 1 details additional precipitation and debris humidity instruments installed on LWS and UWS during the 2010 and 2011 ablation seasons. IWS recorded only air temperature during 2011 (Table 1). LWS and UWS were set up on continuous, mainly granitic, debris cover between 4 June and 13 September 2010, and 6 June and 11 October 2011, and between 6 June and 11 September 2010, and 13 June and 13 September 2011, respectively. IWS was positioned over dirty ice to the west of the central moraine on the upper glacier between 27 July and 11 September 2011, to record air temperature representative of debris-free areas of ice and to calculate the air temperature lapse rate from the continuously debris-covered area down-glacier.

Table 1 Details of instruments installed on meteorological stations at Miage glacier in addition to those listed in Brock and others (2010). ‘L’ is LWS, ‘U’ is UWS and ‘I’ is IWS

To estimate evaporation of rainfall from the debris, a volumetric lysimeter was designed following Sakai and others (2004). The volume of water within the lower container was translated into a depth of water (mm) and subtracted from the rainfall measured by the tipping-bucket rain gauge over the same time period at LWS to give evaporation. Condensation was assumed insignificant.

Debris temperature

Debris temperature data were collected at LWS, where the debris was 0.2 m thick, from 3 August to 16 September 2011. Temperatures were monitored at depths of 0, 0.14 and 0.20 m in the debris. Internal debris temperatures were recorded using HOBO smart sensors (Table 1), while surface temperature was calculated from the upwelling longwave radiation recorded by a down-facing Kipp & Zonen CG3 sensor on a CNR1 radiometer at LWS, assuming a debris emissivity of 0.94 (Brock and others, 2010). The sensor at 0.2 m depth measured temperatures above 0ºC on occasion, so likely did not remain in contact with the ice surface at all times. The mean amplitude of the diurnal temperature cycle was 16.78C at the debris surface, 3.4ºC at 0.14 m depth and 1.0ºC at 0.20 m depth (Fig. 2).

Fig. 2. (a–c) Time series of measured and simulated debris temperatures for 0.2 m thick debris near LWS in 2011. (b) shows simulated debris temperatures at 16 cm because they correspond most closely with that measured at 14 cm; it is likely that the temperature probe moved downwards within the debris. (c) shows simulated debris temperatures at 18 cm because they were assumed constant at 0ºC at 20 cm, whereas the probe at 20 cm was likely not entirely in contact with the ice surface. (d) The measured (Meas) and simulated (Sim) average hourly debris temperatures over the same period.

Discharge measurement

The runoff from Miage glacier was used to evaluate simulated melt from the model. The main Miage glacier proglacial stream emanates from the northern lobe (Fig. 1). Melt from the majority of the glacier (except the southern lobe) drains to the northern lobe stream, based on the results from a series of dye-tracing experiments across the glacier in 2010 and 2011 (Fyffe and others, 2012). Simulated melt includes all of the glaciated area, so measured discharge will slightly underestimate the glacier’s total runoff. The level fluctuations in the northern lobe proglacial stream were measured using either a Druck PTX 1830 pressure transducer (all of summer 2010 and June 2011) or an Onset HOBO pressure transducer (remainder of 2011). Both instruments were installed in a stilling well bolted onto a large boulder, upstream of any confluences. The stage record is complete except for 27 and 28 August 2010, 4–8 September 2010 and 18 June to 3 August 2011. The stage record was converted to discharge through a two-part rating curve derived from 16 dilution gaugings. This rating gives a standard error of the estimate of 0.76 m³ s^–1, which gave a percentage error of 14.6% using the average daily discharge in 2010 of 5.37 m³ s^–1 _.

Ablation measurement

The ability of energy-balance models forced by on-glacier meteorological data to replicate measured melt rates on bare ice and snow surfaces has been demonstrated in previous studies (e.g. Reid and others, 2012). Here an extensive set of sub-debris ablation measurements was collected to assess the performance of the model on debris-covered areas. Low-conductivity white plastic ablation stakes were drilled into the ice after removal of the surface debris, using a Kovacs ice drill. The debris was replaced as naturally as possible after stake installation. Fifteen stakes were installed across the glacier in 2010 and six in 2011 (see Fig. 1 for stake locations). Measured ice melt was converted to water equivalent ablation assuming a glacier ice density of 890 kg m^–3.

Model methods

The model was coded using an hourly time step and a 30 m × 30 m grid over the 2010 and 2011 melt seasons (from 8 June to 10 September 2010 and from 14 June to 12 September 2011). It was forced by meteorological data from LWS, UWS and IWS (when available) and base grids of the glacier’s outline, elevation, debris thickness, surface cover (whether clean ice, dirty ice or continuous debris) and the percentage snow cover for each day. Debris surface temperature data are not required to calculate melt since it is derived iteratively within DEB-Model.

Meteorological data

The meteorological data required as model input are hourly incoming shortwave and longwave radiation, air temperature, wind speed, relative humidity and total rainfall, all of which are measured at LWS and extrapolated across the glacier according to the methods described in Table 2. Surface relative humidity was calculated by the Magnus– Tetens approximation using the dew-point temperature from the HOBO sensor and the radiative surface temperature measured at LWS. The HOBO sensor was within a couple of metres of the CNR1 sensor, and on similar debris. Surface relative humidity is only required by the model to indicate whether surface condensation or evaporation can occur, i.e. when the relative humidity is 100%, and the LWS value is applied to all debris cells. Although not included in DEB-Model, it is quite possible that evaporation and condensation occurs within debris (i.e. below the surface) when surface relative humidity is <100%, as the lower fractions of debris were normally observed to be saturated when pits were excavated. However, the magnitude of the internal debris latent heat flux and effect on melt is likely to be small given the close correspondence of internal debris temperatures simulated by DEB-Model and measured values (Fig. 2). In 2010 and part of 2011 the IWS temperature data were not available, so the UPMET to IWS air temperature lapse rate (and hence IWS air temperatures) were estimated using a statistical regression model developed from the 2011 data, using UWS air temperature and relative humidity as inputs (R ² = 0.58, RMSE = 0.0158C m^–1, with the mean lapse rate –0.03968C m^–1 (equal for measured and simulated values)). There were no rainfall data from UWS in 2010, so the average LWS to UWS rainfall lapse rate of 0.0026 mm h^–1 m^–1 from 2011 was used to extrapolate precipitation.

Table 2 Methods used to distribute the meteorological data over the glacier

Meteorological values were distributed based upon the gridcell elevation, or the relation of the cell elevation to a threshold elevation (Table 2). The threshold elevation for incoming shortwave radiation and wind speed was 2218 m a.s.l., because this marks the transition from the upper glacier where the confining glacier trough plays an important role in shading the glacier surface and modifying the surface wind field (Brock and others, 2010). The wind speed above UWS was extrapolated up-glacier using the average lapse rate for each hour between UWS and Helbronner AWS, a meteorological station situated at 3460 m a.s.l. at Punta Helbronner, 8 km east of UWS (data from ARPA della Valle d’Aosta provided by Fondazione Montagne Sicura, Courmayeur). This station is located at the top of the Monte Bianco ski lift (i.e. not on Miage glacier or on a glacier surface), but was at an appropriate altitude to capture the high-elevation wind speeds which are representative of the upper region of the glacier.

The similarity of the sky-view fraction at points above (mean 0.54, range 0.51–0.56) and below (mean 0.7, range 0.62–0.73) the threshold elevation of 2218 m a.s.l. confirmed that, although a large simplification, the total incoming radiation across most of the lower and upper glacier should be reasonably represented by the measured data from the respective meteorological stations (Fig. 3). Following standard geometrical calculations (e.g. Brock and Arnold, 2000), the incoming shortwave radiation was adjusted for each cell’s slope and aspect. This was done using the solar elevation and azimuth angles for each hour and the slope and aspect of each cell, as derived from a high-resolution digital elevation model (DEM; see next subsection).

Fig. 3. Sky-view factors calculated at a range of elevations on Miage glacier, applying the Sky View function of ArcGIS to the DEM.

Rainfall is assumed to fall as snow if the air temperature falls below 08C and does not contribute to modelled discharge. Evaporation of rainfall is zero unless the cell has a ‘debris’ surface type and a snow cover less than 50%, in which case evaporation is 35% of rainfall, based on measured 2010 lysimeter data. The UWS to IWS air temperature lapse rate is only used when the percentage snow cover on the UWS cell is <50%, since when snow overlies debris at UWS the decrease in air temperature between UWS and IWS is due only to the elevation difference and not to the change from a debris to ice surface. In such cases the air temperature is decreased upwards from UWS at 0.00418C m^–1 (Oerlemans, 2010).

Remotely sensed data

A DEM derived from airborne lidar surveys in 2008 was provided by the Regione Autonoma Valle d’Aosta (VDA DEM hereafter). It has a spatial resolution of 2 m, and a vertical accuracy of better than 0.5 m. As for all other grids, the VDA DEM was resampled to a 30 m cell size and projected onto the WGS84 UTM 32N coordinate system (Fig. 4c).

Fig. 4. Grids used as input for the distributed energy-balance melt model: (a) glacier outline grid (a value of 1 signifies a glaciated cell); (b) surface cover grid (a value of 2 is a debris-covered cell, 3 is a clean-ice or snow cell, and 4 is a dirty-ice cell); (c) elevation grid;

(d) debris thickness grid (the maximum debris thickness has been constrained to 1 m in this figure); (e) an example snow grid (8 June 2010); and (f) catchment outline grid (a value of 1 signifies a cell within the catchment).

A 2005 orthorectified aerial photograph of Miage glacier was used to produce the glacier outline grid (Fig. 4a). Snow cones at the margins of the glacier tongue were included even if it was not known whether they have glacier ice beneath them. At higher altitudes, snow was only included if it was attached to the glacier. The catchment outline followed the mountain ridge surrounding the glacier, and the lower glacier outline beyond the ridge. Separate streams drain the area outside the lower glacier lateral moraines, one of which collects water from Marmot lake and Glacier du Breuillat to the north of the glacier’s bend, while another drains Lac du Combal to the southwest of the lower glacier moraines, so these areas do not contribute to measured proglacial stream runoff.

The surface cover types were delimited by visual inspection of the 2005 orthophotograph and defined as either debris-covered ice, clean ice or dirty ice (Fig. 4b). A debris-covered surface includes continuous debris (even if thin) and areas of snow which were otherwise within the debris-covered area. Clean ice includes snow on the upper glacier. Dirty ice encompasses discontinuous debris. The upper limit of dirty ice was difficult to determine, so the elevation limit at which the next elevation zone contained mostly clean ice (2600 m a.s.l.) was used as the upper boundary. Dirty-ice melt was modelled as a distinct surface type due to the observed high melt rates where there is a considerable but incomplete debris cover.

The MODIS10A snow-cover product, which has 500 m tiles, was used to provide daily percentage snow-cover grids. MODIS (Moderate Resolution Imaging Spectroradiometer) snow-cover tiles can be unsuitable because of cloud or erroneous data. Therefore tiles were chosen on days on which snow-cover data were available for the entire study area. Due to the variability in percentage snow cover between adjacent days, an average of two cloud-free tiles from within a few days of each other were used as a reference day (assigned to the middle day between the two measured days), with the percentage snow cover for each day interpolated linearly between reference days.

A debris thickness map for Miage glacier was produced by Foster and others (2012) using Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) AST08 thermal imagery and meteorological data. Evaluation of the debris thickness values produced was hampered by the difference in scale between field measurements of debris thickness (<1 m) and the 90 m × 90 m ASTER pixels used to produce the map. Nevertheless, the range of debris thicknesses measured in the field within corresponding areas was similar and the spatial distribution of debris thickness corresponded well with an empirically derived debris thickness map produced by Mihalcea and others (2008b). These data extend up to ˜2370 m a.s.l. along the centre, and to 2330 and 2350 m a.s.l. on the western and eastern sides of the glacier, respectively (Fig. 4d). In the present paper, debris thickness values above these elevations were modelled using a normal distribution with the mean (0.16 m) and standard deviation (0.097 m) from a transect of debris thickness measurements made at 2334 m a.s.l. in 2006 (Foster and others, 2012).

Energy-balance equations

The energy-balance equations for debris-covered ice were those used by Reid and Brock (2010). The parameters used in their model (DEB-Model) are given in Table 3. DEB-Model calculates sub-debris melt by determining the conductive heat flux at the base of the debris, which depends on the temperature gradient within the debris. The debris cover is split into ten calculation layers, with the conductive heat flux calculated for each, and the heat flux between the lowest debris layer and the ice used to calculate the melt for that debris thickness. The use of several calculation layers allows the assumption of a linear temperature gradient within each layer, but the nonlinearity of the overall temperature gradient is maintained, meaning melt can be modelled at an hourly time step.

Table 3 Parameters used in distributed energy-balance model

For clean ice, dirty ice and snow, the ablation, a (m w.e.), for each model time step is found from the sum of different energy sources at the surface:

(1)

where S ↓ is the incoming shortwave radiation, L is net longwave radiation, H is sensible heat transfer, LE is latent heat transfer, P is the heat flux due to precipitation, ∆t is the model time step, ρ _w is the water density (999.8 kg m^–3 at 0ºC), L _f is the latent heat of fusion of water (3.34 10⁵ J kg^–1) and_i is the clean-ice/snow/dirty-ice albedo. All fluxes are in W m^–2. The conduction of heat or penetration of shortwave radiation into snow or ice was not modelled, and the snow or ice surface was assumed to remain at 08C. It is acknowledged that this is a simplification since the ice or snow surface may fall below freezing at night or at high elevations, but it was deemed outside the scope of the model to simulate subsurface ice temperatures.

DEB-Model does not account for refreezing of meltwater within the debris, which could occur if debris temperatures become lower than freezing, reducing total melt. Below-freezing debris temperatures have been found in winter (Nicholson and Benn, 2012), but the model was run here for a mainly temperate glacier in the summer, when debris temperatures remained primarily above freezing. Temperatures were <0ºC just 0.6% of the time at 0 m, 3.0% of the time at 0.14 m and 0.5% of the time at 0.2 m within 0.2 m thick debris as measured at 10 min intervals by Onset HOBO temperature and relative humidity probes between 3 August and 15 September 2011. Minimum temperatures recorded were –1.51ºC at 0 m, –0.16ºC at 0.14 m and –0.16ºC at 0.2 m.

The different surface types (clean ice, dirty ice and snow) modelled by the clean-ice model were assigned different emissivity, surface roughness length and albedo values (Table 3). On the crevassed regions of the steep tributary glaciers, Dome, Mont Blanc, Tête Carrée and Bionnassay glaciers, the surface roughness length was increased from 0.007 m to 0.05 m (found over ice hummocks 1–1.5 m high on Breiðamerkurjökull, Iceland (Obleitner, 2000)). The crevassed regions were delimited using elevation bands derived from inspection of the 2005 aerial photograph and the VDA DEM. As the surface roughness increases as snow cover reduces and ice surfaces are exposed, the 0.05 m roughness value was only applied when the percentage snow cover was <50%.

Sensible and latent heat fluxes were calculated by the bulk aerodynamic method, using the Richardson number (R _b) to calculate the stability of the surface layer (Brock and others, 2007). To prevent unrealistically large or small R _b values, it was set equal to zero if R _b became greater than 0.2 (turbulence has ceased and airflow is laminar (Andreas, 2002)), or less than –1 (free convection, conditions very unstable (Oke, 1978)). In total, 10% of values were affected by these constraints. Only the negative latent heat flux was calculated and it was constrained over debris to when the surface relative humidity was 100%.

Model structure

For each model time step and gridcell the model assesses first the percentage snow cover, as derived from MODIS data, and then the surface cover type. If the snow cover is 0% then only the model for the surface cover type is run, either for clean ice (clean-ice model with clean-ice parameters), dirty ice (clean-ice model with dirty-ice parameters) or debris (DEB-Model). If the snow cover is 100% then only the snowmelt model is run (clean-ice model with snow parameters). If the snow cover is between these values then both the snowmelt model and the surface cover type model are run, with the melt apportioned accordingly. Snowmelt outside the glacier area is not calculated by the model. Although melt for debris-covered cells will usually be modelled by DEB-Model, where debris is reported to be very thin DEB-Model will predict very high melt rates. This is probably unrealistic since the debris is most likely ‘patchy’ rather than composed of a constant layer a few mm thick. Therefore if the debris thickness of a cell is <0.01 m, the dirty-ice model is run. A total of 89 cells (out of 5042) within the debris surface type had a debris thickness less than 0.01 m. Rainfall is calculated across the entire catchment, with effective rainfall (rainfall minus evaporation) added to melt to give simulated discharge. Total discharge is simply the sum of melt in all cells with no groundwater component and no routing algorithms used.

Results

Meteorological conditions

There were contrasting snow conditions at the start of the 2010 and 2011 ablation seasons. In 2010, following a very cold spring, snow cover was continuous, apart from the crests of moraine ridges, down to 2300 m a.s.l., with significant patches of snow down to 1900 m a.s.l. By contrast, following a very warm spring, the continuous debris-covered zone was largely snow-free by early June 2011. Average conditions over the June–August period were very similar in each year, with 2010 being slightly warmer (mean air temperature 11.1ºC vs 10.5ºC in 2011) based on values for LWS. Mean incoming shortwave radiation was identical in each year (260 W m^–2), and mean wind speed and relative humidity were very similar (2.9 m s^–1 and 65% in 2010, and 3.1 m s^–1 and 62% in 2011). While the month of June was cool in both years (mean temperature 9.4ºC in 2010 and 9.1ºC in 2011) there were contrasts in July and August. Notably, July 2011 was a very cool and cloudy month, with a mean temperature of 9.4ºC compared with a warm 13.1ºC in 2010, while the pattern was reversed in August, with 2010 being cool (mean temperature 10.5ºC) and 2011 warm (mean temperature 12.6ºC).

A full description of the usual meteorological conditions on Miage glacier is provided by Brock and others (2010). In this study the additional measurement of air temperature over dirty ice at IWS is notable. The mean lapse rate with elevation between UWS and IWS was 0.0396ºC m^–1, an order of magnitude larger than a standard clean glacier air temperature lapse rate of 0.0041ºC m^–1 (Oerlemans, 2010), and also much larger than the mean lapse rate between LWS and UWS (0.007ºC m^–1) both of which are situated on debris. This large air temperature decrease is mainly due to a step change in air temperature between the debris-covered and debris-free regions, showing the strong influence of surface conditions on the air temperature regime. Figure 5 shows the air temperature difference is largest during the day, when the debris temperature is greatest, whereas the night-time lapse rate is relatively low and fairly constant. The warm daytime debris surface leads to large upwelling fluxes of sensible heat and longwave radiation (Reid and others, 2012, Fig. 9), which increases the temperature of the overlying air, whereas the opposite tends to happen over bare ice, i.e. cooling of air overlying a cold ice surface.

Fig. 5. Average hourly UWS (2357 m a.s.l.) minus IWS (2411 m a.s.l.) air temperature, measured between 27 July and 11 September 2011.

Model evaluation

Stake measurements

Simulated ablation for each stake was calculated for the cell containing the named stake, and over the same time period (Table 4). For debris-covered ice stakes, the debris thickness in the corresponding model cell was modified to match that of the stake. The daily average ablation for all debris cells was also plotted against debris thickness and compared with season-long measured ablation data from 25 stakes in 2005 (Brock and others, 2010) (Fig. 6). Data from 2005 were used due to the season-long ablation data collected; stakes in 2010 and 2011 generally had shorter record lengths.

Fig. 6. Scatter graphs of average daily simulated ablation against debris thickness, (a) with 2010 simulated data and 2005 measured data, and (b) with 2011 simulated data and 2005 measured data. The x –axis has been limited to focus on measured data; there are therefore some simulated points not graphed.

Table 4 Measured (meas) and simulated (sim) ablation (a) for all stakes in 2010 and 2011

Melt of debris-covered stakes was on average slightly underestimated by the model, with the average relative difference between measured and simulated melt –0.002 m d^–1, or a 9% underestimation. The difference between measured and simulated melt was generally greatest for stakes with short (<1 week) measurement periods (stakes 8, 9, 13 and 14), possibly due to an error in measured ablation. Removing these short-term measurements reduces the average relative difference between measured and simulated melt to –0.001 m d^–1 (5%). This is a reasonable error value given the model assumptions and the fact that stake measurement is only thought accurate to 0.01 m. There was no significant relationship between the model error and debris thickness, and no other relationships were found between model error and the stake elevation or its distance to the relevant meteorological station. This suggests the distribution of meteorological variables was not the main cause of errors between measured and simulated ablation. Although the relationship between debris thickness and ablation is matched well by the model, simulated data have less variability than measured data (Fig. 6), due to differences in debris properties which remain constant in the model (debris thermal conductivity, emissivity, volumetric heat capacity and surface roughness length for momentum) but vary in reality. These debris properties may also change at a sub-cell scale, which can result in variations in ablation between the stake and gridcell scale. Model errors are therefore most likely due to differences between modelled and actual debris properties, suggesting future model development should account for the spatial variation in debris lithology. An underestimation of melt at a certain stake could be due to the thermal conductivity being greater than modelled (0.96 W m^–1 K^–1), due to more conductive rock material or a higher moisture or air content in voids (the conductivity of water is ˜0.58 W m^–1 K^–1 and for air is ˜0.025 W m^–1 K^–1). Measured thermal conductivity at Miage glacier varies from 0.71 W m^–1 K^–1 to 1.37 W m^–1 K^–1 (Brock and others, 2010), which could result in a 26% decrease or 35% increase in melt, respectively, relative to the imposed value of 0.96 W m^–1 K^–1. Non-conductive effects (e.g. evaporation or condensation from within the debris or the percolation of snowmelt through the debris) which were not explicitly modelled may also result in differences between simulated and measured ablation.

Debris temperature data

To confirm that the model sufficiently captures the rate of heat diffusion through the debris, the within-debris temperatures were extracted for the LWS cell with a 0.2 m debris thickness, and compared with measured temperature at the surface and at depths of 14 and 20 cm between 3 August and 15 September 2011 (Fig. 2). Using the average hourly temperature cycle, the lag time of the temperature peak from the surface to 14 cm was 2 hours 30 min (velocity of 1.6 × 10^–5 m s^–1) and from the surface to 20 cm was 2 hours 40 min (velocity of 2.1 × 10^–5 m s^–1). Simulated and measured hourly surface temperatures match well (RMSE = 1.8ºC, R ² = 0.95), although the 14 cm probe was matched most closely by the simulated temperatures at 16 cm (RMSE = 1.7ºC, R ² = 0.95). This was possibly due to the probe moving down slightly within the debris after installation, and is within the expected error when positioning such a probe within debris. The 20 cm probe matched temperatures between those simulated at 18 and 20 cm (assumed 0ºC in the model); this is reasonable given the difficulty of ensuring the probe remained in contact with the ice and the tendency for probes to migrate within the debris matrix over time. At the surface, simulated and measured average peak temperature and the timing of the peak were the same (to 0.18C), and the 16 cm (simulated) and 14 cm (measured) data had very similar average peak temperatures (4.58C and 4.28C, respectively), although mean simulated 18 cm peak temperatures were 1.08C higher than those of the measured 20 cm probe. Therefore, with the data available, the model replicates the temperature cycle through the debris well.

Discharge data

Measured and simulated daily proglacial discharge is given in Figure 7. Mean proglacial discharge during the measurement period was 5.25 m³ s^–1 in 2010, and 6.59 m³ s^–1 in 2011. In 2010 the RMSE between measured and simulated discharge was 3.26 m³ s^–1, and in 2011 was 3.77 m³ s^–1, with the model bias –0.94 in 2010 and –2.60 in 2011. However, as the model has no routing routines to account for the timing of different components of runoff, the comparison is only given to verify that the total melt produced over the glacier was realistic, even if there are differences in timing and magnitude of peak runoff between the data and model. On 12 and 16 June 2010, simulated discharges were much greater than measured, with peaks composed mainly of rainfall. This was not replicated later in both seasons when peaks in simulated runoff were usually less than measured (except on 16 July 2010). This suggests a proportion of rainfall in June refroze within the snowpack, which was deep and continuous above 2300 m a.s.l. at the beginning of the 2010 field season. Refreezing occurs when melt (and presumably rainfall) percolates through the cold firn in springtime (Reijmer and Hock, 2008), especially when liquid water contacts freezing glacier ice at the base of the snowpack (Bøggild, 2007). Both refreezing and superimposed ice formation reduce immediate runoff. Water could also be stored beneath the glacier, either in the distributed subglacial system or in a till aquifer beneath the glacier (Jansson and others, 2003). Later in the season, underestimation of measured discharge is more likely, due to rain-gauge undercatch (Marsh and Dixon, 2012), which ranges from 10% to 90% for selected catchments in the Swiss Alps (Farinotti and others, 2012). Rain fell more frequently at the UWS compared to LWS gauge, so rainfall extrapolated from LWS may underestimate the rainfall frequency. Furthermore, both gauges are at a relatively low elevation, so rainfall amounts over extensive areas of the catchment above 2357 m a.s.l. are unknown. In spite of these potential errors the magnitude of simulated discharge was similar to that measured overall.

Fig. 7. Measured and simulated runoff, and simulated melt (excluding rainfall), together with measured rainfall at LWS in the lower plot, for (a) 2010 and (b) 2011. The y –axis in (b) has been constrained to 20 m³ s^–1. There are no discharge data for 27 and 28 August 2010, 4–8 September 2010 or 18 June to 3 August 2011, due to either logger issues or damage to the gauging station.

Spatial variation in ablation

The average daily simulated ablation is shown in Figure 8, with the relationship between ablation and elevation given in Figure 9. Average daily ablation of all modelled cells mentioned includes the influence of snow cover.

Fig. 8. Average daily simulated ablation (m w.e. d^–1) for (a) 2010 and (b) 2011. Ablation includes snowmelt for each cell.

Fig. 9. Average daily simulated ablation against elevation for each glacier cell in the distributed model, split by surface cover type for (a) 2010 and (b) 2011. The melt shown includes snowmelt.

Ablation is small at the top of the accumulation zone, but all cells experience some melt over the season (up to 4673 m a.s.l.). Figure 8 shows that the rate of clean-ice ablation on the tributary glaciers increases down-glacier, likely due to higher availability of energy, as well as a decrease in snow cover. Ablation rates at 3500 m a.s.l. were ˜0.02–0.03 m w.e. d^–1, but a short distance down-glacier, on the icefalls of Dome and Mont Blanc glaciers, melt rates are higher (˜0.05 m w.e. d^–1), due to a favourable south-westerly aspect for shortwave radiation receipts and high surface roughness. Parts of Bionnassay and Tête Carrée glaciers also had regions with a higher surface roughness, but as their elevation (Bionnassay) and aspect (Tête Carrée) were not as favourable to ablation, the enhancement of ablation was less marked. On the lower parts of Dome and Mont Blanc glaciers, where the ice becomes dirty, this decrease in albedo combines with the high roughness and favourable aspect to give the region the highest ablation rates on the glacier (0.05–0.06 m w.e. d^–1 on average). These simulated values compare well with mean measured ablation rates at three stakes installed on the same area of dirty ice, close to IWS, throughout the 2005 and 2007 seasons (Brock and others, 2010). High melt rates (identified as regions of thinning) in the sparsely debris-covered regions just up-glacier of the zone of continuous debris cover have been observed on glaciers such as Khumbu Glacier, Nepal (Nakawo and others, 1999), six glaciers in the Adylsu valley, Caucasus (Stokes and others, 2007), and on Tasman Glacier, New Zealand, where annual ablation was 17.2 m w.e. 10 km above the terminus, but <2 m w.e. close to the terminus (Kirkbride and Warren, 1999).

The mid-part of the glacier has a large range of ablation rates which are determined more by the surface cover type and debris thickness than elevation or aspect (Fig. 9). Debris-covered regions following the lateral and medial moraines have low ablation rates of 0.02 m w.e d^–1, but individual cells with very thin debris have high ablation rates (e.g. one cell with a debris thickness of 0.014 m had an average ablation of 0.05 m w.e. d^–1). Similarly, the dirty-ice areas have ablation rates around 0.05 m w.e. d^–1, although this hides the great variation in small-scale ablation rates which were not modelled due to the complexity of the processes involved (e.g. boulder table and cryoconite formation).

Down-glacier of 2300 m a.s.l. the glacier has an almost continuous debris cover, and ablation is determined by debris thickness variations. The increase in debris thickness with distance down-glacier leads to decreasing ablation near the terminus; on the lowest part of the northern lobe, average daily melt is only 0.002 m w.e. d^–1. Figure 9 shows that the along-glacier variation in ablation rate is much more complex than is assumed by a simple ‘reverse ablation gradient’ for a debris-covered glacier, with spatial variability in ablation rates increasing markedly in the middle section of the glacier. This high variability results from contrasting surface types, including areas of dirty and clean ice and areas of debris of varying thickness.

Table 5 shows that the lowest individual ablation rates were seen on clean-ice cells (at high altitude). On average, debris-covered cells have the lowest ablation, with dirty ice consistently highest.

Table 5 Minimum, mean and maximum simulated values of average daily ablation (m w.e. d^–1) for each surface cover type in 2010 and 2011. These values include the influence of snow cover for each cell

Timing of ablation

The time of maximum melt occurs when the temperature gradient within the lowest layer of debris maximizes. To investigate how debris thickness affects the timing of the peak melt rate across the debris-covered zone, the average hourly melt cycle was calculated for all debris-covered cells between 8 June and 14 July 2010. The time of peak melt for each cell is plotted against debris thickness (Fig. 10), showing the timing of peak melt is delayed by the debris, by an amount that increases with the debris thickness. Thinner debris (<0.5 m) shows greater lag-time variability, due to variations in the meteorological conditions across the glacier. Very thick debris (>0.85 m) gives a lag time of 24 hours, because there is no longer a daily cycle of temperature at this depth but instead temperature increases continuously, implying the melt signal is so damped it reflects only the gradually increasing warmth as the season progresses, i.e. the seasonal cycle rather than the diurnal cycle is dominant. For the average debris thickness on the glacier (0.25 m; Foster and others, 2012) the average time of peak melt is ˜16:30.

Fig. 10. Time of peak melt against debris thickness, using the average hourly cycle of all cells modelled as debris-covered for the first 39 days of the 2010 data.

If one compares the time of maximum total melt from each surface type (Fig. 11), peak melt of snow and clean ice occurs at 13:00, of dirty ice at 14:00 and of debris-covered ice at 15:00. The debris melt peak is skewed to thin debris which has a short lag time and produces the bulk of the runoff from this surface type.

Fig. 11. Average total hourly simulated ablation for each surface cover type in 2010.

As well as inducing a lag to peak melt, the debris attenuates sub-daily variations in melt compared to clean-ice melt. Furthermore, ablation beneath debris continues longer into the evening due to the lag caused by the conduction of heat to the base of the debris. These effects are more pronounced beneath thicker debris.

These delays have implications for the rest of the hydro-logical system, because they reduce the diurnal amplitude of supraglacial stream discharges and result in a later input discharge peak. This could be translated into a later and less pronounced daily peak in the proglacial runoff record.

Seasonal fluctuations of melt sources

At the beginning of June 2010, when the glacier was snow-covered down to 2300 m a.s.l., total melt was composed mainly of snowmelt (>80%), with the rest sub-debris melt, and a very small quantity from dirty-ice melt (Fig. 12). The proportion of snowmelt decreased over the season, being replaced firstly with sub-debris melt (up to 30%), and from 9 July 2010 onwards with dirty-ice and clean-ice melt. At the beginning of the 2011 season the snow cover was less extensive than in 2010, with the early June contribution of snowmelt being just 44%, the majority of the remainder of melt being composed of sub-debris melt (30%) and the rest split between clean (14%) and dirty ice (12%) (Fig. 12). There was still a decrease in the proportion of snowmelt due to the exposure of clean ice. In both years there were specific days when the proportion of snowmelt was decreased and replaced by sub-debris, clean-ice and dirty-ice melt (e.g. on 20 June, 30 July and 14 August 2010). On these days, low air temperatures reduced the elevation of the 08C isotherm, decreasing the contributing area of snowmelt and increasing the proportion of clean, dirty and sub-debris melt to runoff. Figure 7 shows these days were associated with relative runoff lows. There was a particularly strong anomaly on 8 August 2011, when sub-debris melt constituted 62% of total melt, with clean ice 24%, dirty ice 13% and snowmelt only 1%, caused by low air temperatures.

Fig. 12. Average daily variations in the proportions of total melt, of debris-covered ice, dirty ice, snow and clean ice: (a) 2010 and (b) 2011.

Clean and dirty ice and snow contribute a greater proportion of runoff than their area, while the debris-covered region contributes a smaller proportion of runoff than its area due to relatively low melt rates. The contribution of debris-covered melt is fairly consistent at 30% both within and between seasons.

Sensitivity analysis

Although the methods used to distribute the meteorological variables across the glacier were carefully considered, and the parameters used were as far as possible from measured data, it is nevertheless important to consider the sensitivity of the model results to these methods and values. This subsection also considers the sensitivity of the model results to air temperature and debris thickness, in an attempt to understand the likely impact of future changes on the quantity of melt produced.

Sensitivity to meteorological lapse rates

In the model, the downwelling longwave radiation is distributed across the glacier from the measured value at LWS using the lapse rate of 0.031 W m^–2 m^–1 given in Marty and others (2002). However, this lapse rate was not measured on the glacier, and the suitability of its application to debris-covered glaciers is unknown. The gradient of the down-welling longwave lapse rate was therefore varied as a test by ±10%. A 10% increase in the downwelling longwave lapse rate resulted in a 1.2% decrease in melt and vice versa, showing that simulated melt was not particularly sensitive to the value chosen. There is uncertainty in the use of the Helbronner data to derive a wind speed lapse rate that applies above UWS, since Helbronner station was not situated on the glacier. The model was run alternately with the wind speed at LWS applied at all cells with an elevation less than 2218 m a.s.l., and the wind speed at UWS applied at all cells with an elevation above this threshold. This resulted in a small (1.6%) decrease in melt, because the lower wind speeds at higher elevations reduced the sensible heat flux. It is likely, however, that the wind speeds do increase at higher elevations since the sheltering effect of the glacial trough will decrease with distance up-glacier. It is known that there is a significant difference in air temperature between the debris-covered and debris-free areas (Fig. 5); however, there is uncertainty regarding how the temperature varies over the 1.5 km region of different surface covers between UWS and IWS. In the base model, the UWS to IWS lapse rate is applied to this region, but it may be that surface cover type is more important than cell elevation in determining air temperature. An alternative method was applied so that on debris-covered cells the LWS to UWS lapse rate was applied from the UWS air temperature, and on clean- and dirty-ice areas the standard clean glacier lapse rate (–0.00418C m^–1; Oerle-mans, 2010) was applied from the IWS air temperature. This resulted in a very small increase in melt (0.01%).

Sensitivity to parameters

The influence of varying the model parameters on the average daily melt is shown in Table 6. This shows that simulated melt is insensitive to changes in the aerodynamic roughness length, probably because the turbulent fluxes tend to occur in opposite directions over debris when compared to ice and snow (e.g. the sensible heat flux would tend to be negative over debris but positive over ice during the day), so the effect on melt averaged over the glacier cancels out. Varying the debris thermal conductivity also only caused a small variation in melt, but this small percentage change will be partly due to sub-debris melt only accounting for ˜30% of the total (Table 7). Simulated melt is moderately sensitive to albedo, and since albedo is highly variable both spatially and temporally (especially over ice and snow), incorporating a more detailed parameterization of albedo in a future version of the model would be useful. Although the debris albedo is changed by the debris lithology, the range is relatively small (from 0.12 for areas of schist to 0.16 for areas of granite and quartz-rich rocks (Brock and others, 2010)) and the debris is generally composed of a mixture of rock types, although the albedo does decrease when the debris is wet. Changing the emissivity ±0.02 (equivalent to a change of ˜2%) resulted in a small change in total melt (2.8%). Therefore, although the model is sensitive to the value of emissivity chosen, the small natural range of emissivity values means the effect on total melt will likely be small.

Table 6 Effect of variations in model parameters on total melt. The values for all surface types were applied within each model run

Table 7 Proportion of melt from different surface cover types in 2010 and 2011

Sensitivity to air temperature and debris thickness

The melt model presented in this paper could be used to estimate changes in glacial ablation caused by future climate change. A first step towards this is to identify the sensitivity of simulated melt to air temperature and debris thickness variations. The predicted air temperature change for 2080–99 for the southern Europe and Mediterranean region is in the range 2.2–5.1ºC (mean 3.5ºC), with the increase larger in June–August (2.7–6.5ºC), according to the A1B emissions scenario (Solomon and others, 2007). An increase in the thickness and areal coverage of debris has been reported for several debris-covered glaciers (e.g. Stokes and others, 2007; Kellerer-Pirklbauer and others, 2008), but future rates of debris thickness change are yet to be identified. Air temperature was varied by changing the hourly record of the LWS and UWS air temperature, and the debris thickness was varied by changing the debris thickness of all cells equally. The model was run for 2010 when the IWS air temperature was modelled from UWS meteorological variables. Cells with a debris thickness less than 0.01 m were modelled with the dirty-ice model, even if they were within the debris-covered surface area: thus, increasing the debris thickness by 0.01 m will have the effect of switching some dirty-ice cells to debris-covered, essentially increasing the debris-covered area (and vice versa for a decrease in debris thickness). Apart from this, the boundaries between different surface cover types remained constant. The combined influence of air temperature and debris thickness was also investigated. Of course, since higher air temperatures increase melt, this would in turn increase the debris thickness, at a rate determined by the quantity of melt and englacial debris concentration (Glazyrin, 1975; Bozhinskiy and others, 1986). This influence was not investigated here, but it could be tested if the englacial debris concentration were known. Variations in debris thickness will become more important over longer (decadal) timescales, and should be included if changes in melt rates in the future are to be assessed.

A 1ºC increase in air temperature resulted in a 6.25% increase in melt (Table 8). However, the relationship between the air temperature change and melt is not linear but very slightly exponential. This was not expected since sub-debris melt is generally less sensitive to air temperature variations than clean ice (Brock and others, 2010), but was mostly caused by the nonlinear response of the dirty-ice cells to increased air temperature. Increasing the air temperature resulted in the smallest increases in ablation where the debris was thickest (the melt increase was <0.05 m over the season), and at very high altitudes. The largest increase in ablation was found in areas with a high surface roughness and favourable aspect, i.e. Dome and Mont Blanc glaciers. The areas of clean ice, even at high elevations, had a larger increase in melt than the majority of the debris-covered tongue. When snow overlay debris in the spring, it counteracted the insensitivity of the sub-debris melt, meaning increases in ablation caused by higher temperatures will be greater in spring. This suggests an increase in spring air temperatures may result in a larger ablation increase than an equivalent increase in summer temperatures.

Table 8 Sensitivity of simulated average daily glacier melt (m³ s^–1) for the 2010 season to air temperature and debris thickness variations, with the % difference given in parentheses

Increasing the debris thickness by 0.01 m resulted in a 1.9% decrease in melt (Table 8). The change in melt with debris thickness was not linear (e.g. the decrease in melt resulting from a 0.01 m debris thickness increase was greater than that for an increase in debris thickness from 0.01 m to 0.02 m). The melt response to a change in debris thickness of ±0.01 m is affected by the switching of some cells from dirty-ice to debris-covered and vice versa. An increase in debris thickness decreased melt the most where the debris was thinner, because of the more sensitive relationship between debris thickness and melt at low thicknesses. If a dirty-ice cell was switched to debris-covered this would cause a large increase in melt, since very thin covers give high melt rates in DEB-Model. However, if a cell had a debris thickness just greater than 0.01 m originally then increasing the thickness resulted in a large decrease in melt. This does not exactly replicate the relationship between melt and debris thickness described by the Østrem curve, but does show how the effect of a debris thickness change on ablation is determined by the original thickness and its relation to the effective thickness (Kirkbride and Dugmore, 2003).

When the +1ºC and 0.01 m scenarios were combined, for regions of thick debris, the influence of the increased temperature overcame that of an increased debris thickness, to produce a melt increase amounting to <0.04 m w.e. over the season. Conversely, for thinly debris-covered regions, the decreased melt due to thicker debris dominated over the effect of the increase in temperature. On cells where the increase in debris thickness increased ablation, melt was further increased by higher air temperature. After the addition of the first 0.01 m of debris, all cells within the debris-covered area are modelled with DEB-Model, so the change in ablation given a 0.01 m increase but no change in cell type gives a 1.76% decrease in melt. The debris thickness therefore needs to increase by ˜0.035 m on average to counteract a 18C rise in air temperature.

The finding that melt rates are highest on the lower parts of the southwesterly-facing tributary glaciers corresponds with the recent evolution of Miage glacier, with both tributaries thinning markedly and Mont Blanc glacier becoming detached from the main glacier tongue. This mirrors Diolaiuti and others’ (2009) orthophotographic analysis of Miage glacier that revealed recent mass loss was greatest in areas of thinner debris.

Discussion and Conclusions

The model described in this paper has provided a fuller understanding of the spatial and temporal variations in ablation on an alpine debris-covered glacier. Although the models used to calculate melt over the debris-covered and debris-free regions had already been developed, to produce a distributed melt model it was necessary to delimit the different surface types, calculate snow cover from satellite data and extrapolate meteorological variables. The latter presented the biggest challenge, with three stations required to account for variable lapse rates on, and between, debris-covered and debris-free areas. Although resulting in an inevitable simplification of the distribution of meteorological variables, the availability of these data meant that lapse rates were developed using on-glacier data. Despite the lack of calibration of parameters, all of which were obtained from measured or published data, the model replicated measured melt reasonably well, and can go some way to modelling proglacial discharge, although uncertainties in catchment rainfall amounts and the lack of water routing hinder a direct comparison.

The highest rates and variability of melt occurred on the mid-part of the glacier, which is characterized by a mixture of surface types including dirty ice and areas of both scattered and thick debris cover. While melt rates decrease towards both the terminus (due to increasing debris thickness) and the upper glacier limit (due to increasing elevation), the along-glacier variation in melt rate on a debris-covered glacier is more complex than the sometimes-assumed ‘reverse ablation gradient’. Melt rates were also high on the crevassed regions of Dome and Mont Blanc glaciers situated above the main glacier tongue, due to high surface roughness and their southwesterly aspect which enhanced the net shortwave radiation flux. The debris-covered zone contributed 27–30% of total runoff, a relatively small amount considering it occupies 42% of the total area of Miage glacier and is concentrated at the lowest (i.e. warmest) elevations. The contribution of sub-debris melt to total runoff remains fairly constant over the ablation season, in contrast to bare ice and snow whose contributions over the season increase and decrease, respectively.

The amplitude of the daily cycle of ablation is attenuated with increasing debris thickness, and for debris >0.85 m thick there is no longer a daily cycle; instead the melt rate responds to long-term (seasonal) temperature cycles. This agrees with the measured attenuation of the temperature signal through debris (Nicholson and Benn, 2012), and contrasts with the more rapid melt response of bare snow and ice surfaces to energy input at the glacier surface. The attenuation of the melt signal means that although ablation occurs into the evening, the quantity of peak daily melt is reduced. This has consequences for supraglacial meltwater production, which could influence the structure and development of the englacial and subglacial hydrological system and possibly glacier dynamics (Fyffe and others, 2012).

Model sensitivity experiments revealed that using alternative strategies to distribute the air temperature in the transitional area between the debris-covered and clean ice caused only a small change in total glacier melt, and likewise not replicating the likely increase in wind speed at elevations above UWS had only a small effect. Melt was not found to be particularly sensitive to the value of the downwelling longwave lapse rate chosen. It was found that the values of albedo used caused a reasonable change in modelled melt, so improving how albedo is parameterized should be a future model refinement. Sub-debris melt is also fairly sensitive to the value of the debris thermal conductivity chosen.

The sensitivity of the model to changes in air temperature and debris thickness was investigated. This revealed that a uniform increase in debris thickness of 0.035 m would be required to offset the effects of a 1ºC rise in air temperature, assuming no change in the total area of debris cover. The increase in ablation for a given temperature increase is lower for thicker debris. Thickly debris-covered regions are therefore less sensitive to air temperature increases. The model results show that debris deposition on previously debris-free ice, due to englacial meltout, will increase the surface melt rate if the deposition thickness is less than the debris critical thickness. Hence, up-glacier expansion of debris covers results in an initial acceleration of melt at the upper end of the debris zone, before increasing debris thickness leads to the insulating effect dominating. The above findings are not surprising given the current knowledge of ablation on debris-covered glaciers, but the ability of the model to sufficiently replicate melt suggests it could be used to understand the melt of other debris-covered glaciers, or used within future distributed runoff or mass-balance models. Climate change could imply changes in other meteorological variables (e.g. wind speed and cloudiness) not accounted for in our simple experiment. To account for these, a valid future goal would be to integrate the model into a physical coupled land–atmosphere model (e.g. Mölg and Kaser, 2011; Collier and others, 2013) to allow the effects of future changes in climate to be more closely replicated. The sensitivity of the surface melt rate to variations in the thickness of thin debris covers means that areas of debris in the range 0 to ˜10 cm need to be mapped with greatest accuracy. It is therefore encouraging that techniques utilizing thermal band satellite imagery to map debris thickness distribution on glaciers are most applicable to thin debris (e.g. Foster and others, 2012).

Physically based modelling of a debris-covered glacier does require large amounts of site-specific data, of which multi-station meteorological data are the least likely to be available. The wind field and air temperature distribution are probably strongly influenced by the debris distribution, but the relationship is poorly understood. To understand the debris’s effect on these variables, a detailed study of the variation in air temperature and wind speed across a debris-covered glacier would be required. Quantifying the degree of influence of the debris on the vertical temperature gradient above the glacier (e.g. by measuring temperature at greater heights above the surface) may also aid future modelling. If the model was to be run outside the melt season it might be necessary to account for the refreezing and thawing of water within the debris (Lejeune and others, 2013) and the melt of snow outside the glacier area, neither of which was simulated in the present study. The distribution of dirty ice and the thickness and distribution of debris are also required model inputs, along with values of the debris thermal properties, which were treated simply in this model. The occurrence of ice cliffs was neglected, although they are known to result in locally high melt rates (Diolaiuti and others, 2006; Reid and Brock, 2014). Future model development should focus on both correctly replicating the increase in melt from clean ice to the effective debris thickness and adding routing routines to account for the timing of runoff (obtained from hydro-logical routing studies of a debris-covered glacier). The model could then be used to understand more fully the input meltwater hydrograph on different parts of the glacier and the resulting proglacial runoff signal.