Sensitivity of modelled mass balance and runoff to representations of debris and accumulation on the Kaskawulsh Glacier, Yukon, Canada

Abstract

Runoff contributions from glacierized catchments are changing in response to accelerating mass loss. We reconstruct the 1980–2022 mass balance, runoff and water budget of the ∼70% glacierized Kaskawulsh River headwaters in Yukon, Canada, using an enhanced temperature-index model driven by downscaled and bias-corrected reanalysis data. Debris is treated using melt-scaling factors based on site-specific measurements of the critical debris thickness. Accumulation is estimated from downscaled precipitation bias corrected based on in situ measurements. Model tuning incorporates observations of the 2007–18 geodetic mass balance and seasonal snowline positions on the Kaskawulsh Glacier. We assess model sensitivity to the representation of supraglacial debris and accumulation, including treatments of these processes that can be applied in the absence of in situ data. Different representations of debris produce <1% variation in the catchment-wide runoff and water budget. In contrast, accumulation estimates that omit in situ data produce 33–40% variations in modelled runoff relative to those that use these data. This work identifies site-specific measurements of accumulation as critical to accurate estimates of mass balance and runoff for the Kaskawulsh Glacier, in contrast to site-specific characterization of the effects of debris which influence estimated thinning rates at the glacier terminus but have little impact on the glacier-wide runoff.

1. Introduction

The downstream hydrological effects of glacier mass loss impact important river systems around the world (e.g. Huss, 2011; Bliss and others, 2014; Huss and Hock, 2018; Chesnokova and others, 2020). In glacierized basins, ice melt exerts an influence on the timing and magnitude of downstream discharge (e.g. Neal and others, 2010; Farinotti and others, 2012; Addor and others, 2014; Valentin and others, 2018) and the physical and chemical characteristics of proglacial streams (e.g. Hood and Berner, 2009), impacting freshwater and near-shore marine ecosystems (e.g. Pitman and others, 2021). Concern for water resources is also mounting in many regions of the world as thinning rates of glaciers outside of the Antarctic and Greenland ice sheets have doubled in recent decades (Hugonnet and others, 2021), and current mass-loss rates suggest that many small glaciers, especially those at mid-latitudes, may disappear entirely by the end of the century (Zemp and others, 2019; Rounce and others, 2023). Quantifying the contributions of glacier melt to catchment-wide water budgets and assessing long-term trends in glacier melt are therefore important, especially as discharge regimes change in response to sustained mass loss (Huss and Hock, 2018). Reconstructing long-term glacier runoff records is challenging in part due to the fact that many catchments in remote, mountainous environments are ungauged. In the absence of in situ discharge measurements, observations of glacier mass change derived from remote sensing products such as Digitial Elevation Models (DEMs) (e.g. Berthier and others, 2010; Foy and others, 2011; Moore and others, 2020; Young and others, 2021a) can be used to estimate the meltwater produced by glacier wastage (La Frenierre and Mark, 2014). Others have employed distributed glacier mass-balance and hydrological models (e.g. Farinotti and others, 2012; Immerzeel and others, 2012; Bliss and others, 2014; Li and others, 2020) to partition sources of runoff and estimate the glacier contribution to catchment-wide discharge. Model challenges persist, however, and generally include high uncertainties in input data as well as observations insufficient to constrain model parameters (van Tiel and others, 2020).

Here, we use a distributed mass-balance model to reconstruct the runoff and water budget of a highly glacierized, ungauged catchment in southwest Yukon. We examine how the use of in situ observations to parameterize and tune the mass-balance model influences the estimated runoff and water budget compared to alternative parameterizations that omit glacier-specific information and could be applied in data-scarce catchments. In particular, we assess model sensitivity to (1) the representation of supraglacial debris and (2) the accumulation bias correction. Debris on a glacier surface can either enhance or inhibit melt, depending on the critical debris thickness (Østrem, 1959). The representation of debris in mass-balance models has been shown to influence estimated sub-debris ablation rates and mass-balance gradients (e.g. Juen and others, 2014; Rounce and others, 2021; Compagno and others, 2022). Accumulation inputs also generally represent large sources of uncertainty in glacier mass-balance models (e.g. Machguth and others, 2009; Tarasova and others, 2016), with model performance depending strongly on the availability of observational data (e.g. Immerzeel and others, 2014). We further assess the sensitivity of the estimated water budget to sources of tuning data including the glacier-wide geodetic mass balance and distributed snowlines delineated from satellite images.

2. Study area

The Kaskawulsh Glacier catchment, which we refer to as the Kaskawulsh River headwaters (Figure 1), is a highly glacierized region located within the Traditional Territories of the Kluane, Champagne & Aishihik, and White River First Nations, in the St. Elias Mountains of Yukon, Canada. The catchment is 1704 km², and ∼70% glacierized over an elevation range of approximately 750–3500 m a.s.l. The Kaskawulsh Glacier itself is a 70 km-long valley glacier representing ∼9% of the glacier-ice volume in the Yukon (Farinotti and others, 2019). The debris-covered terminus marks a drainage divide between the Yukon and Alsek River watersheds and is the site of a recent drainage reorganization in which meltwater that previously drained to the Bering Sea was abruptly rerouted to the Gulf of Alaska, resulting in decreased discharge to the Ä’äy Chù (Slims River) and reduced water levels in Łhù’ààn Mân (Kluane Lake) (Shugar and others, 2017). Recent estimates suggest the Kaskawulsh Glacier lost mass at an average rate of −0.46 ± 0.17 m w.e. a⁻¹ between 2007 and 2018 (Young and others, 2021a), nearly matching the regional mass loss rate estimated for the St. Elias Mountains as a whole (Berthier and others, 2010). Mass loss in the catchment is expected to accelerate in the future as temperatures rise in southwest Yukon, which has already experienced more warming than nearly all other regions in Canada (Bush and Lemmen, 2019). Even under a stable climate, however, estimated ice fluxes on the Kaskawulsh Glacier suggest that the glacier is still in the early stages of dynamic adjustment to sustained mass loss over the last several decades, with a minimum committed terminus retreat of 23 km estimated under the 2007–18 climate (Young and others, 2021a).

Figure 1. Study area (blue star, inset upper right) located within the Traditional Territories of the Kluane, Champagne & Aishihik, and White River First Nations. Blue shading indicates the glacierized area, with major tributaries of the Kaskawulsh Glacier labelled: North Arm (NA), Central Arm (CA), Stairway Glacier (SW), South Arm (SA). The regional inset at bottom left shows the locations of two Environment and Climate Change Canada (ECCC) weather stations (cyan circles) located in Burwash Landing (BL) and Haines Junction (HJ). Basemap sources: Esri, Maxar, Earthstar Geographics, and the GIS User Community.

3. Mass-balance model

The distributed mass-balance model used in this study is adapted from Young and others (2021a), and described only briefly here. Changes to the model introduced in this study include an annually adjusted surface-elevation scheme and use of distributed snowline observations in the model tuning procedure (see Robinson, 2024). We also introduce revised parameterizations of debris-covered ice ablation and snow accumulation, described in Section 4 and Section 5, respectively.

3.1. Model description

The mass-balance model calculates the distributed climatic mass balance $\dot{b}_{\rm sfc}(x,y)$ on a 200-m grid spacing with a 3-h timestep as

(1)\begin{equation} \dot{b}_{\rm sfc}(x,y) = \dot{c}_{\rm sfc}(x,y) - \dot{a}_{\rm sfc}(x,y), \end{equation}

where $\dot{c}_{\rm sfc}(x,y)$ is the distributed surface accumulation and $\dot{a}_{\rm sfc}(x,y)$ is the distributed surface ablation. For the accumulation component, this study builds on the work of Young and others (2021a) who developed an elevation-dependent accumulation bias correction for the Kaskawulsh Glacier based on in situ data from the Kaskawulsh River headwaters and neighbouring catchments, which is refined in this study to improve accuracy for this specific catchment (Section 5).

Ablation is approximated as the surface melt (M; m w.e.), calculated using the enhanced temperature-index model of Hock (1999),

(2)\begin{equation} M(x,y) = \begin{cases} (\text{MF} + a_{\rm snow/ice}I(x,y))T(x,y) & \text{if}\ T \gt 0{^\circ}\textrm{C} \\ 0 & \text{if}\ T\leqslant 0{^\circ}\textrm{C}, \end{cases} \end{equation}

where $T(x,y)$ ( $^\circ$C) is the distributed air temperature and $I(x,y)$ is the distributed potential direct clear-sky solar radiation (W m⁻²). MF (m w.e. 3 hr⁻¹ $^\circ$C⁻¹), a _snow and a _ice (m w.e. 3 hr⁻¹ $^\circ$C⁻¹ m² W⁻¹) are, respectively, the melt factor and radiation factors for snow and ice that are empirically determined during the tuning process. These radiation factors differ for snow and ice and vary inversely with albedo, such that a _ice > a _snow. While physically based energy-balance modelling approaches have been previously applied to both the Kaskawulsh Glacier (e.g. Hill and others, 2021) and other small glaciers in the St. Elias mountains (e.g. MacDougall and Flowers, 2011), these methods are generally data-intensive and limited to short time periods with point-scale calibration and validation data. In contrast, this study calculates surface melt using an enhanced temperature-index model, which has less extensive data requirements and is better suited for fully-distributed modelling over a multi-decadal period in the data-scarce Kaskawulsh River headwaters.

The refreezing process is accounted for using a thermodynamic parameterization to estimate the total amount of liquid water (from snowmelt or rainfall) that can be retained by percolation and refreezing in the snowpack, referred to as the total potential retention mass $P_{\tau}(x,y)$ (m w.e.) (Janssens and Huybrechts, 2000). P_τ in each gridcell is approximated as a proportion ( $P_{\rm r}(x,y)$) of the distributed annual precipitation in a given hydrological year ( $P_{\rm annual}(x,y)$; m w.e.):

(3)\begin{equation} P_{\rm r}(x,y) = \frac{c}{L}|{\rm min}(T_{\rm mean}(x,y),0)|\frac{d}{P_{\rm mean}(x,y)}, \end{equation}

where c (2097 J kg⁻¹ K⁻¹) is the specific heat capacity of ice, L (333.5 kJ kg⁻¹) is the latent heat of fusion Cuffey and Paterson, 2010, $T_{\rm mean}(x,y)$ is the local mean annual air temperature for a given hydrological year, $P_{\rm mean}(x,y)$ (m w.e.) is the mean annual precipitation over the whole study period (1980–2022) and d is a prescribed thickness of the thermal active layer, set to 2 m (Janssens and Huybrechts, 2000; Young and others, 2021a). The maximum allowable value of the retention fraction P_r is 1, and therefore the maximum possible potential retention mass P_τ is equal to the annual precipitation (P _annual), since

(4)\begin{equation} P_{\tau}(x,y) = P_r(x,y)\,P_{\rm annual}(x,y). \end{equation}

While P $_{\tau}(x,y)\, \gt \,0$, any melt that occurs is assumed to refreeze, and therefore the maximum amount of refreezing that can occur is capped at $P_{\tau}(x,y)$. Once the upper limit of P $_{\tau}(x,y)$ has been reached, any additional snowmelt or rainfall is assumed to run off (Huybrechts and De Wolde, 1999; Janssens and Huybrechts, 2000) until $P_{\tau}(x,y)$ is renewed at the beginning of the next hydrological year. Therefore, the amount of water that is refrozen ( $R(x,y)$; m w.e.) is related to the available meltwater ( $M_{\rm snow}(x,y)$) and the potential retention mass ( $P_{\tau}(x,y)$) in each gridcell and at each 3-hourly timestep by

(5)\begin{equation} R(x,y) = \begin{cases} M_{\rm snow}(x,y) & \text{if}\ P_{\tau}(x,y)\,\geq\,M_{\rm snow}(x,y) \\ P_{\tau}(x,y) & \text{if}\ 0\,\leq\,P_{\tau}(x,y)\, \lt \,M_{\rm snow}(x,y). \end{cases} \end{equation}

We follow Bliss and others (2014) in defining glacier runoff, Q_g, as the sum of all sources of runoff over the glacierized area:

(6)\begin{align} Q_g(x,y) & = M_{\rm glacier\,\,ice}(x,y) + M_{\rm snow}(x,y)\nonumber \\ & \quad + M_{\rm refrozen\,\,snowmelt/rain}(x,y) + P_l(x,y) - R(x,y), \end{align}

including glacier ice melt ( $M_{\rm glacier\,\,ice}$), snowmelt (M _snow), ice melt from the refrozen snowmelt/rain layers formed during a previous refreezing event ( $M_{\rm refrozen\,\,snowmelt/rain}$) and rainfall (P_l) minus the snowmelt and rainfall that is refrozen (R). The total catchment runoff is the sum of glacier runoff and runoff from the non-glacierized area. Snowmelt, rainfall and refreezing are treated the same over the non-glacierized area as the glacierized area. Losses from groundwater infiltration and evapotranspiration are neglected. We make the simplifying assumption that all runoff instantaneously exits the catchment, and do not incorporate a meltwater routing module (e.g. Farinotti and others, 2012; Finger and others, 2015). Modelled discharge therefore does not account for runoff transit times, groundwater, supraglacial ponding, or englacial storage, which would delay or reduce the estimated discharge. However, for our purpose of examining how the use of in situ observations to parameterize and tune the mass-balance model influences the estimated runoff and water budget, this simple estimation of runoff is sufficient.

3.2. Catchment geometry

Delineation of the glacierized area within the catchment is based on outlines from the Global Land Ice Measurements from Space inventory (GLIMS) Randolph Glacier Inventory (RGI 6.0) (RGI Consortium, 2017) (Kaskawulsh Glacier RGI ID: 60-01.16201). The use of a constant glacier area over time means that the impact on runoff caused by the competition between declining glacier area and accelerating mass loss intensity (e.g. Huss and Hock, 2018) is neglected. However, since the Kaskawulsh Glacier has undergone minimal changes in area in the recent past, with a 1.5% reduction glacier area between 1977–2007 (Foy and others, 2011), neglecting changes in glacier area over 1980–2022 likely has a minimal impact on modelled runoff.

Dynamic surface lowering is accounted for by annually updating the surface elevation of the glacierized area based on a distributed estimate of the average annual elevation-change rate between 1977 and 2018. To generate this estimate, we use DEMs of the study area from 1977, 2007 and 2018 (Berthier and others, 2010; Young and others, 2021a). We calculate the time-weighted average annual elevation change on the Kaskawulsh Glacier between the periods 1977–2007 and 2007–18. We generate a smoothed annual elevation-change map for 1977–2018 by fitting a curve to the time-weighted mean elevation change between the two periods in 200-m elevation bins (Figure S1). The resulting distributed estimate of annual elevation-change is applied to all glaciers in the catchment to get the distributed surface elevation for each year in the study period prior to 2018. In the absence of DEMs after 2018 we assume that the surface is fixed for the remainder of the study period (2018–22).

3.3. Input data

The temperature and precipitation data used to drive the mass-balance model are obtained by downscaling and bias correcting the North American Regional Reanalysis (NARR) dataset (Mesinger and others, 2006). NARR data are available beginning in 1979 and include gridded outputs for a suite of meteorological variables at 3-hourly timesteps on a 32 km 32 km grid, downscaled to a 200-m grid over the catchment. Potential direct clear-sky solar radiation (I in Equation 2) is calculated using the Hock (1999) Distributed Enhanced Temperature-Index Model, which accounts for the effects of topographic shading, slope, and aspect.

3.3.1. Temperature

We downscale and bias correct NARR temperature data following the approach of Young and others (2021a). Temperature downscaling involves an interpolation scheme from Jarosch and others (2012) in which a linear regression is used to correlate NARR air temperature and geopotential height within the lower layer of the atmosphere. The slope and intercepts of the linear regression are taken as the local lapse rate and sea-level air temperature, respectively, for each NARR grid point. These lapse rates and air temperatures are then bilinearly interpolated across the model domain at the 200-m grid spacing and used to calculate 2-m air temperature at the gridcell elevation. We adopt monthly temperature bias correction factors from Young and others (2021a) based on air temperatures measured on or proximal to the Kaskawulsh Glacier.

3.3.2. Precipitation

Following Young and others (2021a), NARR precipitation is downscaled using a regression-based approach from Guan and others (2009) that relates NARR surface precipitation to the Easting, Northing and elevation of the coarse NARR gridcells (Figure S4). Downscaled precipitation is partitioned into rain and snow using a prescribed temperature threshold of 1 $^\circ$C. Snow accumulation is bias corrected by multiplying downscaled accumulation ( $c_{\rm ds}(x,y,t)$) by an elevation-dependent correction factor C(z):

(7)\begin{equation} c_{\rm bc}(x,y,t) = c_{\rm ds}(x,y,t) \, C(z). \end{equation}

The accumulation bias-correction C(z) is determined from the ratio between measured and downscaled accumulation as a function of elevation (see Section 5).

4. Site-specific treatment of supraglacial debris

4.1. Debris thicknesses on the Kaskawulsh Glacier

We use a distributed estimate of debris thickness (100-m gridcell size) for the Kaskawulsh Glacier from a global dataset (Rounce and others, 2021) (Figure S5) but discard the associated critical debris thickness of 13 cm. Studies that have measured the critical debris thickness (e.g. Østrem, 1959; Khan, 1989; Mattson, 1993; Juen and others, 2014) have found values <5 cm, including a 1966 study on the Kaskawulsh Glacier where measurements indicated a critical debris thickness of approximately 4 cm (Loomis, 1970). Thus, the estimated critical thickness of 13 cm in the global dataset is likely too high and would suggest enhanced melt along the medial moraines (Figure 2d), which are instead observed to be raised above the adjacent clean-ice surface. We use in situ measurements of melt on clean and debris-covered ice to determine a site-specific critical debris thickness with which to correct the sub-debris melt-scaling factors from the global dataset (Rounce and others, 2021). Sub-debris melt-scaling factors are unitless, multiplicative factors that enhance or inhibit the clean-ice melt (Equation 2) depending on the debris thickness.

Figure 2. Overview of field experiment to measure the critical debris thickness and resulting sub-debris melt-scaling factors. Ablation stakes were installed in dirty ice (DI00) and debris-covered ice (DB01–DB04) on 19 July 2022 (a) and measured on 31 August 2022 (b). Measured debris thicknesses and net ablation are listed in Table S1. (c) Relationship between debris thickness and ablation on the Kaskawulsh Glacier. (d) Original sub-debris melt-scaling factors for the Kaskawulsh Glacier from Rounce and others 2021 with a critical thickness of 13 cm. (e) New site-specific sub-debris melt-scaling factors generated using a critical thickness of 1.9 cm, determined from the curve in panel (c).

4.2. Field experiment

Seven ablation stakes were installed on or proximal to the medial moraine at the North Arm–Central Arm confluence (Figure 1): one in clean ice, one in dirty ice (DI00) and five in debris-covered ice (DB01–DB04) (Figure 2a). Circular frames with a diameter of 1 m were installed around the ablation stakes and filled with fine-grained sediment (Figure S7) to control the debris thickness (between 1- and 4-cm-thick debris), with the exception of one stake which was installed on the nearby medial moraine in debris approximately 7 cm thick. Debris thicknesses and stake heights were measured on 19 July 2022 when the stakes were installed and again on 31 August 2022. Stake DB01 had formed a depression in the surface approximately 5±3 cm deep, while stakes DB02, DB03 and DB04 had developed ice-cored debris cones ranging in height from 40±10 cm to 110±30 cm (Figure 2b).

Over the course of the ∼6-week experiment, debris cover within the framed areas thinned due to washout from surface streams and downslope redistribution as the cones developed. Average debris thicknesses from July 19 to August 31 2022 were estimated using a positive degree-day weighted average of the initial and final debris thickness measurements (Table S1). Data from the field experiment were interpolated using a cubic spline to construct a site-specific ‘Østrem curve’, which we then apply to the whole Kaskawulsh Glacier to generate new sub-debris melt-scaling factors (Figure 2c). From this curve, the critical debris thickness was determined to be 1.9±0.7 cm, with maximum melt occurring at a debris thickness of 0.6±0.3 cm. For debris thicknesses outside our measurement range (>5 cm), we adopt the same debris thickness–ablation relationship as Rounce and others 2021 (Figure S8).

4.3. Impact of site-specific sub-debris melt-scaling factors

Our estimate of the critical debris thickness represents a substantial reduction from the estimate of 13 cm in the global debris dataset (Rounce and others, 2021). The new site-specific sub-debris melt-scaling factors predict differential ablation that is more consistent with the observed morphology of the medial moraines. Sub-debris melt is inhibited over roughly 82% of the debris-covered area, compared to 37% melt-inhibited area estimated by Rounce and others 2021. For debris thicker than 35 cm (∼10% of the debris-covered area), the site-specific melt-scaling factors and the melt-scaling factors from the global debris dataset (Rounce and others, 2021) are nearly identical.

5. Site-specific accumulation bias correction

5.1. In situ accumulation measurements

In April/May from 2007 to 2022, 27 sets of measurements of snow depth and density were made at 18 different locations within the Kaskawulsh River headwaters between 1220–2670 m a.s.l. (Figure 3a, Table S2). At each site, snow water equivalent was calculated by integrating discrete density measurements, made with a wedge sampler, over the snowpack depth (e.g., Pulwicki and others, 2018). The mean depth-integrated snow density within the catchment between 2007 and 2022 was 338 kg m⁻³ with a standard deviation of 38 kg m⁻³. Additional estimates of seasonal snow accumulation are available from NASA’s Operation IceBridge (NASA-OIB) airborne radar campaign, which surveyed large portions of the North Arm, Central Arm, and South Arm of the Kaskawulsh Glacier on May 10 2021 (Li and others, 2023). We convert these radar-derived snow depths to snow water equivalent using the mean measured snow density of 338 kg m⁻³.

Figure 3. Overview of the accumulation bias correction. (a) Downscaled, uncorrected NARR annual accumulation for 1980–2022, with in situ measurements from snowpits shown by circles. (b) NARR annual accumulation bias corrected with the site-specific elevation-dependent correction based on the ratio between measured and downscaled accumulation (Equation 7) shown in (c). (d) Comparison of co-located accumulation measurements from NASA’s Operation IceBridge and downscaled NARR accumulation with no bias correction (grey), the new site-specific bias correction in (b) (purple), and a bias correction based on ECCC precipitation-gauge data (blue). Mean Absolute Error (MAE) between measured and modelled accumulation is reported for each.

5.2. Selection of elevation-dependent bias-correction function

The elevation-dependent accumulation bias correction C(z) (Equation 7) is determined from the ratio of observed seasonal snow accumulation to downscaled NARR accumulation (Figure 3a). We generate a suite of potential functional forms for the bias correction by linearly interpolating between values of observed to downscaled accumulation averaged over a range of elevation bins (Figure S9). Co-located measurements of accumulation from the NASA-OIB survey of Kaskawulsh Glacier in May 2021 are compared with downscaled and bias-corrected NARR accumulation on the same date to select the precise functional form of the bias correction (Figure S10): averaging over 450 m elevation bins produced the minimum root mean square error between NASA-OIB-measured accumulation and the downscaled and bias-corrected NARR accumulation (Figure 3c). The resulting elevation-dependent bias-correction function C(z) ranges from 1.27 to 2.43, indicating an underestimation of measured accumulation at all elevations by the downscaled NARR data. For elevations outside the range covered by the in situ data, the value of C(z) is kept uniform and equal to the nearest interpolated value. We take the in situ measurement biases spanning 2007 to 2022 as indicative of the long-term pattern and apply this bias correction to the period 1980–2022 (Figure 3b).

5.3. Bias correction with precipitation-gauge data

We also evaluate the changes in modelled mass balance and runoff under the assumption that no in situ accumulation data exists for the Kaskawulsh River headwaters. In this scenario, we could drive the model with uncorrected downscaled NARR data (Figure 3a) or develop an alternative bias correction based on publicly available precipitation gauge data from ECCC stations. The two closest ECCC stations to the Kaskawulsh River headwaters are ‘Burwash A’, located at 820 m a.s.l. approximately 65 km northwest of the Kaskawulsh Glacier terminus, and ‘Haines Junction YTG’, located at 596 m a.s.l. approximately 59 km east of the terminus (Figure 1). NARR precipitation is downscaled at each of the station locations following the approach described in Section 3.3.2 and compared to measured monthly precipitation at both stations (Figure S13). Monthly correction factors for each gridcell in the model are calculated as the distance-weighted average of the correction factors from the two stations. Downscaled NARR precipitation generally overestimates precipitation measured at the two stations (Figure S14), in contrast to the biases within the catchment where NARR generally underestimates the observed accumulation.

5.4. Impact of accumulation bias correction

The site-specific accumulation bias correction based on snow depth and density measurements from within the catchment increases the catchment-wide mean annual accumulation from 1980 to 2022 by 80% compared to downscaled, uncorrected NARR accumulation (Figure 3a, 3b). This reduces the mean absolute error between the in situ snowpit observations and NARR accumulation from 0.36 m w.e. for the uncorrected data (Figure 3a) to 0.18 m w.e. for the site-specific bias corrected data (Figure 3b). Conversely, the alternative bias correction based on regional precipitation gauge data reduces mean annual accumulation by 25% relative to the uncorrected data. The performance of each representation of accumulation (uncorrected, corrected based on catchment-specific accumulation measurements, corrected based on regional precipitation gauge data) is evaluated for the 2021 accumulation season by comparing against the co-located airborne radar-derived measurements. Relative to uncorrected data, the site-specific bias correction improves the spatial distribution of accumulation in the catchment, reducing the mean absolute error between measured and modelled accumulation by 67% (Figure 3d). The precipitation-gauge bias correction exacerbates the mismatch between measured and modelled accumulation, resulting in a 33% increase in mean absolute error relative to uncorrected data.

6. Model tuning procedure

6.1. Mass balance and snowline targets

The melt model (Equation 2) is tuned to two empirical targets: (1) the 2007–18 glacier-wide geodetic mass balance (Young and others, 2021a) and (2) the observed snow cover determined by snowline positions delineated from satellite imagery. The geodetic mass balance was determined by Young and others (2021a) using DEMs of the glacier surface in 2007 and 2018 derived from SPOT5/6/7 satellite observations.

Snowline positions were delineated by eye from over 50 Landsat-8 and Sentinel-2 satellite images from May to September from 2013 to 2019, with the majority of cloud-free images in June–August. Snowlines were categorized as either upper bounds, marking the boundary above which the surface is continuously snow covered, or lower bounds, marking the boundary below which the surface is completely snow-free (Figure 4a). We delineated separate upper and lower bounds on each of the major tributaries to the Kaskawulsh Glacier for a total of 223 individual snowlines. A rasterized version of the observed snow cover in each satellite image was generated by categorizing each model gridcell as a snow-covered surface, snow-free surface or an intermediate transition zone, depending on the elevation of the gridcell relative to the mean elevation of the upper and lower bounds on each tributary (Figure 4b). An individual image score is calculated for each satellite image by comparing the rasterized observed snow cover (Figure 4b) to modelled snow cover on the model date that matches the date of the satellite image. Individual image scores are calculated as ${N_{\rm matching}}/{N_{\rm gridcells}}$, where N _matching is the number of gridcells where the modelled surface type (snow or ice) matches the rasterized observed surface type on the corresponding date and N _gridcells is the total number of gridcells. Gridcells in the transition zone between upper and lower bounds are excluded from these counts, since the model does not resolve partially snow-covered surfaces. A final ‘snowline score’ is then calculated for each simulation based on a temporally weighted average of individual image scores for each satellite image. The final snowline scores, which indicate how well observed snow coverage in space and time is replicated in the model, are normalized by the score representing a perfect match between modelled and observed snow cover in every satellite image, such that the maximum score is 1.

Figure 4. Snowline delineation and rasterization. (a) Sentinel-2 satellite image of the Kaskawulsh Glacier on 17 July 2016, one of the 51 such satellite images used in snowline delineation. Lower bounds (orange) and upper bounds (blue) of the snow are delineated for each major tributary. (b) Rasterized version of the snow cover in (a), showing bare ice (brown, below the lower bound), snow (blue, above the upper bound) and transition zone (green, between the upper and lower bounds).

6.2. Parameter selection procedure

We initially perform 10,000 simulations using randomly selected combinations of the melt-model parameters MF, a _snow and a _icesampled from independent normal distributions (Young and others, 2021a) (Figure 5a–c). Simulations where a _ice < a _snow are discarded (e.g. Hock, 1999, 2003; Young and others, 2018), since snow generally has a higher albedo than bare ice (e.g. Warren, 2019) and should therefore have a smaller radiation factor. Of the remaining simulations, only those with a modelled mass balance that falls within three standard deviations of the 2007–18 geodetic mass balance, −0.46 ± (3 0.17) m w.e. a⁻¹, are retained and are binned according to their modelled 2007–18 mass balance (Figure 5d). A normal distribution defined by the mean and standard deviation of the geodetic mass balance is imposed on the binned results and scaled such that it encompasses exactly 100 simulations, which are then selected from each bin as those with the highest snowline scores (Figure 5e). This procedure ensures that simulations with the top snowline scores comprise the final ensemble of model simulations and that the ensemble yields a mean modelled 2007–18 average glacier-wide mass balance identical to the observed.

Figure 5. Overview of model tuning procedure. (a–c) 10,000 combinations of a _ice, a _snow and MF (grey bars) are randomly selected from truncated normal distributions (black curves). Parameter combinations that yield a modelled 2007–18 mass balance ( $\dot{B}_{\rm mod}$) within 3 standard deviations of the the 2007–18 geodetic mass balance ( $\dot{B}_{\rm obs}$) (red and light blue bars) and have a _ice ≥ a _snow (light blue bars only) are retained. (d) Simulations that meet the criteria described above are binned according to $\dot{B}_{\rm mod}$ (number of bins is square root of sample size, bin size = 0.041 m w.e. a⁻¹). A normal distribution (black curve) defined by the mean and standard deviation of $\dot{B}_{\rm obs}$ is scaled such that it encompasses exactly 100 simulations, which are selected from each bin on the basis of their snowline scores (navy bars), resulting in the distribution shown in panel (e). Note that the values of a _ice, a _snow and MF shown here are divided by 8 to run with the 3-hourly model timestep, and have units of m w.e. 3 hr⁻¹ $^{\circ}$C⁻¹ m² W⁻¹ ( $a_{\rm ice/snow}$) and m w.e. 3 hr⁻¹ $^{\circ}$C⁻¹ (MF) in the model.

We refer to the tuned mass-balance model with site-specific representations of debris and accumulation (described in the previous sections) as the reference model. The mass-balance model is then retuned following the same procedure to explore alternative treatments of debris or accumulation. These are (1) a debris-free case, (2) using sub-debris melt-scaling factors from a global debris dataset (Rounce and others, 2021), (3) using downscaled, uncorrected NARR accumulation and (4) using a bias correction based on ECCC precipitation-gauge data from outside the catchment (Table S4). In each of the retuned models, only one parameterization (debris or accumulation) is changed at a time.

6.3. Value added analysis

Finally, we test the model sensitivity to the tuning procedure by excluding each of the tuning targets in turn. In each of these tests, we run the mass-balance model with the site-specific representation of debris and accumulation and select the 100 simulation ensemble as described below:

1. (1) Test 1 removes the constraint a _ice ≥ a _snow, but otherwise follows Section 6.2.

2. (2) Test 2 excludes the observed 2007–18 glacier-wide mass balance as a constraint and selects the 100 simulations with the highest snowline scores from those where a _ice ≥ a _snow.

3. (3) Test 3 excludes snowline observations as a constraint. From the simulations where a _ice ≥ a _snow, we randomly sample from the normal distribution on the binned mass balance rather than sampling according to the highest snowline scores.

7. Model results

7.1. Reference mass balance and water budget

From the reference model we estimate that the average 1980–2022 mass balance for the glacierized area was −0.38 ± 0.15 m w.e. a⁻¹ with a mean equilibrium line altitude (ELA) of about 2100 m a.s.l. Modelled thinning rates exceed 9.5 m w.e. a⁻¹ on the northern edge of the Kaskawulsh Glacier terminus where thin debris produces a slight melt enhancement. The distributed mean mass balance (Figure 6a) shows the melt-inhibiting effect of debris over a large portion of the terminus region where lighter shades of orange (debris-covered ice) can be seen adjacent to darker shades of red (debris-free ice). Sinuous patterns corresponding to medial moraines originate at the confluence of Stairway Glacier with the main trunk, and at the confluence of South Arm with the trunk, extending to the debris-covered region of the terminus. The medial moraines are approximately 200–400 m across and exhibit less melt than the surrounding clean ice due to the shielding effect of debris thicker than the estimated critical thickness.

Figure 6. The reference model (a) mass balance (Equation 1) (b), refreezing (Equation 5) and (c) runoff (Equation 6) from 1980 to 2022.

We estimate that the average annual runoff from the Kaskawulsh River headwaters over 1980–2022 was 1.89 ± 0.70 Gt a⁻¹, with peak daily discharge rates of approximately 300 m³ s⁻¹ in early to mid July. A total of 61% of catchment-wide runoff originates from glacier ice melt, while snowmelt contributes 31% (Table 1). Refreezing (Figure 6b) plays an important role in reducing runoff early in the melt season, with approximately 20% of the annual snowmelt refrozen. A fraction of the ice that forms as a result of refreezing snowmelt/rain (∼28%) is later remelted, contributing ∼2% of the annual runoff. At high elevations (> 2900 m a.s.l.), all surface melt is refrozen, and thus no runoff occurs from this zone (Figure 6c), while at lower elevations, the refreezing potential (Equation 4) is generally reached by early August, after which all subsequent snowmelt contributes directly to runoff. Rainfall contributes 6% of the annual runoff and occurs primarily at low elevations in late July and early August.

Table 1. Glacierized area-wide mass balance and catchment-wide discharge for 1980–2022 from the reference model and alternative debris-treatment and accumulation bias-correction models (two each). Uncertainties reported are the standard deviations of the 100 simulations comprising each model ensemble

7.2. Model sensitivity to debris

The modelled glacier-wide mass balance over 1980–2022 is independent of debris treatment, a product of retuning the model to match the geodetic mass balance from 2007 to 2018. However, local ablation rates of both debris-covered and debris-free ice differ considerably at low elevations, with differences decreasing with elevation (Figure 7). The sub-debris ice ablation rate averaged over the debris-covered area is 3.90 m w.e. a⁻¹ using the reference model, increasing to 4.72 m w.e. a⁻¹ for the debris-free model and 5.49 m w.e. a⁻¹ for the model with sub-debris melt-scaling factors from Rounce and others 2021. These differences produce variations in the modelled glacier topography, including inverted moraines that exhibit higher melt rates than the surrounding ice when using sub-debris melt-scaling factors from Rounce and others 2021. Using the site-specific sub-debris melt-scaling factors yields ablation rates up to 3.7 m w.e. a⁻¹ higher over clean ice compared to the medial moraines at similar elevations.

Figure 7. Annual ablation (1980–2022) on the main trunk of the Kaskawulsh Glacier estimated using the reference model (a), debris-free model (b) and Rounce and others 2021 debris model (c). Differences in modelled ablation are shown for the reference model minus the debris-free model (a)−(b) in (d) and the reference model minus the Rounce and others 2021 debris model (a)−(c) in (e).

Widespread debris cover over the south lobe of the terminus (Main and others, 2023) leads to reduced ablation compared to the surrounding clean ice for both the reference model and the model with sub-debris melt-scaling factors from Rounce and others 2021, as both treatments of sub-debris melt are similar over the 20- to 50-cm-thick debris (Rounce and others, 2021) in this zone. Compared to the reference model, neglecting debris produces increased ablation over the debris-covered part of the south lobe by up to 6.5 m w.e. a⁻¹. Despite the local variations in ablation rates between debris treatments, adjustments to the melt-model parameters from retuning compensate for differences in ablation across the catchment. As a result, the catchment-wide runoff and water budget vary by <1% (Table 1).

7.3. Model sensitivity to accumulation bias correction

The reference model has a 1980–2022 average winter balance of 0.74 m w.e a⁻¹ at the end of the accumulation season, while the model with uncorrected accumulation and the model bias corrected with ECCC precipitation-gauge data have, respectively, winter balances of 0.38 m w.e a⁻¹ and 0.29 m w.e a⁻¹ (Figure 8a–c). As a result, net ablation and runoff differ significantly across the three models to compensate for differences in accumulation and achieve the same mass balance as enforced through the tuning procedure. Relative to driving the model with downscaled uncorrected NARR precipitation, bias correcting with site-specific data increases the annual catchment-wide runoff by 44%, while bias correcting with precipitation gauge data reduces runoff by 19%. Peak annual discharge is also sensitive to the accumulation bias correction, varying from ∼200 m³ s⁻¹ in the model with uncorrected accumulation to ∼300 m³ s⁻¹ in the reference model and ∼170 m³ s⁻¹ in the model bias corrected with ECCC precipitation-gauge data (black lines in Figure 8d–f).

Figure 8. Comparison of modelled mass balance and runoff from the reference model (a, d), the model with uncorrected accumulation (b, e) and the model bias corrected with ECCC precipitation-gauge data (c, f). (a–c) Glacier-wide annual accumulation (blue), ablation (red) and cumulative mass balance (black) averaged over 1980–2022. The date where $\dot{B}$ = 0 (printed) is the average onset of net ablation. (d–f) Catchment-wide melt-season daily discharge (m³ s⁻¹) averaged over 1980–2022. Pie chart and percentages represent the fractional contributions to total runoff from each source in legend. Bars on the right y-axis show the annual runoff (Gt a⁻¹) from each source (listed in Table 1). Shading on the time series and annual totals show ± 1 σ of variability in the 100 simulations that comprise each model ensemble.

The estimated water budget across all representations of accumulation varies by < 10% for each component, despite significant changes in runoff magnitude. The tuning procedure ensures the best match between modelled and observed snow cover, leading to little variation in the duration of accumulation/ablation seasons between models and thus little variation in the modelled water budget. Similarly, the ELA and accumulation area ratio (AAR) vary by < 2% across accumulation models.

7.4. Value added analysis

7.4.1. Test 1: Excluding a _ice ≥ a _snow constraint

Retaining simulations where a _ice < a _snow increases the number that fall within the geodetic mass-balance target by 130% (+893) out of the initial 10,000 parameters combinations (Figure 5). However, following the tuning procedure, none of the simulations with a _ice < a _snow are selected for the model ensemble since they yield consistently lower snowline scores than simulations where a _ice ≥ a _snow (Figure 9a). This constraint therefore adds no value beyond what the delineated snowlines offer, as the final ensemble for Test 1 is identical to the reference ensemble. Excluding simulations where a _ice < a _snow (and thus excluding generally lower snowline scores) is a simple means of model improvement in the absence of snowline data.

Figure 9. Summary of results from value added analysis Test 1 (a, d), Test 2 (b, e) and Test 3 (c, f). Note the difference in y-axes scales in panels (a–c). (a–c) Final simulation ensembles (navy blue dots) selected for each test based on the tuning criteria described in Section 6.3. (d–f) Catchment-wide melt-season daily discharge (m³ s⁻¹) averaged over 1980–2022. Pie chart and percentages represent the fractional contributions from each source to total discharge. Bars on the right y-axis show the annual runoff (Gt a⁻¹) from each source in the legend (listed in Table 2).

7.4.2. Test 2: Excluding the geodetic mass balance

Without the 2007–18 mass-balance constraint, the mean snowline score in the final ensemble for Test 2 is the same as the mean snowline score in the reference ensemble, but the modelled mass balances are considerably different, ranging from −4.50 to +0.36 m w.e. a⁻¹ (Figure 9b). Modelled snow cover is well constrained by choosing the best snowline scores, such that the mass balance and runoff differences between the reference model and Test 2 are negligible above the ELA, with catchment-wide snowmelt just 5% less than the reference model (Table 2). Parameters a _snow and MF, which together control snow melt and thus the distributed snow cover, occupy a much narrower range compared to the reference ensemble (Figure 10). Without tuning the model to the observed glacier-wide mass balance, a _ice and thus ice ablation are completely unconstrained, leading to a 103% increase in ice ablation and a mean 1980–2022 mass balance of −1.38 ± 1.15 m w.e. a⁻¹ (Table 2). Mass balance data are thus a critical part of the tuning procedure.

Figure 10. Histograms of the melt-model parameters (a) a _ice, (b) a _snow, and (c) MF that comprise the final ensembles for each value added test. Note that Test 1 is identical to the reference ensemble. The values of a _ice, a _snow, and MF shown here are divided by 8 in the model to be compatible with the 3-hourly model timestep and have units of m w.e. 3 hr⁻¹ $^{\circ}$C⁻¹ m² W⁻¹ ( $a_{\rm ice/snow}$) and m w.e. 3 hr⁻¹ $^{\circ}$C⁻¹ (MF) in the model.

Table 2. Glacier-wide mass balance and catchment-wide discharge for 1980–2022 from the reference model and Tests 2 and 3 of the value added analysis. The results of Test 1 (not shown) are identical to the reference model. The AAR and ELA are also reported

7.4.3. Test 3: Excluding snowline observations

Randomly selecting simulations to populate the normal distribution on the observed mass balance, rather than selecting them based on snowline scores, leads predictably to a greater spread in scores (Figure 9c) and in the range of melt-model parameter values, especially for a _snow and MF (Figure 10). While differences in the long-term glacier-wide mass balance and runoff are minimal between Test 3 and the reference model, neglecting snowline scores produces a 17% increase in discharge from snowmelt and a 4% decrease in discharge from glacier ice melt compared to the reference model. Compared to Test 2, which we assume leads to the best representation of observed snow cover, excluding snowline data from tuning yields a higher mean ELA (+110 m), and a smaller AAR (0.58 vs 0.63) (Table 2). The primary value of including snowline observations in tuning is thus to constrain snowmelt and other parameters related to snow cover, which, in turn, influences the mass balance.

8. Discussion

8.1. Low catchment-scale sensitivity to debris

The site-specific treatment of debris includes a substantial reduction in the critical debris thickness, resulting in widespread reductions in the sub-debris melt-enhancement factors compared to those of Rounce and others 2021. At local scales, the choice of debris parameterization produces considerable variations in modelled ablation and surface topography, particularly in the terminus region (e.g. Compagno and others, 2022). At glacier termini, thick insulating debris can result in inverted ablation gradients (e.g. more ablation upglacier compared to at the terminus) (Rounce and others, 2021) and can inhibit retreat compared to the debris-free scenario (e.g. Compagno and others, 2022). Thick debris in the terminus region of the Kaskawulsh Glacier may be contributing to observed stagnation (e.g. Main and others, 2023) and minimal retreat (e.g. Foy and others, 2011). The complicating effects of debris argue in favour of realistic and glacier-specific representations of debris in models, particularly for future projections of glacier evolution (e.g. Rounce and others, 2021; Compagno and others, 2022).

Despite local variations in ablation on the Kaskawulsh Glacier as a function of debris treatment, the net effect of changing the debris treatment on the modelled water budget is minimal, and tuning the models to the geodetic mass balance forces the net ablation across each debris model to be identical and reduces model sensitivity. In this case, the low sensitivity of the modelled water budget to changes in the debris treatment is due in part to the relatively small fraction of debris cover on the Kaskawulsh Glacier. Debris-covered ice represents 7% of the glacierized area, which is within the typical range for glaciers in the Yukon–Alaska region (5–15%) (Scherler and others, 2018). For a more heavily debris-covered glacier, we would expect the modelled water budget to be more sensitive to the treatment of debris. Supraglacial debris on the Kaskawulsh Glacier could have a more significant influence on mass balance and runoff in the future, as the fraction of debris-covered ice is expected to increase through time due to the lateral expansion of medial moraines, the progressive up-glacier appearance of debris as the ELA rises, and local debris thickening over stagnant termini (e.g. Stefaniak and others, 2021; Compagno and others, 2022).

Tuning to the geodetic mass balance can enable the model to compensate for the inclusion or exclusion of debris, reducing the sensitivity of the estimated mass loss. For instance, Compagno and others 2022 showed that for all glaciers across High Mountain Asia (12–13% debris covered), retuning a glacier-evolution model with and without debris changed the projected mass loss in 2100 by just 1–3%. However, the difference in projected mass loss becomes much more significant for individual glaciers with > 50% debris cover. Conversely, Rounce and others 2021 tune a global glacier evolution model with regional mass-balance data for the debris-present scenario and then conducted simulations without retuning the model for the debris-free scenario, resulting in a 37% reduction in sub-debris ablation globally. While retuning a model when the model structure or physics changes (as is done in this study) reduces model sensitivity, applying a model without retuning (as was done by Rounce and others 2021) facilitates a better process-based understanding of the impact of debris on glacier runoff and mass balance.

8.2. Importance of catchment-specific accumulation data

Gridded reanalysis precipitation products often perform poorly in topographically complex, high-elevation terrain (e.g. Immerzeel and others, 2015; Bannister and others, 2019; Hunter and others, 2020). For the Kaskawulsh Glacier, we find that NARR data generally underestimate accumulation, especially at high elevations. Machguth and others 2009 showed that driving a glacier mass-balance model of the Swiss Alps with downscaled, uncorrected regional climate-model precipitation led to underestimating the mass balance of four Swiss glaciers by 0.25–0.75 m w.e due to systematic biases in the underlying accumulation data. Hydrological models are also frequently driven by interpolated local station data (van Tiel and others, 2020). This study demonstrates that low-elevation station data should be used with caution to estimate precipitation in mountainous terrain, as these stations are often not representative of climate in nearby glacierized catchments and may misrepresent biases in reanalysis precipitation. While our tuning approach reduces model sensitivity to the accumulation bias correction with respect to the net mass balance, there are still significant differences in modelled mass-balance gradients, winter balances, and ablation. These sensitivities necessitate careful treatment of accumulation, especially for studies of glacier dynamics and evolution.

Correctly estimating the total volume of precipitation is one of the most important controls on modelled runoff (e.g. Tarasova and others, 2016), especially for glacierized catchments like the Kaskawulsh River headwaters where most precipitation falls as winter accumulation. More spatially and temporally extensive in situ accumulation observations would thus help improve the accuracy of modelled runoff in this catchment. Here, we assumed a constant relationship between downscaled and measured accumulation over time; however, repeat surveys of accumulation using airborne radar would help quantify the interannual variability in seasonal accumulation and examine the time dependence of the biases in NARR data. Additional observations are also needed to characterize the relationship between accumulation and elevation where observations are sparse (e.g., in the southern tributaries). More broadly, improving estimates of snow water equivalent derived from spaceborne remote-sensing products (e.g. Eppler and Rabus, 2021) is an important avenue for future work, as ground measurements of snow density are still needed in combination with remotely-sensed snow depth to estimate snow water equivalent.

8.3. Value of observational targets in model tuning

Tuning the model to the geodetic mass balance integrates both accumulation and ablation processes (Konz and Seibert, 2010), while the snow lines serve to constrain the timing of runoff from snow and ice melt. Our results highlight, unsurprisingly, the high value that the geodetic mass balance adds to model tuning. Indeed, excluding the geodetic balance from tuning produces ice ablation rates that are largely inconsistent with observations. By contrast, when snowlines are excluded, total ice ablation differed by <5%. However, tuning to the geodetic balance can also lead to compensating errors in modelled ablation if the estimated accumulation is incorrect (e.g. Konz and Seibert, 2010; van Tiel and others, 2020). Including other observational datasets in model tuning, such as point measurements of ablation (e.g. Young and others, 2021a) and accumulation (e.g. Young and others, 2021b), streamflow data (e.g. Konz and Seibert, 2010; Tarasova and others, 2016), and glacial melt extents (e.g. Scher and others, 2021) in addition to the geodetic balance may help reduce compensating errors in the net ablation (e.g. Finger and others, 2015).

An advantage to our tuning approach is that it only uses remote-sensing-derived data, making it more applicable to in situ data-scarce catchments. If data from detailed local studies are not available, however, regional mass-balance datasets (e.g. Hugonnet and others, 2021) can fill this gap (e.g. Rounce and others, 2021; Compagno and others, 2022).

9. Conclusion

This study quantifies the multi-decadal mass balance and runoff from a hydrologically important, highly glacierized ungauged catchment in southwest Yukon, with particular attention to assessing model sensitivity to (1) the treatment of sub-debris melt and (2) the accumulation bias correction. We include in our investigation treatments of these processes that can be applied in the absence of in situ or catchment-specific data.

Treating debris using site-specific sub-debris melt-scaling factors produces variations <1% in the catchment-wide discharge and water budget, compared to neglecting debris or using sub-debris melt-scaling factors from a global dataset. Differences in local ablation rates with various debris treatments are significant, however, over the extensively debris-covered terminus region of the Kaskawulsh Glacier where ablation rates are highest. Although debris-cover represents a small fraction of the glacierized area in the Kaskawulsh River headwaters, accounting for it using site-specific observations may improve estimates of glacier surface evolution and retreat, especially as the terminus nears stagnation and debris cover increases over time (e.g. Stefaniak and others, 2021).

In contrast to the treatment of debris, catchment-wide discharge varies considerably as a function of the accumulation bias correction. Accumulation inputs that omit site-specific observations reduce catchment-wide discharge by 33–40% compared to the site-specific accumulation bias correction. Despite tuning the model to the observed mass balance, major model challenges still include high uncertainties in the input precipitation data which can produce compensating errors in modelled ablation. Improving the spatial coverage of accumulation measurements should thus be a high priority for future in situ data collection efforts in this area and similarly glacierized catchments. Measurements spanning large elevation ranges and multiple accumulation seasons will be of particular help in characterizing the spatial and temporal stability of any bias correction.

Glacier runoff estimates can be critical for understanding downstream changes in water availability, impacts to aquatic ecosystems and landscape evolution. In the case of the Kaskawulsh River headwaters, local and regional glacio-hydrological changes are already producing shifts in the timing and magnitude of freshwater that is delivered to the Gulf of Alaska. There is thus a need for coupled mass-balance and ice-dynamics model projections of the Kaskawulsh Glacier in response to its recent climatic imbalance (Young and others, 2021a). The treatment of debris and accumulation impact important mass-balance parameters that will influence these projections, and our work highlights the value of catchment-specific data in this pursuit.