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 Reference Young, Flowers, Berthier and Latto(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 Reference 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, Reference Hill, Dow, Bash and Copland2021) and other small glaciers in the St. Elias mountains (e.g. MacDougall and Flowers, Reference MacDougall and Flowers2011), 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, Reference Janssens and Huybrechts2000). 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, Reference Cuffey and Paterson2010, $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, Reference Janssens and Huybrechts2000; Young and others, Reference Young, Flowers, Berthier and Latto2021a). 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, Reference Huybrechts and de Wolde1999; Janssens and Huybrechts, Reference Janssens and Huybrechts2000) 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 Reference Bliss, Hock and Radić(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, Reference Farinotti, Usselmann, Huss, Bauder and Funk2012; Finger and others, Reference Finger, Vis, Huss and Seibert2015). 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, Reference Huss and Hock2018) 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, Reference Foy, Copland, Zdanowicz, Demuth and Hopkinson2011), 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, Reference Berthier, Schiefer, Clarke, Menounos and Rémy2010; Young and others, Reference Young, Flowers, Berthier and Latto2021a). 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, Reference Mesinger2006). 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 Reference Hock(1999) Distributed Enhanced Temperature-Index Model, which accounts for the effects of topographic shading, slope, and aspect.

3.3.2. Precipitation

Following Young and others Reference Young, Flowers, Berthier and Latto(2021a), NARR precipitation is downscaled using a regression-based approach from Guan and others Reference Guan, Wilson and Xie(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.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, Reference Rounce2021) (Figure S5) but discard the associated critical debris thickness of 13 cm. Studies that have measured the critical debris thickness (e.g. Østrem, Reference Østrem1959; Khan, Reference Khan1989; Mattson, Reference Mattson, Gardner and Young1993; Juen and others, Reference Juen, Mayer, Lambrecht, Han and Liu2014) 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, Reference Loomis1970). 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, Reference Rounce2021). 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 Reference Rounce2021 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 Reference Rounce2021 (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, Reference Rounce2021). 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 Reference Rounce2021. 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, Reference Rounce2021) are nearly identical.

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, Reference Pulwicki, Flowers, Radić and Bingham2018). 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, Reference Li2023). 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.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, Reference Young, Flowers, Berthier and Latto2021a) and (2) the observed snow cover determined by snowline positions delineated from satellite imagery. The geodetic mass balance was determined by Young and others Reference Young, Flowers, Berthier and Latto(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, Reference Young, Flowers, Berthier and Latto2021a) (Figure 5a–c). Simulations where a _ice < a _snow are discarded (e.g. Hock, Reference Hock1999, Reference Hock2003; Young and others, Reference Young, Arendt, Hock and Pettit2018), since snow generally has a higher albedo than bare ice (e.g. Warren, Reference Warren2019) 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, Reference Rounce2021), (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.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 Reference Rounce2021. 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 Reference Rounce2021. 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 Reference Rounce2021 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 Reference Rounce2021 debris model (a)−(c) in (e).

Widespread debris cover over the south lobe of the terminus (Main and others, Reference Main2023) 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 Reference Rounce2021, as both treatments of sub-debris melt are similar over the 20- to 50-cm-thick debris (Rounce and others, Reference Rounce2021) 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