Application of an improved surface energy balance model to two large valley glaciers in the St. Elias Mountains, Yukon

Abstract

Glacier surficial melt rates are commonly modelled using surface energy balance (SEB) models, with outputs applied to extend point-based mass-balance measurements to regional scales, assess water resource availability, examine supraglacial hydrology and to investigate the relationship between surface melt and ice dynamics. We present an improved SEB model that addresses the primary limitations of existing models by: (1) deriving high-resolution (30 m) surface albedo from Landsat 8 imagery, (2) calculating shadows cast onto the glacier surface by high-relief topography to model incident shortwave radiation, (3) developing an algorithm to map debris sufficiently thick to insulate the glacier surface and (4) presenting a formulation of the SEB model coupled to a subsurface heat conduction model. We drive the model with 6 years of in situ meteorological data from Kaskawulsh Glacier and Nàłùdäy (Lowell) Glacier in the St. Elias Mountains, Yukon, Canada, and validate outputs against in situ measurements. Modelled seasonal melt agrees with observations within 9% across a range of elevations on both glaciers in years with high-quality in situ observations. We recommend applying the model to investigate the impacts of surface melt for individual glaciers when sufficient input data are available.

Introduction

The St. Elias Mountains, located along the Yukon–Alaska border, contain ~33 170 km² of glacial ice (Pfeffer and others, 2014), and are home to one of the largest icefields outside of the polar regions (Clarke and Holdsworth, 2002). Glaciers in the Yukon/Alaska region are experiencing consistently negative mass balances (−72.5 ± 8 Gt a⁻¹ from 2002 to 2019; Cirací and others, 2020), high thinning rates (e.g. −0.4 to −0.6 m w.e. a⁻¹ from 2000 to 2007; Foy and others, 2011)and rapid reductions in areal extent (Barrand and Sharp, 2010). These high thinning rates and strongly negative mass balance have caused the region to become one of the most significant contributors to global sea level rise, with only Arctic Canada expected to surpass the region's contributions in the 21st century (Hock and others, 2019).

The retreat of mountain glaciers has significant local hydrological implications. As glaciers melt and retreat, meltwater runoff from glacierised basins typically increases until a maximum runoff value is reached (‘peak water’), beyond which runoff decreases (Huss and Hock, 2018). The timing and magnitude of peak water relative to current runoff varies globally and regionally. Huss and Hock (2018) showed that ~50% of glacierised basins globally have already passed peak runoff, whereas Chesnokova and others (2020) showed that large, heavily glacierised basins in the Yukon have likely not reached peak water, but some smaller basins likely have. The timing and volume of peak water is important for downstream communities who depend on the water resources provided by mountain glaciers. Therefore, it is important to be able to accurately model past, current and future glacier melt rates to predict water resource availability.

Glacier surface melt rates are commonly modelled using surface energy balance (SEB) models (e.g. Rye and others, 2010; Ebrahimi and Marshall, 2016; Shaw and others, 2016; Fitzpatrick and others, 2017; Noël and others, 2018; Bash and Moorman, 2020). SEB models may be applied at both regional and local (i.e. individual glacier) scales. At the regional scale, SEB models can be used to predict ablation in regional mass-balance models (e.g. Noël and others, 2018). Local scale models are often used to assess the sensitivity of glacier mass balance to climatic variations (e.g. Oerlemans and Fortuin, 1992; Engelhardt and others, 2015; Ebrahimi and Marshall, 2016), and to study in-depth the mass-balance characteristics and melt volumes of individual glaciers (e.g. MacDougall and Flowers, 2011; Wheler and Flowers, 2011; Fitzpatrick and others, 2017; Bash and Moorman, 2020). These glacier-scale models are an important tool to examine the availability of meltwater runoff on the surface of glaciers, which has been linked to ice dynamics (Iken and Bindschadler, 1986; Willis, 1995; Herdes, 2014).

Glacier-scale SEB models have generally been successful at modelling melt rates when validated against in situ observations, typically agreeing within ~10–30%, including those applied to two small alpine glaciers in the Donjek Range of the St. Elias Mountains (MacDougall and Flowers, 2011; Wheler and Flowers, 2011). However, recent models have four primary limitations:

1. (1) The lack of a readily available high-resolution (<500 m) surface albedo product. Models often rely on highly-averaged parameterisations based on modelled snow cover, elevation and mean clean ice albedo (e.g. Brock and others, 2000; MacDougall and Flowers, 2011). Surface albedo is one of the most sensitive parameters in SEB models (Oerlemans and Fortuin, 1992; Ebrahimi and Marshall, 2016), leading to large uncertainties in modelled ablation (more than 35% uncertainty in the SW radiation component; Brock and others, 2000).

2. (2) Models sometimes neglect the possibility that the glacier surface may be shaded by adjacent valley walls and steep topography, or fail to properly account for cast shadows (Olson and Rupper, 2019). For example, this mechanism has been suggested by Thomson and Copland (2017) to be responsible for non-uniform surface lowering on White Glacier, Arctic Canada. Shading has been included by a few energy-balance studies (e.g. Arnold and others, 1996, 2006; Aubry-Wake and others, 2018), but has not yet been combined with a high-resolution, non-parameterised surface albedo to achieve the most accurate representation of SW radiation absorption.

3. (3) Models often neglect to account for supraglacial debris. Sufficiently thick debris cover drastically lowers surface albedo while locally insulating the glacier surface and reducing melt rates (Reznichenko and others, 2010; Steiner and others, 2018). This is an area of considerable current research, with much focus on debris-covered glaciers in the Himalayas (e.g. Kraaijenbrink and others, 2017; Nicholson and Stiperski, 2020; Mölg and others, 2020). Steiner and others (2018) found that debris cover on Lirung Glacier, Nepal, locally reduced energy fluxes transferred to the glacier surface by 10-100%. Without careful treatment of debris cover, SEB models which take spatially variable surface albedo into account will compute elevated melt rates over debris compared to clean ice due to the lower surface albedo, leading to significant model errors over these regions. Melt models for debris covered glaciers highlight the large difference in melt rates over debris-covered and dirty ice. For example, Fyffe and others (2014) found that supraglacial debris cover reduced total melt volumes on Miage Glacier in the French Alps by 60%

4. (4) Models that assume the ice surface is at the melting point may overestimate melt rates by ~10% (Greuell and Konzelmann, 1994; Pellicciotti and others, 2009; Wheler, 2009). A subsurface model is necessary to account for heat conduction through ice to compute surface temperature. Of the models that include a subsurface heat conduction model (SSM), the heat conduction is usually solved independently from the energy balance (e.g. Greuell and Konzelmann, 1994; Pellicciotti and others, 2009; Wheler, 2009; MacDougall and Flowers, 2011). However, these processes are highly coupled, as the energy balance at the surface depends on the surface temperature, controlled by upwards longwave (LW) radiation and sensible heat flux, whereas the subsurface temperature depends on the energy available for warming at the surface along with the deeper thermal gradient driven by long-term ice temperature conditions. More complicated models that account for supraglacial lake formation and firn densification have solved these processes simultaneously (e.g. Buzzard and others, 2018), but this approach has not been adopted by most energy-balance models.

We present an improved distributed SEB model that addresses these four limitations, and therefore allows for more accurate quantification of the spatial distribution of melt. We apply the model to Kaskawulsh Glacier and Nàłùdäy (Lowell) Glacier in the St. Elias Mountains, Yukon, using in situ meteorological data and validate outputs against measured surface ablation. The improved model can be used to investigate the meltwater volumes of individual glaciers with greater accuracy, and has the potential to be upscaled to the entire St. Elias Mountains to quantify the controls on melt at regional scale. Future application of the melt model in combination with a supraglacial meltwater routing model may provide important evidence for constraining the drivers of glacier surges in this region.

Study area

The St. Elias Mountains are a high elevation (up to 5959 m a.s.l.) mountain range located in southwest Yukon and southeast Alaska (Fig. 1). Glaciers in the St. Elias Mountains are diverse, representing a range of sizes, elevations and dynamic behaviours (Clarke and Holdsworth, 2002), including a high concentration of surge-type glaciers (Meier and Post, 1969; Clarke and others, 1986). Surge-type glaciers are characterised by a semi-regular oscillation between two dynamic regimes. In the quiescent phase, glacier flow is slower than the balance velocity and mass builds up in a reservoir zone. In the active phase, flow speeds dramatically increase, typically by more than an order of magnitude, and mass is rapidly transferred down-glacier (Cuffey and Paterson, 2010). Although the exact mechanisms are not fully understood, surges are suggested to be controlled by subglacial hydrology, bed characteristics, thermal conditions (Cuffey and Paterson, 2010) and enthalpy budgets (Sevestre and Benn, 2015; Benn and others, 2019).

Fig. 1. Nàłùdäy (Lowell) Glacier and Kaskawulsh Glacier in the St. Elias Mountains, Yukon. (a) Location of the St. Elias Mountains on the border of Alaska, USA and Yukon, Canada. (b) Location of Kaskawulsh Glacier (KG) and Nàłùdäy Glacier (NG) within the St. Elias Mountains. Google Earth, Image: Landsat/Copernicus, IBCAO. (c) False colour composite image of Nàłùdäy Glacier from 3 August 2015 Landsat 8 scene. Blue dot indicates the location of weather station AWSN1. (d) False colour composite image of Kaskawulsh Glacier from 3 August 2015 Landsat 8 scene. Blue dot indicates the location of co-located weather stations AWSK1 and AWSK2. Red dots in (c) and (d) indicate the locations of surface ablation measurements and on-ice HOBO temperature and relative humidity sensors; orange dots indicate the locations of albedo measurements made by Williamson and others (2016).

Glaciers in this region show a range of surge characteristics. Nàłùdäy Glacier (Bevington and Copland, 2014) and Dań Zhúr (Donjek) Glacier (Kochtitzky and others, 2019) have rapid surge cycles (~10–13 year quiescent phase, ~1–2 year active phase over the past 50 years), whereas Trapridge Glacier has surge phases lasting decades (Clarke and Blake, 1991). Moreover, successive surges on a single glacier show distinct characteristics; Kochtitzky and others (2019) showed that surges on Donjek Glacier have varying patterns of surge initiation and termination. Although more research is necessary to understand the full details of surge mechanisms in the St. Elias Mountains, the quantity and distribution of surface meltwater delivered to the bed is potentially an important control on the surge dynamics (Meier and Post, 1969; Bevington and Copland, 2014; Kochtitzky and others, 2019).

We apply our improved SEB model to Kaskawulsh Glacier (60° 44′37″ N, 138° 57′9″ W) and Nàłùdäy Glacier (60° 19′21″ N, 138^° 27^′42^′′ W), and validate the predicted melt rates to in situ measurements. These glaciers are located ~50 km apart on the eastern slope of the St. Elias Mountains (Fig. 1) and experience a similar climatic regime. Kaskawulsh Glacier is a ~70 km long valley glacier ranging in elevation from ~2500 to ~820 m a.s.l. at the terminus (Flowers and others, 2014) with a complex network of tributary glaciers. Nàłùdäy Glacier is a similarly sized (~65 km long) valley glacier, with an elevation range from ~1500 m a.s.l. in the St. Elias Icefields to ~500 m a.s.l. at the terminus.

Despite their physical and climatic similarities, the glaciers differ in their dynamic regimes. Nàłùdäy Glacier is a surge-type glacier (Clarke and Holdsworth, 2002; Bevington and Copland, 2014), illustrated by extensive looped moraines in the lower regions (Fig. 1). Recently, surges have been occurring more frequently (~12 years compared to ~15 year historical average) but with lower velocities than past surges, with the net result of less terminus displacement (Bevington and Copland, 2014). The last surge of Nàłùdäy Glacier was observed in 2009–10, and with the observed quiescent phase of ~12 years it is expected that another surge is nearly due (in ~2022; Bevington and Copland, 2014). In contrast, there is no evidence that Kaskawulsh Glacier is a surge-type glacier, although some of its tributaries do surge (Clarke and Holdsworth, 2002; Foy and others, 2011). The dynamic behaviour of Kaskawulsh Glacier is driven by meltwater inputs, but not in the cyclical way of Nàłùdäy Glacier. Instead, Kaskawulsh responds similarly to the majority of valley glaciers where seasonal velocity patterns are directly driven by meltwater inputs, with high velocities in the spring (up to 2–3× winter values) during periods of rapid melt of the winter snowpack (Iken and Bindschadler, 1986; Willis, 1995; Herdes, 2014; Altena and Kääb, 2017).

Both of these glaciers exert significant controls on regional hydrology. As recently as 1909, Nàłùdäy Glacier advanced during a surge phase to partially block drainage of the Alsek River, causing a major flood downstream when the ice dam broke (Clague and Rampton, 1982). In spring 2016, Kaskawulsh Glacier switched from draining north through the Slims River to the Bering Sea to a southerly drainage through the Alsek River to the Gulf of Alaska (Shugar and others, 2017).

We have chosen to apply the SEB model to these two glaciers in order to compare meltwater inputs and how this may contribute to their dynamic regime. This study is part of a larger project to explain the dynamic regimes of Kaskawulsh and Nàłùdäy, and by modelling and comparing their surface melt we are taking an important step for this. Comparing the quantity and spatial distribution of meltwater on Kaskawulsh and Nàłùdäy Glaciers may provide evidence to constrain the sources of enthalpy (Benn and others, 2019) that contribute to their contrasting dynamics.

Data and methods

Data

Digital elevation model

We drive the SEB model (described below) with digital elevation models (DEMs) extracted from 32 m resolution v3.0 ArcticDEM mosaic elevation data, where each mosaic tile is constructed from imagery acquired between 25 June 2011 and 17 March 2017 by the DigitalGlobe WorldView-1, WorldView-2 and WorldView-3 satellites (Porter and others, 2018). The DEMs are used to calculate local slope and aspect in order to distribute solar radiation, and to distribute temperature by elevation according to the temperature lapse rate. The DEMs were cropped to contain the ablation zone of both glaciers as well as neighbouring peaks and ridges that are likely to shade the glacier surface at low solar angles (Fig. 2).

Fig. 2. Study area and surface albedo. (a) Kaskawulsh Glacier DEM (ArcticDEM mosaic v3.0) covering the glacier as well as surrounding topography. (b) Albedo map of Kaskawulsh Glacier derived from the mean of five Landsat 8 scenes acquired on 15 July 2014, 18 July 2015, 3 August 2018, 18 August 2018 and 30 August 2019. Red pixels have been identified as debris-insulated by our debris algorithm. (c) Nàłùdäy Glacier DEM (ArcticDEM mosaic v3.0). (d) Albedo map of Nàłùdäy Glacier derived from five Landsat 8 scenes acquired on 15 July 2014, 3 August 2015, 23 July 2017, 8 August 2017 and 30 August 2019, with debris cover as in (b). All panels use UTM zone 7N projection.

We use glacier outlines from the Randolph Glacier Inventory Version 6.0 (RGI Consortium, 2017), and include only the area up to the equilibrium line altitude (ELA) as our model considers bare ice only. This approach neglects any meltwater transport downglacier from above the ELA, and so our melt volumes represent only in situ melt below the ELA. By assessing the location of snowlines in Landsat 8 scenes used to derive surface albedo, we define the ELA to be 2100 m a.s.l. on Kaskawulsh Glacier (previously identified as 1958 m a.s.l. in 2007; Foy and others, 2011), and 1750 m a.s.l. on Nàłùdäy Glacier (previously identified as between 1520 and 1700 m a.s.l. in 2010; Bevington and Copland, 2014).

Albedo

The SEB is highly sensitive to surface albedo (Oerlemans, 1991; Ebrahimi and Marshall, 2016). Parameterisations of surface albedo are able to represent average albedo sufficiently, but lack spatial variability (Brock and others, 2000). Recent SEB models have modelled albedo as a function of snow cover (Rye and others, 2010; Noël and others, 2018), or derived albedo from high-resolution UAV imagery (Bash and Moorman, 2020). Although both approaches can capture variability in surface albedo, the latter has the advantage of being derived from spatially distributed imagery. The use of UAV imagery, however, requires in situ measurement and is not easily scalable. Instead, we derive surface albedo from Landsat 8 imagery.

Naegeli and others (2017) derived surface albedo for two glaciers in the Swiss Alps from both Landsat 8 and Sentinel-2 imagery, and validated the satellite-derived albedo with in situ measurements. Here, we apply a simplified version of their method to derive surface albedo from Landsat 8 imagery, and describe our spatially distributed surface albedo maps. We convert the spectral reflectance from Landsat 8 to a broadband reflectance using the narrow to broadband conversion from Liang and others (2003) in order to approximate the surface albedo:

(1)$$\alpha = 0.356 B_2 + 0.130 B_4 + 0.375 B_5 + 0.085 B_6 + 0.072 B_7 - 0.0018\comma\; $$

where B _n is the surface reflectance in band n. Naegeli and others (2017) showed that this conversion neglects the anisotropy of reflection from snow and ice, and introduces up to a 10% bias in derived albedo values. Due to the complexity of accounting for the anisotropic reflections, we use this albedo directly as an input to the SEB model and account for the bias when reporting our uncertainty in modelled melt volumes. We have chosen to derive surface albedo from Landsat 8 scenes rather than applying the MODIS albedo product directly due to their respective spatial resolutions. With Landsat 8 scenes, we derive surface albedo at 30 m resolution, whereas MODIS is only available at 500 m resolution.

We computed surface albedo for Kaskawulsh and Nàłùdäy from five snow- and cloud-free Landsat 8 scenes. For Kaskawulsh, we used scenes acquired on 15 July 2014, 18 July 2015, 3 August 2018, 18 August 2018, and 30 August 2019. Landsat 8 scenes acquired on 15 July 2014, 3 August 2015, 23 July 2017, 8 August 2017 and 30 August 2019 were used to derive albedo for Nàłùdäy. Although seasonal trends in albedo have been found in other studies (e.g. Brun and others, 2015), where albedo decreases throughout the season as more dirt and debris becomes exposed, we found no evidence of this pattern in the albedo derived from these Landsat scenes. Therefore, we averaged the five scenes to obtain an albedo more representative of the mean. Variation between albedo maps was up to 0.1 (~30%) on Kaskawulsh and 0.03 (~10%) on Nàłùdäy, with no clear annual trend.

Meteorological data

The SEB model is forced by in situ meteorological data from a combination of automated weather stations (AWSs) and shielded Onset HOBO temperature and relative humidity (RH) U23 Pro sensors during the portion of the 2010–14 and 2018 melt seasons when the glacier surface is snow-free in the ablation area. From 2006 to 2018, an AWS (AWSK1) on a nunatak adjacent to Kaskawulsh Glacier (60^° 44^′32^′′ N, 139^° 09^′57^′′ W; Fig. 1) provided temperature, pressure, incoming SW radiation, wind speed and wind gust data. AWSK1 did not measure incoming LW radiation. From 29 August 2017 until the station was decommissioned on 26 July 2018 the wind speed sensor malfunctioned, preventing application of the SEB model during the 2017 melt season. From 2018 to 2020, AWSs on nunataks adjacent to Kaskawulsh Glacier (AWSK2, located beside AWSK1) and Nàłùdäy Glacier (AWSN1; 60^° 18^′45^′′ N, 138^° 33^′36^′′ W) provided temperature, pressure, RH, incoming and outgoing SW and LW radiation, wind speed, wind gust and wind direction data (Fig. 1). Unfortunately, the motherboard of AWSK2 failed on 11 May 2019 and was replaced on 5 September 2019, preventing us from modelling the 2019 melt season.

From 2010 to 2014 four HOBOs on Kaskawulsh Glacier (Lower, Middle, Upper, South Arm; Fig. 1, Table S1) recorded temperature and relative humidity near the glacier surface. However, no data were recorded from 20 September 2010 to 13 August 2011. From 2017 to 2019, three HOBOs on Kaskawulsh and Nàłùdäy provided temperature and RH data (labelled Lower, Middle and Upper stations on each glacier; Fig. 1). Table S1 summarises the temporal availability of HOBO measurements, and Table S2 summarises the available weather station data.

We assume that surface elevations and the positions of on-ice HOBO sensors are static. These assumptions neglect that melt itself changes surface elevations, and that the on-ice sensors are advected down-glacier. Each of the HOBO sensors was co-located with a dual-frequency global positioning system (dGPS) receiver, which recorded mean horizontal velocities of 141.3 to 164.3 m a⁻¹, and vertical velocities of − 0.6 to − 4.5 m a⁻¹, over the period 2010–14 at the Kaskawulsh stations (Herdes, 2014). The stations were reset to their original locations every few years, so were always within ~600 m horizontally and ~20 m vertically of their starting position over our 2010–18 study period. Horizontal and vertical velocities are similar on Nàłùdäy Glacier. The dGPS data from the middle station from 25 August to 11 October 2017 show a mean horizontal velocity of 134 m a⁻¹ and vertical velocity of − 7 m a⁻¹. Over our study period (27 July–15 September, or 51 days), these velocities translate to a horizontal displacement of 19 m and a vertical displacement of − 0.98 m.

Temperature and RH from the HOBOs are used to force the model as they provide measurements close to the ice surface, whereas the weather stations are on nunataks adjacent to the glacier and ~100 m above it. However, the AWS measurements of incoming SW and LW radiation and wind speed are used when available.

From 2010 to 2014 when AWSK1 did not measure incoming LW radiation, we model incoming LW radiation using the Stefan–Boltzmann law:

(2)$$LW_{\rm in} = \varepsilon_{\rm a} \sigma T_{\rm a}^4\comma\; $$

where σ is the Stefan–Boltzmann constant, ɛ _a is the atmospheric emissivity and T _a is the air temperature. Following Ebrahimi and Marshall (2016), we parameterise atmospheric emissivity as a linear function of the relative humidity (RH) and the cloud fraction f. We calibrated the parameterisation using the emissivity derived from incoming LW radiation and air temperature from AWSK2, finding

(3)$$\varepsilon_{\rm a} = 0.6007 + 0.0021{\rm RH} + 0.2510f.$$

We apply the method of Crawford and Duchon (1999) to compute the cloud fraction based on measured incoming SW radiation. Our parameterisation captured most of the variability in the atmospheric emissivity in 2018. We found the standard deviation of the residuals (the difference between the measured LW radiation and modelled radiation using the parameterisation) to be 23 W m^-2, only slightly larger than the range 9–20 W m^-2 reported by Ebrahimi and Marshall (2015) for similar parameterisations.

Lapse rates

The air temperature in the SEB model is corrected for elevation by computing temperature lapse rates (Γ) on each glacier using data from the 2018 on-ice temperature sensors. On Kaskawulsh Glacier, we found that the lapse rate was − 3.98°C km⁻¹ (from 23 July 2018 to 17 March 2019; R ² = 0.87), whereas on Nàłùdäy Glacier the lapse rate was − 3.26^°C km⁻¹ (from 27 July 2018 to 20 November 2018; R ² = 0.89).

In addition to temperature, pressure is modelled by assuming hydrostatic equilibrium, and density is computed according to the ideal gas law for moist air, accounting for variations in temperature and pressure with elevation. The resulting expressions are

(4)$$T_{\rm a} = T_{{\rm a}\comma 0} - \Gamma\Delta z$$

(5)$$P = P_0 - \rho_{\rm a} g \Delta z$$

(6)$$\rho_{\rm a} = {P_0 - p_{\rm v}\over R_{\rm d} T_{\rm a} + g\Delta z} + {\,p_{\rm v}\over R_{\rm v} T_{\rm a}\lpar z\rpar + g\Delta z}\comma\; $$

where Δz is the elevation difference between a point in the DEM and the elevation that the reference state (T _a,0,  P ₀) was measured at, R _d and R _v are the dry air and water vapour ideal gas constants, T _a is the temperature in kelvin and p _v is the temperature dependent vapour pressure.

Surface ablation measurements

We evaluate the performance of the SEB model by comparing modelled ablation rates to in situ surface ablation measurements. Surface ablation measurements are available from two separate periods. From 2010 to 2014, two Judd Communications LLC Ultrasonic Depth Sounders (UDSs) were installed on Kaskawulsh Glacier to automatically record changes in surface height (Fig. 3). The UDSs were initially installed at the upper and south arm stations; the south arm UDS was moved to the lower station in August 2013. The UDS data show significant noise in 2011, 2012 and 2013. In 2010, we have high-quality UDS data at the upper and south arm stations, and in 2014 we have high-quality data at the upper and lower stations. These years are where we have the best ability to quantify model performance.

Fig. 3. Modelled (blue) surface ablation compared to ultrasonic depth sounder (UDS) measurements (red) on Kaskawulsh Glacier, 2010–14 at the (a–e) upper station, (f–h) south arm station and (i, j) lower station. Note the difference in modelled period between years.

Since 2017, time lapse cameras recorded hourly images of ablation stakes marked with stripes every 5 cm on both Nàłùdäy and Kaskawulsh. Images were used to calculate surface ablation at several locations on each glacier. These surface ablation data are used to evaluate model performance (Fig. 4, Table S1).

Fig. 4. Modelled (blue) surface ablation (m) compared to timelapse camera ablation stake measurements (red) from 27 July to 15 September 2018 at the (a) Kaskawulsh lower station, (b) Kaskawulsh middle station, (c) Kaskawulsh upper station, (d) Nàłùdäy lower station and (e) Nàłùdäy middle station. The dashed line in (e) shows the extended melt record for Nàłùdäy Glacier computed from melt observations at the lower station.

On Nàłùdäy Glacier, the middle time lapse camera rotated so that the ablation stake was out of the frame from 8 August 2018 onwards. To obtain a complete melt record, we therefore compute the mean ratio of melt at the middle and lower stations when we have both melt observations, and multiply the melt measurements at the lower station by this ratio. We found the mean melt ratio was 0.91 with standard deviation 0.14, and with correlation between the measurements at the middle and lower stations of 0.998, indicating that this approach provides a robust estimate of ablation at the middle stake location. The extended melt record is shown in Figure 4, and is used in Tables 2 and S4 to quantify model performance.

Glacier surface shading

Shading of the glacier surface by valley walls and neighbouring topography is an important component of high-resolution SEB models (Olson and Rupper, 2019). This mechanism is inherently non-local, requiring the DEM to cover any nearby prominent ridges, and especially up to the ridges on the valley walls containing the glacier, highlighting one advantage of the complete spatial coverage provided by the ArcticDEM mosaic data. It is important to note the distinction between calculating shadows cast by neighbouring topography, which depend on the time of day and day of year, and approximating shading by the sky view factor (SVF). The SVF is the ratio of the sky area that is visible (e.g. unobstructed by surrounding terrain) to the complete half-hemisphere area. Since the SVF is constant over time and does not depend on the solar geometry, using the SVF to compute shading neglects the time dependence of shading and the importance of the aspect of the topography obstructing the view (e.g. the difference between a north- and south-facing cliff adjacent to a glacier). Instead of using the SVF, we directly compute shading for our DEMs.

We implement an algorithm to shade DEMs based on the description by Corripio (2003) and Olson and Rupper (2019). The algorithm computes shadows cast by high-relief topography for a given solar position, which is determined by the time of day, day of year and latitude. The algorithm is fully described by Corripio (2003), so here we provide only a brief description. The algorithm traces solar rays across the DEM, computing the projection of each cell onto a plane perpendicular to the incoming solar rays, labelled the solar plane (see Fig. 6 of Corripio, 2003). Starting from the edge of the DEM closest to the sun, the algorithm traces solar rays across the DEM. A cell is shaded when its projection onto the solar plane is less than any of the previously computed projections along the solar ray. For cells that are determined to be in the shade, the direct incoming SW radiation is set to zero.

Debris insulation

In this section, we present an automated algorithm to delineate regions on the glacier that are insulated by debris cover. Our distributed albedo maps are sufficiently high resolution to capture low albedo over regions of dirty ice and debris cover, and so if we were to neglect the spatial distribution of thick debris on the glacier surfaces we would overestimate melt in regions of thick debris cover.

When debris reaches a critical thickness (e.g. ~0.05 m in the laboratory experiments of Reznichenko and others, 2010), the debris acts to insulate the glacier surface and reduce melt rates (e.g. Reid and Brock, 2010; Reznichenko and others, 2010; Steiner and others, 2018). During the day when energy balance is positive, the debris surface warms and builds a steep thermal gradient within the debris layer, maintaining near freezing temperatures at the debris–ice interface and reducing energy transferred to the glacier surface. Overnight, the debris surface efficiently releases the heat absorbed during the day back to the atmosphere. Kraaijenbrink and others (2018) observed this diurnal cycle on Lirung Glacier in Nepal, showing that debris surface temperatures increase to above 15^°C during the day, but was thick enough to prevent heat conduction through to the ice surface and reduce glacier melt rates directly under the debris by 10–100% (Steiner and others, 2018). High surface temperatures invert near-surface temperature profiles and introduce atmospheric instability over debris during the day, significantly altering turbulent heat fluxes over the debris as well as atmospheric circulation patterns over the entire glacier basin when enough of the surface is covered by debris (Collier and others, 2015; Nicholson and Stiperski, 2020). SEB models have been adapted to accurately model melt rates of debris-covered glaciers (e.g. Reid and Brock, 2010; Fyffe and others, 2014; Shaw and others, 2016), showing good agreement with ablation stake measurements. Since our glaciers are not heavily debris covered (≤7%) and we do not have data on the thickness of debris to apply a heat conduction model within the debris layer, we make the approximation that no melt occurs where the surface is insulated by debris as identified by the algorithm below.

Our algorithm is based on analysing the DEMs and surface albedo maps to find medial moraines, which are the dominant type of debris cover on our study glaciers. The albedo maps allow us to find regions with debris cover or dirty ice, and locally elevated regions in the DEMs show where this debris cover is thick enough to reduce melt rates.

For each point in the DEM, we calculate the second derivative of the surface elevation in the across glacier direction, which provides information on how the surface slope changes in the across-glacier direction (surface curvature):

(7)$$\partial_{{\bf yy}} \eta = {\bf y}^T H {\bf y}\comma\; $$

where H is the Hessian matrix of the elevation η and y is a unit vector in the across-glacier direction. We expect negative (downwards) curvature where the surface has been insulated, representing convex surface, since the debris-insulated surface is higher than its surroundings (Mölg and others, 2020).

Therefore, we identify a point as covered by debris sufficiently thick to insulate the glacier surface if the following three conditions are met:

1. (1) $\alpha \lt \bar \alpha$,

2. (2) α < α _max,

3. (3) ∂_yyη < γ.

Condition (1) finds regions with albedo lower than the average albedo $\bar \alpha$, computed using a 1 km moving average, whereas condition (2) ensures we only consider regions with albedo less than a specified maximum debris albedo threshold α _max = 0.125. Together, these conditions find dirty and debris-covered regions. Condition (3) is used to distinguish regions of dirty ice with enhanced melt from regions of debris cover that are thick enough to insulate the surface and reduce melt.

The surface curvature threshold γ controls the typical magnitude of surface curvature of debris cover. The curvature threshold is always negative, as a negative second derivative implies the region is elevated compared to its surroundings. We use γ as a tuning parameter while keeping the value of γ within the expected range. For a medial moraine of width 100 m and height 10 m, we expect γ ≈ −4 × 10⁻⁴ m⁻¹. γ is tuned so that debris cover matches visible regions of debris in the Landsat 8 scenes and the and the 32 m resolution ArcticDEMs, as well as field observations. In the case of Kaskawulsh and Nàłùdäy the curvature parameter likely differs due to their different dynamic regimes. Each time Nàłùdäy surges, the surface becomes heavily crevassed and fractured, effectively erasing any elevation difference between clean ice and debris-covered ice. Following the surge, differential ablation slowly builds up the elevation difference again. This results in small elevation differences between debris-covered and clean ice compared to Kaskawulsh, where elevation differences are continually enhanced.

When applied to Nàłùdäy Glacier, the algorithm identified the primary medial moraine extending from the junction with Dusty Glacier to the terminus, along with several longitudinal ridges of debris cover on the southern half of the main trunk near the lower and middle stations (Fig. 2). Several debris patches and longitudinal features were also identified along both margins almost up to the elevation of the upper station. Significant debris cover was found to be distributed throughout the terminus region below the lower station, which matches with field observations. Overall, we found 11.9 km² of the 363 km² glacier surface (3.3%) to be insulated by debris. We found the debris algorithm had the best ability to classify debris cover when using a curvature threshold γ = −1.5 × 10⁻⁴ m⁻¹ on Nàłùdäy Glacier. This is in line with our expectation of lower elevation differences between debris-insulated and clean ice due to the surging behaviour of Nàłùdäy.

On Kaskawulsh Glacier, we found that a surface curvature threshold value of γ = −4 × 10⁻⁴ m⁻¹ clearly identified the primary medial moraines originating from the junction of the north and central arms, and the junction with south arm. Similar to Nàłùdäy Glacier, we found debris patches distributed across the terminus region (Fig. 2). In total, we found 26.9 km² (7.0%) of the 385 km² surface of Kaskawulsh Glacier to be insulated by debris. Automatically derived moraine locations matched with the locations of moraines visible in Landsat 8 imagery and with field observations.

Subsurface model

The ice surface temperature is an important input to the SEB model, as neglecting subsurface heat flux may lead to overestimation of total melt by 0.8–10.4%, especially at high elevations (Pellicciotti and others, 2009). Therefore, we include a simple 1-D SSM based on that of Greuell and Konzelmann (1994). The model is further simplified by assuming all energy absorption and melting is in the surface layer (Wheler, 2009; Wheler and Flowers, 2011; Buzzard and others, 2018).

Under these assumptions, the subsurface model may be written in the conservation law form

(8)$$\rho_{\rm i} c_{\rm pi} {\partial T_{\rm i}\over \partial t} + {\partial q\over \partial z} = 0\comma\; $$

where ρ _i is the density of ice, c _pi is the specific heat capacity of ice, T _i is the subsurface ice temperature and z is the depth below the surface. The heat flux q is given by

(9)$$q = - k_{\rm i} {\partial T_{\rm i}\over \partial z}\comma\; $$

where k _i is the heat conductivity of ice. The top boundary condition is the heat flux at the surface, which is equal to the energy available to warm the surface. We partition the total heat flux Q _net into energy used to melt ice, Q _M, and energy used to warm the surface layer, Q _T. The bottom boundary condition is a combination of requiring the bottom temperature to be equal to the 12 m ice temperature and a zero heat flux condition, so that:

(10)$$-k_{\rm i} \left.{\partial T_{\rm i}\over \partial z}\right\vert _{z = 0} = Q_{\rm T}$$

(11)$$T\lpar z = H\rpar = T_{12\, {\rm m}}$$

(12)$$q\lpar z = H\rpar = 0\comma\; $$

where we use the 12 m depth ice temperature T _12 m = −3^°C measured on a small tributary glacier adjacent to Kaskawulsh in September 2008 by Wheler and Flowers (2011), as no more recent temperature measurements exist for Kaskawulsh or Nàłùdäy. Q _T is the energy used to warm the surface, computed according to

(13)$$Q_{\rm T} = \left\{\matrix{ Q_{\rm net} \hfill & T_{\rm s} + \Delta T \lt 0 \hfill \cr \vskip4pt - {h \rho_{\rm i} c_{\rm pi} T_{\rm s}\over \Delta t_{\rm SSM}} \hfill & T_{\rm s}\lt 0\comma\; T_{\rm s} + \Delta T \gt 0 \hfill \cr 0 \hfill & T_{\rm s} \geq 0\comma\; \hfill }\right.$$

with surface temperature T _s, layer thickness h and time step Δt _SSM. The maximum warming potential ΔT is defined as the amount the surface layer would warm within a model time step of length Δt _SSM with total heat flux at the surface Q _net:

(14)$$\Delta T = {\Delta t_{\rm SSM} Q_{\rm net}\over \rho_{\rm i} c_{\rm pi} h}.$$

The energy available for melting is then computed as

(15)$$Q_{\rm M} = Q_{\rm net} - Q_{\rm T}.$$

SEB model

Our energy-balance model is based on that of Ebrahimi and Marshall (2016) and Bash and Moorman (2020). The net energy flux at the surface (Q _net) is computed from the net balance of SW (Q _SW) and LW (Q _LW) radiation, and latent (Q _E) and sensible (Q _H) turbulent heat fluxes:

(16)$$Q_{\rm net} = Q_{\rm SW} + Q_{\rm LW} + Q_{\rm E} + Q_{\rm H}.$$

SW radiation

Incoming SW radiation is modulated by surface slope, aspect and solar geometry. The method to distribute solar radiation is based on, and is functionally equivalent to, the model employed by Bash and Moorman (2020), but has been reformulated in terms of local unit vectors to better integrate with the DEM shading algorithm (Corripio, 2003). The unit solar vector s is computed based on the time of day and day of year (Corripio, 2003). The local incident SW radiation I ^′ is then computed from the global incident radiation I ₀, where shaded cells have I ₀ = 0, as

(17)$$I^\prime = \lpar {\bf n}\cdot{\bf s}\rpar I_0\comma\; $$

where n is the upward unit normal perpendicular to the glacier surface. Following Bash and Moorman (2020), we add a diffuse radiation component to all cells:

(18)$$I_{\rm diff} = 16 \psi^{1/2} - 0.4 \psi\comma\; $$

where ψ is the solar angle of elevation in degrees. The net SW radiation flux is

(19)$$Q_{\rm SW} = \lpar 1-\alpha\rpar I^\prime + I_{\rm diff}\comma\; $$

where α is the local surface albedo derived from Landsat 8 scenes.

LW radiation

Outgoing LW radiation is computed from the surface emissivity and surface temperature according to the Stefan–Boltzmann law:

(20)$${\rm LW}_{\rm out} = \epsilon_{\rm s} \sigma T_{\rm s}^4\comma\; $$

where σ is the Stefan–Boltzmann constant, ɛ _s is the ice emissivity and T _s is the surface temperature modelled by the subsurface model. Incoming LW radiation from weather stations is used where available (AWSK2 and AWSN1), and is modelled for AWSK1 (Eqn (2)). Incoming LW is distributed according to air temperature:

(21)$${\rm LW}_{\rm in} = {\rm LW}_{{\rm in}0}{T_{\rm a}^4\over T_{{\rm a}\comma 0}^4}\comma\; $$

where LW _in0 is the LW radiation measured by the AWS, T _a,0 is the temperature recorded by the AWS and T _a is the distributed air temperature (Eqn (4)). The relative change in incoming LW radiation is quite small since the temperature only varies by a few degrees across the glacier. With the largest lapse rate that we calculated (− 3.98^°C km⁻¹), the difference in temperature across Kaskawulsh Glacier is ~ 3^°C, corresponding to a change in incoming LW radiation of $\sim \!4.5\%$ or ~ 13 W m⁻². The net LW radiation is simply

(22)$$Q_{\rm LW} = {\rm LW}_{\rm in} - {\rm LW}_{\rm out}.$$

Turbulent heat fluxes

The sensible heat flux Q _H and latent heat flux Q _E are

(23)$$Q_{\rm H} = \rho_{\rm a} c_{\rm p} k^2 U \left({T_{\rm a}\lpar z\rpar - T_{\rm s}\over \ln\lpar z/z_0\rpar \ln\lpar z/z_{0H })} \right)$$

(24)$$Q_{\rm E} = \rho_{\rm a} L_{\rm v} k^2 U \left({q_{\rm a}\lpar z\rpar - q_{\rm s}\over \ln\lpar z/z_0\rpar \ln\lpar z/z_{0E})}\right)\comma\; $$

where c _p is the constant pressure heat capacity of air, k is von Karman's constant, L _v is the latent heat of vapourisation of water and z ₀, z _0H, z _0E are the momentum, heat and moisture roughness lengths. The parameter z is the height above the glacier surface of the air temperature and specific humidity measurements (Table S3). T _a(z), q _a(z) are the air temperature and specific humidity measured by on-ice HOBOs at a height z above the glacier surface, after distributing quantities by elevation according to the lapse rates, U is the wind speed measured at the AWSs, q _s is the specific humidity at the glacier surface and ρ_a is the density of air.

The SEB (Eqn (16)) depends on the surface ice temperature, while the temperature evolution within the ice depends on the SEB through the boundary conditions (Eqn (10)). We avoid an iterative scheme (e.g. Greuell and Konzelmann, 1994; Wheler, 2009; Wheler and Flowers, 2011) by formulating the models as a coupled system of differential equations (e.g. Buzzard and others, 2018). Differentiating the SEB model (Eqn (16)) with respect to time, assuming that external forcing is constant within each time step, we find an evolution equation for the net heat flux at the surface that depends on the surface temperature. We therefore combine the SEB and SSM models into the coupled system

(25)$$0 = \rho_{\rm i} c_{\rm pi} {\partial T_{\rm i}\over \partial t} + {\partial q\over \partial z}$$

(26)$${{\rm d}Q_{\rm net}\over {\rm d}t} = -4 \epsilon _{\rm s} T_{\rm s}^3 {\partial T_{\rm s}\over \partial t} - {\rho_{\rm a} c_{\rm p} k^2 u\over \ln\lpar z/z_0\rpar \ln\lpar z/z_{0H}\rpar }{\partial T_{\rm s}\over \partial t}.$$

The complete system is given by Eqns (25) and (26), combined with (9) and boundary conditions (10)–(12).

Model implementation

The model is implemented with two different time steps. The first time step, Δt _AWS, is set by the frequency of the AWS data (Table S3). For AWSK1, this time step is 1 hour; for AWSK2 and AWSN1 the time step is 2 hours. The second time step, Δt _SSM, corresponds to the subsurface model for which we use 15 min. During each time step Δt _AWS we assume the meteorological variables are constant. We then solve the coupled SEB and subsurface model equations (Eqns (25) and (26)) with a time step Δt _SSM, and compute the average surface temperature, heat used for warming the ice, and heat used for melting ice over the long time step Δt _AWS. These average quantities are used to compute the total amount of melt during the time step Δt _AWS. For instance, for AWSK2 and AWSN1 we take eight 15-min time steps in the subsurface model for each 2 hour SEB interval. At the end of each AWS time step, the average melting heat flux Q _M is used to melt ice in the surface layer.

Spatial derivatives are implemented using finite differences on a uniform grid. The grid extends down to 12 m depth with a uniform layer thickness of 1 m. Following Buzzard and others (2018) we assume that all SW radiation is absorbed by the surface layer. This approach neglects the fact that SW radiation exponentially decays with depth, warming and melting the subsurface layers (Greuell and Konzelmann, 1994). However, with our layer thickness of 1 m, only ~5% of SW radiation would penetrate to the second layer (Greuell and Konzelmann, 1994). Therefore, we believe that this approximation is valid considering the simplifications it allows in the model formulation.

We tested the sensitivity to the vertical grid spacing by reducing the layer depth from 1 to 0.5 m and reducing the time step from 15 to 5 min. The maximum absolute difference in melt volumes was <1%, and so we believe our time step and layer depth are sufficiently small to resolve the subsurface thermal structure of these two glaciers.

We model surface melt on Kaskawulsh Glacier starting from the time that the glacier surface is snow-free at the upper station (~1700 m a.s.l.) By analysing 12 Landsat 8 images from June through August 2014–19, we determined that Kaskawulsh Glacier is typically snow free at the upper station by 23 June (Fig. S1). The model is run from this date until the end of the melt season, which we define as 15 September, for a total season length of 85 days, subject to data availability constraints. This date for the end of the melt season is in agreement with Herdes (2014), who found the melt season ended as early as 11 September at the upper station of Kaskawulsh in 2011. We only have AWS data on Nàłùdäy Glacier beginning in late July of 2018, by which time it was snow free and so we run the model from this date until 15 September. Each 85-day model run takes ~12 hours on an Intel^® Core™ i5-6300U CPU with 8GB RAM.

Model evaluation and uncertainty estimation

We evaluate model performance quantitatively by computing the total model error (ME) between modelled and measured ablation (using UDS data or time-lapse ablation stake measurements) at the end of the melt season and the mean ablation rate error (ARE) at each station. The ARE is defined as the difference between the slope of the best-fit lines through the modelled and observed melt timeseries. The ME and ARE are complementary metrics in that they are each robust with respect to different types of errors. ME is not impacted by errors in the middle of the melt season, and is controlled by the total error at the end of the model run. In contrast, ARE measures the melt trend throughout the entire modelled period, and is not significantly impacted by small variations in the first and last few days of the season. We report these metrics at each measurement location by extracting the model values in the cell containing the station.

We estimate uncertainty in modelled melt volumes based on uncertainties in the input data. In particular, we account for uncertainties in our derived surface albedo maps, our modelled incoming LW radiation for AWSK1, our wind speed measurements (since they do not reflect measurements directly on the ice) and our sub-surface ice temperature.

The uncertainty in surface albedo is a result of our simple method to derive albedo from surface reflectance data. Naegeli and others (2017) showed that neglecting the anisotropy correction when deriving surface albedo from Landsat 8 surface reflectance can result in up to a 10% bias in derived surface albedo values, depending on surface slope and aspect, as well as solar geometry. Therefore, we perturb our derived albedo values by 10% to derive an uncertainty in modelled melt totals related to our neglecting the anisotropy correction. As albedo varies between DEM cells, the resulting albedo uncertainty is uniquely defined in each cell.

We include an uncertainty contribution from LW radiation when we have had to model incoming radiation. We tuned a parameterisation of emissivity using data from AWSK2, where there was a standard deviation of 23 W m⁻² in the residuals between our modelled and measured incoming LW radiation. Therefore, we assume that there is an uncertainty of 23 W m⁻² in our modelled incoming LW radiation during this time period, and we compute the corresponding change in surface melt, assuming the uncertainty is uniformly distributed across the glacier surface.

Wind near the surface of a melting glacier is typically dominated by a thin layer (< 100 m) of katabatic winds (van den Broeke, 1997; Oerlemans and Grisogono, 2002), and so we expect a non-negligible difference in wind speed between the AWS measurements and the true wind speed near the glacier surface. However, the wind speed profiles have been shown to be consistent through the melt season (Oerlemans and Grisogono, 2002). We see from Eqn (23) that the wind speed only impacts the turbulent heat fluxes. These are also the fluxes impacted by the momentum roughness length. Thus, we consider tuning the roughness length a proxy for adjusting the wind speed. Since we have tuned the roughness length to achieve the best fit with 6 years of in situ measurements, we have accounted for any systematic bias in the wind speed. Therefore, we do not include an uncertainty contribution due to the wind speed.

The final uncertainty contribution is from the ice temperature at 12 m depth. We have used the value reported by Wheler and Flowers (2011) measured on a nearby glacier as no recent data are available on the temperature structure of Kaskawulsh and Nàłùdäy Glaciers. As the measurement was taken at a higher elevation than our study sites, we do not expect the subsurface temperature on Kaskawulsh and Nàłùdäy to be below − 3^°C, and so we assess the sensitivity of modelled melt totals to perturbing the subsurface ice temperature to the maximum value 0^°C. Averaged over the 2012 melt season, the mean difference in net heat flux over Kaskawulsh Glacier was − 2.3 W m⁻², or − 0.052 m w.e. over the period 23 June–15 September. We use this heat flux to compute the uncertainty in modelled melt due to the subsurface ice temperature in each melt season.

We compute the total uncertainty in modelled melt volumes as the combination of the contributions due to albedo ($\delta _\alpha$), deep ice temperature ($\delta _{T_{\rm 12\, m}}$) and LW radiation (when we have modelled incoming LW radiation; δ_LW):

(27)$$\delta = \sqrt{ \delta_\alpha^2 + \delta_{T_{\rm 12\, m}}^2 + \delta_{\rm LW}^2}.$$

We report modelled melt values as melt ± uncertainty (δ) at each measurement location.

Results

Albedo maps

The surface albedo maps (Fig. 2) agree well with existing in situ measurements and with MODIS surface albedo. Williamson and others (2016) measured in situ surface albedo along several transects on Kaskawulsh Glacier (centred on 60.72^° N, 138.82^° W) on 9 August 2015 and Nàłùdäy Glacier (centred on 60.33^° N, 138.61^° W) on 16 August 2015, including along at least one 500 m long transect, and compared these in situ measurements to the MODIS MOD10A1 snow albedo product. We compare these in situ and MODIS albedo values from Williamson and others (2016) to our albedo maps derived from Landsat 8 scenes acquired on 3 August 2015 (6 days before measurements on Kaskawulsh and 13 days before measurements on Nàłùdäy). To approximate the in situ sampling method, we compute the mean, range and standard deviation of our satellite-derived albedo values within a 500 m square centred on the reported coordinate (Table 1). Note that a single 500 m × 500 m MODIS pixel was used for each glacier, so we are unable to report the range or standard deviation.

Table 1. Summary statistics of surface albedo within a single 500 m × 500 m MODIS pixel footprint on Kaskawulsh and Nàłùdäy derived from Landsat 8 scenes on 3 August 2015 compared with values reported by Williamson and others (2016) from MODIS MOD10A1 snow albedo data and in situ measurements

The mean Landsat-derived surface albedo agrees with the MODIS albedo very well, with differences of 4.8% on Kaskawulsh and 0.3% on Nàłùdäy. The standard deviation of the Landsat-derived albedo also agrees well with the in situ measurements. However, the mean and extreme albedo values disagree between the in situ and Landsat methods. This is not surprising since we have more than 250 Landsat pixels within each of the measured plots, while the in situ measurements had only 11 samples on Kaskawulsh and nine on Nàłùdäy (Williamson and others, 2016). Moreover, physical hazards (streams, crevasses, etc.) influenced the locations of the in situ transect but are captured within the remotely sensed data. Considering these limitations, the difference in mean and range between the satellite-derived and in situ measured albedo is not surprising.

Kaskawulsh Glacier, 2010–14

We applied our SEB model to Kaskawulsh Glacier during the 2010–14 melt seasons (Fig. 3). Modelled melt values for the five melt seasons (including three complete melt seasons: 2012–14) and four locations on the glacier are given in Table 2.

Table 2. Modelled melt with estimated uncertainties (m w.e.)

Uncertainties are computed from estimated uncertainty in surface albedo, 12 m depth ice temperature, and incoming LW radiation (2010–14 only) using Eqn (27.)

At the lower station modelled melt averaged over the 2012–14 full 85 day melt seasons was 3.82±0.63 m w.e., or 0.045 m w.e. d⁻¹. At the middle station modelled melt averaged 3.67± 0.62 m w.e., or 0.043 m w.e. d⁻¹. At the upper station modelled melt averaged 3.16± 0.62 m w.e., or 0.037 m w.e. d⁻¹. At the south arm station modelled melt averaged 3.43 ± 0.60 m w.e., or 0.040 m w.e. d⁻¹. Modelled melt was highest in 2013, reaching 4.12 ± 0.63 m w.e. at the lower station, or 0.048 m w.e. d⁻¹. Distributed melt totals for the entire modelled region (Fig. 2) from 2012 to 2014 averaged 3.15 ± 0.61 m w.e. with a standard deviation of 0.84 m w.e., equivalent to a daily melt rate of 0.037 m w.e. d⁻¹.

Agreement between modelled melt and the UDS surface elevation data from 2010 to 2014 is quantified by the ME and ARE in Table S4. We have the best ability to quantify model performance in 2010 and 2014 when the UDS data are high quality. For the 2010 melt season, ME is $-0.08\, {\rm m}\ \lpar -3.0\% \rpar$ and ARE is ${-0.0037}\, {\rm m}\, {\rm w.e.}\ {\rm d}^{-1}\ \lpar -6.8\% \rpar$ at the upper station, whereas ME is 0.52 m (16%) and ARE is $0.0094\, {\rm m}\, {\rm w.e.}\ {\rm d}^{-1}$ (17.5%) at the south arm station. In 2014, ME is ${-0.13}\, {\rm m}\ \lpar -3.6\% \rpar$ and ARE is ${-0.0056}\, {\rm m}\, {\rm w.e.}\ {\rm d}^{-1}\ \lpar -11.7\% \rpar$ at the upper station, whereas ME is 0.29 m (7.1%) at the lower station and ARE is 0.0017 m w.e. d⁻¹ (3.2%) at the lower station. The high ARE at the upper station is at least partially due to the noisy UDS data in mid-July which may be a result of instrument error. In these years, ME is within ~0.30 m (~7%). Modelled melt agrees with measurements within our estimated uncertainty at the lower, upper and south arm stations for all years from 2010 to 2014, with the exception of 2010 along the south arm and 2013 at the upper station.

Kaskawulsh and Nàłùdäy, 2018

The modelled surface melt for the period 27 July–15 September 2018 is shown for both Kaskawulsh and Nàłùdäy Glaciers in Figure 5 and is compared to measurements in Figure 4. Unfortunately, due to various failures in the collection of the in situ ablation measurements we do not have a continuous surface ablation record on Kaskawulsh in 2018. On Nàłùdäy, only the lower station has a continuous melt record, which we have used to approximate measured melt at the middle station.

Fig. 5. Modelled surface melt from 27 July to 15 September 2018 on (a) Kaskawulsh Glacier and (b) Nàłùdäy Glacier. Interior dark regions with no surface melt are those that we have classified as insulated by debris cover.

Modelled melt rates on Kaskawulsh in 2018 fall within the range of melt rates from 2012 to 2014. Modelled melt rates averaged over the period 27 July–15 September 2018 were 0.046 m w.e. d⁻¹ at the lower station, 0.041 m w.e. d⁻¹ at the middle station, 0.035 m w.e. d⁻¹ at the upper station and 0.038 m w.e. d⁻¹ at the south arm station. Modelled melt rates on Nàłùdäy were 0.048 m w.e. d⁻¹ at the lower station, 0.042 m w.e. d⁻¹ at the middle station and 0.034 m w.e. d⁻¹ at the upper station.

The model results are quantitatively compared to time-lapse ablation stake measurements in Table S4. ME is <0.25 m ($\lt \!9\%$) and the ARE is <0.005 m w.e. d⁻¹ ($\lt \!10\%$) at all locations on both glaciers. Modelled ablation agrees with timelapse camera ablation stake measurements within our uncertainty estimates at all locations except at the middle station on Nàłùdäy, where we have had to use melt at the lower station as a proxy to estimate melt at the higher elevation middle station, making it difficult to make conclusions about ME at this location.

Shading by adjacent topography

We expect the relative importance of shading of the glacier surface to be highly sensitive to glacier geometry, aspect and the surrounding topography. We quantify the impact of shading on Kaskawulsh and Nàłùdäy Glaciers in 2018 by calculating the difference in modelled melt at the end of the modelled period (27 July–15 September 2018) with and without including surface shading in the model. Averaged over the entire surface of each glacier, the effect of shading is very small (a total difference of − 0.028 m w.e. or $-1.6\%$, and − 0.020 m w.e. or $-1.2\%$, for Kaskawulsh and Nàłùdäy, respectively; Fig. 6). Due to aspect and neighbouring topography, the impact is generally greater on Kaskawulsh.

Fig. 6. Effect of shading the glacier surface for the period 27 July to 15 September 2018. Difference in modelled surface melt (m w.e.) with and without surface shading on (a) Kaskawulsh Glacier and (b) Nàłùdäy Glacier. Negative values indicate a reduction in surface melt due to shading of the glacier surface. Note the logarithmic colour scales.

Along the main trunk of Kaskawulsh the change in modelled surface melt over the period 27 July–15 September 2018 due to shading ranges from 0 to − 0.05 m w.e. The impact at the lower, middle and upper stations is − 0.02, − 0.001 and − 0.003 m w.e., respectively. This small change is due to the stations being located near the centre of the glacier where the glacier has a predominantly easterly aspect. In comparison, we find large differences in modelled melt along the south arm. The combination of a northerly aspect and a narrower valley leads to differences in melt up to − 0.35 m w.e. (~20%) along the western margin. However, the difference at the south arm HOBO location was only − 0.04 m w.e. (2.0%) since it is located near the centre of the glacier, to the east of the shading footprint.

In comparison, shading of the glacier surface is relatively less important on Nàłùdäy Glacier. The difference is negligible across lower elevations, with a maximum of ~ 0.025 m w.e. where the glacier has a more northerly aspect. There is a large difference (~ 0.25 m w.e., or ~15%) on one of the tributaries, but the difference is more spatially localised than on Kaskawulsh. On Nàłùdäy, the difference at the lower, middle and upper stations is − 0.002, − 0.001 and − 0.02 m w.e., respectively. The upper station has a larger difference since it is located where the glacier has a northeasterly aspect, allowing the surface to be shaded late in the afternoon, whereas the small (< 1 cm w.e.) difference at the lower and middle stations is typical of the main trunk of Nàłùdäy.

Subsurface model

We also compare model outputs with and without including the subsurface model to quantify its importance, using a subsurface ice temperature of − 3^°C. We use the 2012 melt season (23 June–15 September) at Kaskawulsh as a case study because (1) we have complete temporal coverage to run and validate the model, (2) modelled melt compares well to measured melt and (3) air temperatures were abnormally cold in early July 2012, leading to significant warming heat fluxes. When we exclude the SSM, we assume the surface layer of ice is a constant 0^°C.

Throughout the 2012 melt season, the subsurface model reduces modelled melt on Kaskawulsh Glacier by ~0.1 to 0.35 m w.e. compared to the SEB model without the subsurface model, with the difference increasing rapidly with elevation (Fig. 7). At the lower, middle, upper and south arm stations the model difference is − 0.14 m w.e. $\lpar -4.0\% \rpar$, − 0.21 m w.e. ($-6.1\%$), − 0.26 m w.e. ($-8.6\%$) and − 0.20 m w.e. ($-6.2\%$), respectively. There is relatively little variation with surface aspect, slope or albedo (Fig. 7). Instead, the spatial variation in melt differences is driven by the lower air temperatures at higher elevations. Colder temperatures allow the surface to cool further below the melting point overnight, resulting in larger warming heat fluxes the following morning, and therefore reduced melt totals. In 2014, which was a warmer year, we found the difference to be only − 0.19 m w.e. ($-6.3\%$) at the upper station.

Fig. 7. Impact of subsurface model for heat conduction throughout the 2012 melt season (23 June–15 September) on Kaskawulsh Glacier. (a) Difference in modelled melt (m w.e.) between SEB model with subsurface model and SEB model without subsurface model. Negative values indicate that the subsurface model reduced surface melt. (b) Absolute value of difference as a function of elevation.

We also investigated the surface temperatures and warming heat flux throughout the beginning of the 2012 melt season (23 June–23 July). From 30 June to 19 July cold overnight low temperatures at the upper station allowed the surface to cool to as low as − 1^°C (Fig. S2). As the SEB became positive the following day, significant heat, up to ~150 W m⁻², was used to warm the ice rather than melt the surface (Fig. S2). Over this 1 month time period, the models diverged by almost 0.1 m w.e. at the upper station.

The impact of the subsurface model on melt volumes persisted when removing the deep thermal gradient by imposing a subsurface ice temperature of 0^°C. In this case, the difference in modelled melt volumes on Kaskawulsh from 23 June to 15 September 2012 ranged from − 0.04 to 0.32 m w.e. At the upper station, melt was reduced by 0.20 m w.e., while melt was reduced by 0.15 m w.e. at the lower station.

Discussion

Model performance

Our model generally shows very good agreement with surface ablation measurements. Total ME is generally <9% at the end of the melt season, and in years with high-quality in situ measurements is often within ~5%. These errors are reduced significantly compared to the 6–29% errors reported in previous studies in the region (e.g. MacDougall and Flowers, 2011; Wheler and Flowers, 2011) as discussed further below.

The primary exceptions to the good agreement are the 2013 melt season and the south arm station on Kaskawulsh Glacier. In these cases, the model diverges from the measurements near the end of the melt season, possibly due to late-summer snowfall events, particularly at the higher elevation stations, or near- and below-freezing air temperatures altering atmospheric stability near the melting glacier surface affecting the turbulent heat fluxes, as the parameterisations we have used assume a stably stratified layer at the glacier surface (Fitzpatrick and others, 2017).

Model performance is similar on both Kaskawulsh and Nàłùdäy Glaciers. However, the model parameters leading to the highest accuracy did not directly transfer between the glaciers. In tuning the model we found that a momentum roughness length of 3 mm provided the best results on Kaskawulsh Glacier, whereas a momentum roughness length of 5 mm provided the best results on Nàłùdäy Glacier. This finding is in line with MacDougall and Flowers (2011), who found their energy-balance model transferred well temporally but not between glaciers. Our approach to tuning the model using the surface roughness parameter is common in SEB models (e.g. Ebrahimi and Marshall, 2016; Bash and Moorman, 2020). However, MacDougall and Flowers (2011) did not tune their energy-balance model using the surface roughness, instead using directly measured values to more accurately assess the spatial transferability of the model. The surface roughness values we derived by tuning our model are similar to values found in previous studies when using the roughness length as a tuning parameter (e.g. 5 mm by Bash and Moorman, 2020), and to values found by direct measurement on glaciers in the Canadian Rockies (0.7–4.5 mm; Fitzpatrick and others, 2019). In our case, the roughness length may be a reflection of the glacier's dynamic behaviour. The surface of Nàłùdäy Glacier is expected to be more heavily crevassed and fractured following its surge phase, potentially leading to a higher surface roughness (Herzfeld and others, 2000; van der Veen and others, 2009). Alternatively, the differing roughness lengths may be a result of a different bias between the wind speed as measured by the AWSs and the wind speed near the glacier surfaces on each glacier.

Our model performs less robustly on the south arm of Kaskawulsh than along the main trunk. The south arm is a northerly flowing tributary confined in a much narrower valley than the generally easterly/northeasterly flowing main trunk (Fig. 1). We found the surface of the south arm to be heavily shaded by neighbouring terrain, reducing seasonal melt totals by up to ~10% (Fig. 6). Combined with the northerly aspect, the south arm experiences significantly less incident SW radiation than the main trunk. Reduced SW radiation inhibits the growth of a weathering crust and consequently reduces the aerodynamic roughness length (Stevens and others, 2018), leading to reduced turbulent heat fluxes. Since we have used a spatially uniform roughness length, we may be overestimating turbulent heat fluxes at the south arm station on Kaskawulsh Glacier.

Two distributed energy-balance studies have been carried out recently in the St. Elias Mountains, both in the Donjek Range located just north of Kaskawulsh Glacier. Wheler and Flowers (2011) modelled 1.80 m w.e. of melt at 2280 m a.s.l. from 1 May to 13 September 2008 on a small south-facing glacier located 2 km north of Kaskawulsh Glacier, whereas we modelled ~1 m w.e. of melt at ~2000 m a.s.l. (Fig. 8) on a glacier with an easterly aspect, where we would expect lower melt rates. MacDougall and Flowers (2011) modelled ~1.2 m w.e. of melt at ~2300 m a.s.l. on the north-facing North Glacier. Considering the variations in elevation, aspect and temporal extent, our modelled melt rates are comparable to these previous studies. Wheler and Flowers (2011) report relative MEs ranging from 6 to 29%, whereas MacDougall and Flowers (2011) report errors ranging from 10 to 18% on two glaciers in the Donjek Range. Compared to these studies, our errors ranging from 5–10% represent a significant improvement in model performance.

Fig. 8. Modelled melt rates and elevation profiles on Kaskawulsh and Nàłùdäy Glaciers for the period 27 July–15 September 2018. (a) Average melt (m w.e.) at 1000, 1250 and 1500 m on Kaskawulsh (solid) and Nàłùdäy (dashed) Glaciers. (b) Centreline elevation (m a.s.l.) on Kaskawulsh and Nàłùdäy. (c) Average melt (m w.e.) on Kaskawulsh and Nàłùdäy.

We found that the subsurface model (Eqn (25)) is important for accurately modelling melt volumes. At the upper Kaskawulsh station in 2012, melt was reduced by 0.25 m w.e. (8.6%). Near the ELA, melt was reduced by up to 0.37 m w.e. (13%). These differences in melt agree well with the results of Wheler and Flowers (2011) and Pellicciotti and others (2009) who found melt rates were reduced by ~10% in the upper ablation zone, and so we conclude that our implementation of the subsurface model as a system of coupled differential equations is functionally similar to the iterative approach used in these previous studies, but is a more consistent formulation as we allow surface temperature and energy balance to evolve simultaneously. Moreover, we found that the subsurface model significantly reduced melt volumes even for a subsurface temperature of 0^°C, indicating the reduction in melt is due to the cooling of near-surface layers overnight more than the deep thermal gradient. The subsurface model will be even more important for yearly mass-balance models where the persistence of a winter cold wave into the spring may drastically slow the onset of extensive melt, and for cold polar glaciers where subsurface heat flux due to the deep thermal gradient will be non-negligible even during the peak of summer.

We found that shading of the glacier surface reduced glacier-wide melt rates by $\lt \!2\%$, but locally by more than 20%. This is likely a strong enough impact to locally alter the spatial distribution of surface meltwater, and therefore has implications for ice dynamics in terms of volumes and locations of meltwater accessing the subglacial system. Surface shading may also help to explain spatial variations in surface elevation, as we found high elevations in our DEM along the western margin of the Kaskawulsh south arm, coincident with where our model predicts lower ablation rates due to shading. Moreover, we may be underestimating the impact of cast shadows on melt volumes due to the resolution of our DEM (32 m). Olson and others (2019) showed that using a 32 m resolution DEM leads to an average SW radiation that is 7–20% greater when compared to an 8 m resolution DEM, since coarser resolutions act to smooth out sharp elevation changes. Thus, shading becomes especially important in high-resolution melt modelling (e.g. Bash and Moorman, 2020).

Insulation by debris also has a small global impact on melt rates on our study glaciers (~7%) but in localised regions the difference may be enough to impact the distribution of meltwater across the surface. Moreover, the global difference is likely to be much larger on more heavily debris-covered glaciers, for example in High Mountain Asia where ~11% of total glacier area is debris covered, and nearly 100% in some lower ablation areas (Kraaijenbrink and others, 2017), or following rock avalanches (e.g. Shugar and others, 2012).

The final major modification we made in our energy-balance model is deriving high-resolution albedo (30 m resolution) from Landsat 8 surface reflectance data (Naegeli and others, 2017). The resulting albedo maps have a similar range and standard deviation to nearly coincident in situ measurements. Moreover, our model has similar performance across a range of elevation on both Kaskawulsh and Nàłùdäy Glaciers, suggesting we have captured the spatial variation in albedo. Therefore, we believe the Landsat 8 surface reflectance appears to be a highly portable and sufficiently accurate proxy for surface albedo that can be used in SEB models.

Model limitations and uncertainties

We have reported large uncertainties in modelled melt and ablation from 2010 to 2014 (up to 20%) as a consequence of uncertainty in the surface albedo and incoming LW radiation, highlighting the importance of measuring these quantities accurately. The uncertainty in surface albedo could be reduced by directly measuring surface albedo on each of the glaciers in the study area and calibrating the albedo maps to these measurements. The in situ measurements of Williamson and others (2016) provide evidence that our albedo maps are reliable enough to be used without calibration, but more thorough comparison and calibration using in situ data could reduce the uncertainty in melt from SW radiation when a tighter uncertainty bound is required. The uncertainty in incoming LW radiation can be eliminated by installing weather stations that measure LW radiation (e.g. AWSK2 and AWSN1).

The UDS data in Figure 4 illustrate some of the limitations of the depth sounder data. Throughout 2010, 2011 and most of 2014, the data are high quality with low noise levels. However, throughout 2012 and 2013, and particularly at the upper station, the data have significant noise levels; consecutive measurements at 30 min intervals vary by as much as 30 cm. We believe this behaviour is due to instrument error, perhaps due to reflections from off-nadir locations as the station melts out unevenly, which makes it difficult to quantify model performance. Although ME is large in 2013 (− 0.68 m or $-14\%$ at the upper station), it is not clear how to partition the bias between suspected instrument error and ME.

Implications for understanding surge mechanisms

Total modelled melt is similar at Kaskawulsh and Nàłùdäy in 2018. We modelled higher melt rates at coincident elevations on Kaskawulsh compared to Nàłùdäy, particularly in the lowest 10 km of the glaciers, but these higher melt rates were largely offset by the higher elevation range of Kaskawulsh (Fig. 8). The higher melt rates we found on Kaskawulsh at coincident elevations are in part a result of a lower surface albedo on Kaskawulsh compared to Nàłùdäy leading to increased absorption of SW radiation (Fig. S3). Lower surface albedo on Kaskawulsh may be a result of the glacier geometry; Kaskawulsh Glacier has several major confluences (Fig. 1), each providing an opportunity for debris to access the middle of the glacier. Furthermore, since Nàłùdäy is a surge-type glacier, its surface is upheaved and fractured following each surge, exposing regions of subsurface clean ice and leaving the surface relatively clean. In contrast, debris on Kaskawulsh slowly melts out and accumulates on the surface, leading to lower albedo even away from medial moraines and confluences.

We can examine whether volumes of surface water input could impact surge dynamics by considering the enthalpy balance suggested by Benn and others (2019). We have shown similar melt volumes on Kaskawulsh and Nàłùdäy, suggesting that differences in enthalpy between these two glaciers must be in the subsurface environment rather than driven by volumes of surface water input. The surface slope is very similar between the two glaciers (Fig. 8), and combined with similar melt rates, we expect a similar rate of surface steepening due to surface melt, and a similar resulting change in driving stress. The surface slope should increase in the upper regions of Nàłùdäy Glacier leading up to a surge, but at present this does not seem to be impacting the lower portions of this glacier. To fully determine the control of meltwater on dynamics, a supraglacial routing model and subglacial hydrology model (e.g. Werder and others, 2013) will be necessary, and with our SEB model we have excellent inputs to drive such models.

Conclusions

We have presented an energy balance model that captures the melt measured at Nàłùdäy and Kaskawulsh Glaciers, both spatially and temporally, to within 9% (range: − 8.4 to 7.1%), a significant improvement in performance compared to other SEB models that have been applied in this region. In addition to improved performance, we have demonstrated the importance of spatially varying albedo, shading by adjacent topography (reducing melt by more than 20% in localised regions), insulation by thick debris (reducing total melt volumes by ~7% averaged across our study area), and subsurface heat flux (improving agreement between modelled and measured melt by up to 8.6% at the upper Kaskawulsh station) in both the volume and distribution of melt on both glaciers. Together, the improvements we have made to existing SEB models significantly advance the spatial representation of surface melt. The improved model is an important tool in understanding meltwater inputs to glacier systems, and how meltwater influences ice dynamics in the St. Elias Mountains.

Seasonal melt volumes are similar between the two glaciers, with higher melt towards the terminus of Kaskawulsh Glacier due to lower albedo and slightly higher turbulent heat fluxes offset by a higher elevation range compared to Nàłùdäy Glacier. The relatively higher albedo on Nàłùdäy may be due to its surge-driven dynamics that cause extreme fracturing and upheaval every 12–15 years, exposing cleaner subsurface ice. Similarly, we find that regions of debris insulated ice are more elevated on Kaskawulsh Glacier compared to Nàłùdäy, again likely due to the surface topography reorganisation when Nàłùdäy surges.

The model we have presented requires extensive input data. Combined with the relative complexity of the model, it is best-suited to catchment-scale applications where high-quality in situ data are already available to better understand the volume and spatial and temporal distribution of meltwater production throughout the melt season. Since we derive distributed surface albedo from Landsat 8 data, the model also has wide applicability to further investigate the impacts of the darkening of mountain glaciers (Naegeli and Huss, 2017; Di Mauro and others, 2020) and ice sheets (Box and others, 2012; Bond and others, 2013; Dumont and others, 2014; Williamson and others, 2019; Tedstone and others, 2020).

We found some of the model components to have less of an impact on total melt volumes than others. In particular, surface shading made little difference to seasonal melt volumes on both Kaskawulsh and Nàłùdäy. However, this finding is primarily a result of the geometry of these particular glaciers. They generally flow to the east and fill relatively wide valleys. We showed that shading was relatively more important along the South Arm of Kaskawulsh, indicating this is an important mechanism for narrower, northerly aspect glaciers. Therefore, it is important to quantify the impact of each of the model components (shading, subsurface model and insulation by debris) before deciding whether to include or exclude them.

Future research to combine our SEB model with a surface routing algorithm and subglacial hydrology model would further constrain the controls on surging in the region and provide evidence to explain the contrasting dynamics of Kaskawulsh and Nàłùdäy.