We calculate surface ice velocity with velocity mapping software, v m a p (Shean, 2019), using the NASA Ames Stereo Pipeline image correlator (Shean et al., 2016). We use coregistered (see 3.5) consecutive 1 m resolution ALS DEMs re-sampled to 3 m resolution and 3 m resolution PlanetScope imagery orthorectified with the DEMs. We use a kernel size of 51 pixels, an output posting of 5 m, a 1σ Gaussian filter and a Bayes expectation–maximization weighted affine adaptive window correlator (Shean et al., 2016). All image or DEM pairs we use to create velocity fields span one year on average (365 ± 21 days) and have an average glacier area coverage of 89 ± 9% (Table 2). We mosaic all velocity products which span a given year (2–4 pairs) to create 25 m resolution x (N–S) and y (E–W) velocity fields for each year, which are less prone to artifacts, holes and bias that may affect a velocity field from an individual image pair. Mass-balance stake locations were surveyed with an accuracy of ± 0.25 m (Pelto et al., 2019), and we use their displacement to evaluate the velocity fields. Conrad Glacier has 64 total stake velocity pairs, which includes three transverse cross-sections of stakes (Pelto et al., 2019). Kokanee Glacier has 11 total stake velocity pairs, Nordic glacier has 13 pairs, Zillmer Glacier has 23 pairs and Illecillewaet Glacier has no stake velocity data. To investigate whether the widely available NASA MEaSUREs ITS_LIVE ice velocity products (Gardner et al., 2020) can be used in lieu of our velocity fields, we also compare annual ITS_LIVE velocity products to our mosaics for each glacier except Kokanee Glacier, which is outside ITS_LIVE coverage.

We use ice penetrating radar (IPR) measurements of ice thickness from Pelto et al. (2020) for the five glaciers, with uncertainty estimated between 5 and 10% depending on the quality of the bed reflection. We assume an uncertainty of 10% for IPR observations (Pelto et al., 2020). We also use modeled ice thickness from Farinotti et al. (2019), Pelto et al. (2020). For the modeled ice thickness estimates, we double the average relative error from Pelto et al. (2020) for these five glaciers (5.1%) as we are using the ice thickness at-a-point rather than glacier-wide, yielding an uncertainty on ice thickness ( σ H ) of 10.2%. For ice thickness in western North America, Farinotti et al. (2019), hereafter referred to as FAR19, uses three models (Huss and Farinotti, 2012; Frey et al., 2014; Maussion et al., 2019). Pelto et al. (2020) optimized distributed ice thickness for all five study glaciers using a cross-validation approach to minimize error between observed and modeled ice thickness with the Open Global Glacier Model (OGGM) (Maussion et al., 2019). We also use ice thickness from Clarke et al. (2013) to derive ice flux for Conrad Glacier, available at 200 m resolution. A portion of the accumulation zone of Illecillewaet Glacier is assigned to the neighboring Geike Glacier within the Randolph Glacier Inventory (RGI) 6.0 (Pfeffer et al., 2014). Pelto et al. (2020) corrected the outline of Illecillewaet Glacier prior to the surface inversion. We combine the thickness of Geike Glacier (RGI60-02.03686) from FAR19 with that of Illecillewaet Glacier to cover the missing portion for the FAR19 bed estimate of Illecillewaet Glacier. Where bounded by ice-free terrain, we manually delineate the glacier boundary using a ALS DEM hillshade of the previous late summer and a height change DEM (Pelto et al., 2019). For ice divides, assuming that glacier flow is more sensitive to surface slope than bed geometry, the boundary of each glacier is determined using a watershed algorithm with our ALS DEMs.

We correct modeled ice thickness for thinning that occurred between the date of the ice thickness estimates, which all used the Shuttle Radar Topographic Mission DEM and our 2016 LiDAR DEMs as described in Pelto et al. (2020). We also correct modeled and observed ice thickness for thinning which occurs during our three-year study using ALS DEMs.

Our velocity fields (Figure 2) and ice thickness measurements (Figure 3) allow us to calculate ice flux at several cross-glacier “gates” approximately transverse to flow. We sample points every 25 m along each gate, effectively breaking each gate into small segments through which the ice flux is calculated and summed. We use our x and y velocity fields to calculate ice flow direction and magnitude. Due to the curvilinear nature of flow at our glaciers (Figure 2) we take the angle of a given flux gate to calculate the velocity perpendicular to that gate for each segment. To convert this velocity from mean-surface to depth-averaged velocity, we compare seasonal (June–September) to annual stake velocity estimates on Conrad Glacier (the only glacier where we measured seasonal stake velocities) and find that average seasonal surface velocities generally varied by less than 10%. Assuming the discrepancy between winter and summer velocity to represent sliding velocity ( u b ) at these points, we calculate (Cuffey and Paterson, 2010, p. 310) the depth-averaged velocity ( u ¯ ): Φ = u ¯ H = u b H + 2 A n + 2 ( ρ g H α ) H 2 (2) which we find to be 85% of the surface velocity for ( A = 2.4 × 10⁻²⁴ Pa⁻³ s⁻¹) the assumed value of the flow-law coefficient for temperate ice, ( n = 3) the flow-law exponent, ( Φ ) the ice flux, (H) the ice thickness, (ρ = 910 kg m⁻³) the density of ice, (g = 9.81 m s⁻²) the acceleration due to gravity, and (α) the glacier surface slope. Without basal sliding, theoretical calculations suggest that the depth-averaged velocity is 80% of the surface velocity (Cuffey and Paterson, 2010). Thus, the ice flux through a given gate, Φ g a t e , is the sum of ice flux for each segment of the gate: Φ g a t e = ∑ i v ⊥ i H i w i (3) where v ⊥ i is the velocity component perpendicular to the ith segment, H i is the ice thickness at the ith segment and w i the width of the ith segment. We manually delineate flux gates so that they: 1) contain GPR measurements for as much of the cross-section as possible; 2) are orthogonal to the glacier centerline; and 3) represent the elevation range of the glacier. We hereafter refer to the area between flux gates as flux bins (or below the lowest gate and above the top gate). Flux in ( Φ i n ) for the uppermost bin on each glacier is assumed to be zero. We produce three sets of ice flux estimates for each year. Each set uses a different ice thickness estimate: IPR and OGGM from Pelto et al. (2020) and FAR19.

To calculate b a for each flux bin from our flux estimates, we quantify height change ( d h / d t ) and vertical ice velocity ( V z , which we use to describe both emergence and submergence velocity), while neglecting firn compaction ( V f i r n = 0.0 m a⁻¹). Each component is expressed in meters of surface displacement per year. We use these components to calculate surface height change from mass balance (h, where h = b a / ρ ), which multiplied by density (ρ) of a given bin yields the b a of a given bin (j): b a , j   =   [ d h j d t − V z , j − V f i r n , j ] ⋅ ρ j (4) We use coregistered 1 m resolution DEMs collected via repeat fixed-wing ALS surveys to generate d h / d t for each year from 2015 to 2018. ALS surveys average 1–3 laser shots m⁻² (Pelto et al., 2019). We assess height change uncertainty ( σ d h / d t ) over stable terrain after correction for effective sample size (Pelto et al., 2019), with an average σ d h / d t of 0.31 ± 0.10 m. Flux estimates (m³ a⁻¹) for each bin (Eq. 3) are normalized relative to the bin’s surface area A j (m²) to calculate V z (m a⁻¹).

While our primary analysis excludes firn compaction, we test including V f i r n in our mass-balance calculation (Eq. 4). Firn compaction is not considered in most geodetic studies because such studies primarily have multi-annual temporal scales (e.g., Tennant and Menounos, 2013; Magnússon et al., 2016; Zemp et al., 2013), and focus on glacier-wide mass balance where dynamic processes can be neglected. Even annual geodetic balance studies often neglect firn compaction (e.g., Pelto et al., 2019; Klug et al., 2018), as the required data input is onerous (e.g., Sold et al., 2013, 2015), the uncertainty high, and the net effect on mass balance small (Belart et al., 2017). We estimate V f i r n following the method of Sold et al. (2013) where over the cycle of one hydrological year and integrated over the whole firn column, one annual accumulation layer is transformed from snow to ice. For each flux bin with retained accumulation, we take the positive glaciological b a observations multiplied by the firn area of the particular bin. Mapping the area over which accumulation is continuously retained is difficult (Sold et al., 2015; Dunse et al., 2009), especially considering that firn area is shrinking on all five study glaciers, with multiple years of firn commonly exposed (Pelto et al., 2019). We merge all recent accumulation area masks (Pelto et al., 2019) from 2012 to 2018 (Figure 4) to best estimate the glacier area overlain with multi-year firn. We classify a given pixel as firn if covered by accumulation in ≥ 50% of the years. We take uncertainty on surface height change due to firn compaction ( σ V f i r n ) to be 10% (Pelto et al., 2019).

To estimate the mean density of each flux bin with retained accumulation (Figure 4), we use a firn model to estimate the density-depth relationship, which combined with the density of ice and mean ice thickness for the associated flux gate at the bottom of the bin, we estimate the mean column density of the flux bin. We use the Nabarro-Herring firn model (Arthern et al., 2010, Eq. 4), which we adapted from (Grinsted, 2021), fed with average annual temperature (Supplementary Figure 1), and average annual accumulation based on the nearest automatic snow weather station to each glacier (Supplementary Table 1) and a environmental lapse rate of −6.0 K km^-1. We scale firn column thickness by firn area for each bin (Figure 4). More details are in Supplementary Material 1.

We assess our modeled b a with our glaciological surface b a observations. The average number of annual balance measurement locations for each glacier are: Conrad, 38; Kokanee, 18; Illecillewaet, 8; Nordic, 22; and Zillmer 24 (Figure 2). See Pelto et al. (2019) for measurement details. We compare modeled flux-bin V z with V z as estimated at ablation stakes following Beedle et al. (2014) using an Eulerian framework and assuming that the measured ablation h at the stake is representative of a fixed position: V z = d h d t − h (5) and d h / d t is ALS height change sampled at the stake location. We compare modeled mass-balance gradients and flux-bin b a to observed gradients and b a . Glaciological surface b a is calculated for each flux bin by taking the mean of all observations within a given bin. If there are less than two b a observations in a given bin, we select observations proximal to the bin with a search range that is 20% larger than the bin elevation range, if necessary incrementally increasing by 20% until the requirement is satisfied. We take uncertainty on glaciological b a as 0.26 m w.e., the average for these glaciers (Pelto et al., 2019).

We use mean error ( M E ) and mean absolute error ( M A E ) to compare flux-bin b a estimates to observations: M E = 1 n ∑ i = 1 n b a , i o b s − b a , i m o d (6) M A E = 1 n ∑ i = 1 n | b a , i o b s − b a , i m o d | (7)

We assess the mass-balance profile by approximating the profile as a single linear function ( d b / d z ), but also compare using a piecewise function comprised of two linear functions above ( d b + / d z ) and below ( d b − / d z ) the equilibrium line altitude (ELA) for Conrad and Nordic glaciers. No change in slope is detected for Kokanee, Illecillewaet or Zillmer glaciers when fit with a piecewise function. The ELA is the breakpoint for our piecewise functions (Furbish and Andrews, 1984; Malone et al., 2019). Uncertainties on balance gradients are taken as the standard error ( S E ) of the slope (Rabatel et al., 2005) for linear and piecewise functions. We chose linear and piecewise linear functions because of their wide application in glaciology (e.g., Kuhn, 1984; Rea, 2009) and their use in many models (e.g., Huss and Farinotti, 2012; Clarke et al., 2013).

We estimate the ELA by regressing modeled mass-balance. We compare these ELA to satellite-observations of the transient snow line (TSL), which we sample every 50 m. We estimate the satellite-ELA as the mean ± 1σ of these TSL samples. Because of the irregular nature of the retained accumulation on our glaciers (Figure 4), the σ of the satellite-ELA can be large.

Systematic error in velocity can be assessed off-ice (e.g., Dehecq et al., 2019), and may stem from the cross-correlation itself, or DEM pre-processing. To assess systematic error, we quantify off-ice velocity, consider the uncertainty of the coregistration, and compare velocity fields against stake velocities. We assume coregistration accuracy of ± 2 pixels (Pelto et al., 2019). In the native resolution of the DEMs this equates to ± 2 m. Given that the average magnitude of off-ice velocity (<1.0 m a⁻¹) and the average difference between velocity fields and stake velocities ( ∼ 1.4 m a⁻¹) are similar, this error and the error of coregistration are summed quadratically to derive systematic error ( σ v s y s t ). Random error inherent in the cross-correlation and due to interannual variability is also estimated. Average displacement is roughly 5 pixels, and maximum displacement is 15–20 pixels in the fastest areas of flow. Thus most features do not show strong deformation and are easy to track. The accuracy of the cross-correlation is estimated to be ± 0.5 pixels, which corresponds to ± 1.5 m a⁻¹. This error, which also encompasses interannual variability, is taken as the random error. Random and systematic uncertainties combine to produce total ice velocity uncertainty ( σ v ): σ v = σ v r a n d 2 + σ v s y s t 2 (8)

While IPR observations are numerous, gaps in IPR data exist along some flux gates (Figure 3). We interpolate between the available measurements using a polynomial function (the degree of the polynomial function varies and is chosen to avoid overfitting i.e., unrealistic bed shapes) (Supplementary Figure 2), and where we interpolate, we assume an additional 10% uncertainty on ice thickness.

To calculate uncertainty in ice flux for a given gate ( σ Φ g a t e ), we propagate the uncertainties of ice velocity and ice thickness in the ice flux equation (Berthier et al., 2003) summed for all segments (i) of the gate: σ Φ g a t e 2 = ∑ i [ σ v 2 ( H i 2   +   σ H 2 ) v ⊥ i 2 ] w i 2 (9) To address uncertainty in our depth-averaged velocity (we approximate depth-averaged velocity as 85% of surface velocity), we consider full slip (depth-averaged velocity is 100% of surface velocity) and no slip (depth-averaged velocity is 80% of surface velocity) to produce a σ Φ g a t e for each velocity, and take the maximum σ Φ g a t e (Young et al., 2020). For a given flux bin σ Φ from the gate at the top of the bin represents the uncertainty on flux in ( σ Φ i n ) and σ Φ from the gate at the bottom of the bin represents uncertainty on flux out ( σ Φ o u t ). Combining uncertainties, we estimate σ Φ for the flux bin: σ Φ = σ Φ i n + σ Φ o u t (10) Uncertainty on ice column density is taken as 10%. Uncertainties on d h / d t , Φ , and V z are summed quadratically to estimate uncertainty on height change from mass balance σ h for each bin j and then to calculate uncertainty on modeled mass balance σ b a : σ b a , j = ( σ h , j ⋅ ρ j ) 2 + ( h j ⋅ σ ρ , j ) 2 (11)

The ice thickness estimates from IPR, OGGM and FAR19 produce different cross-sectional areas and variable bed shapes (Figure 3). However, the ratios of cross-sectional area between adjacent gates have a median difference of 4.8% between the three sources. These ratios play a major role in the estimates of flux in and flux out for a given bin. IPR cross-sectional areas are, on average, 6.1 ± 8.7% larger than OGGM with a range of −2.3% to +22.8%. Relative to FAR19, cross-sectional area from IPR is 16.8 ± 11.2% greater for all five glaciers. The six lowest gates at Conrad Glacier have smaller IPR cross-sections than OGGM and FAR19, and greater IPR cross-sections for the three gates in the accumulation zone. This systematic difference is also evident between the model estimates, OGGM cross-sections for the lower gates are greater than FAR19 cross-sections, and smaller than FAR19 for the higher gates (Figure 3).

Average glacier-wide ice velocity for Conrad Glacier is 18.7 ± 0.4 m a⁻¹ (Figure 5) with a maximum of 95 m a⁻¹ within the upper ice fall (Figure 2). Illecillewaet Glacier average glacier-wide ice velocity is 13.3 ± 0.5 m a⁻¹ with a maximum of 65 m a⁻¹. Average glacier-wide ice velocity is 5.0 ± 0.3 m a⁻¹ for Kokanee Glacier with a maximum of 30 m a⁻¹. Nordic Glacier average glacier-wide ice velocity is 8.0 ± 0.3 m a⁻¹ with a maximum of 48 m a⁻¹. Average glacier-wide ice velocity is 8.1 ± 0.6 m a⁻¹ for Zillmer Glacier with a maximum of 42 m a⁻¹.

The median and interquartile (IQR) range of our glacier-wide annual feature-tracking-derived velocity fields and those same velocity fields as sampled along the flux gates, are all within ± 10%. The ITS_LIVE median velocities fall near the bottom of the IQR range of our feature-tracking-derived velocities (Figure 5). Glacier-wide ITS_LIVE velocities are on average 39% slower than our estimates. For every glacier, the median surface velocity of each annual feature-tracking-derived velocity mosaic as sampled along the flux gates is typically within ± 10%. Average variability between each years’ mosaic is ± 4.0% and between each years’ mosaic and the three-year mosaic is ± 3.5%.

We find a mean difference between our 99 stake velocity pairs and feature-tracking velocity of 1.7 ± 0.51 m a⁻¹, and a median absolute difference of 3.2 m a⁻¹ (Figure 6). For Kokanee Glacier, stake velocities are less than feature-tracking velocities by 0.7 m a⁻¹. Stake velocities at Conrad, Nordic and Zillmer Glaciers are greater than feature-tracking velocities by 2.5 m a⁻¹. At Conrad Glacier, stakes close to the glacier margin along three transverse cross-sections of stakes have a larger median absolute difference (2.8 m a⁻¹) than those near the center (2.0 m a⁻¹).

Flux-bin vertical ice velocities are generally between −2 and +3 m a⁻¹ (Figure 7). Vertical ice velocities estimated at ablation stakes and from modeled ice flux are visually similar, but a quantitative comparison is not attempted. The model estimates are for entire flux bins and not intended to produce point observations. Our greatest observed stake emergence velocity was +3.2 m a⁻¹, with a surface velocity of 47 m a⁻¹. Further details are in Supplementary Material 4.

Modeled d b / d z are on average within ± 6.6% of observed d b / d z , with a median difference between modeled and observed gradients of +7.6% (Table 3; Figure 8). Thirty-six of forty-five modeled d b / d z are within ± 1.0 S E of observed d b / d z (Supplementary Tables 2–6). Six modeled d b / d z are outside ± 1.0 S E at Zillmer Glacier. All modeled d b / d z are less steep than observed d b / d z for Kokanee Glacier, where three model d b / d z outside ± 1.0 S E of observed d b / d z occur. If we instead use three-year mosaics of ice velocity rather than annual mosaics, modeled d b / d z only change by 1.6 ± 5.0% on average. For the four glaciers where it is available, using ITS_LIVE velocities decreases modeled d b / d z by 39% on average compared to our feature-tracking-derived velocities (Supplementary Figure 4). For piecewise linear functions for Conrad and Nordic glaciers, we find that the model estimates of d b − / d z were within ± 2% of observed on average, but that and d b + / d z are within ± 18% for Conrad Glacier and 45% for Nordic Glacier gradients (Figure 9). The median difference reveals that d b − / d z is three-times steeper than d b + / d z (Supplementary Tables 2–5). Variability of d b + / d z is greater than for d b − / d z for Conrad and Nordic glaciers, and the S E values of d b + / d z are on average double those of d b − / d z . Between the b a estimates from the three ice thickness sources, IPR produces d b / d z closer to observed d b / d z than OGGM or FAR19. FAR19 has the lowest d b / d z values for every glacier except Nordic Glacier. However, for the flux-bin b a estimates, FAR19 has slightly lower M E and M A E than OGGM or IPR (Table 3). Using Clarke et al. (2013) ice thickness we underestimate d b / d z at Conrad Glacier by 22.6% (Supplementary Figure 5) vs. ± 8% for modeled d b / d z from the other three sources. Regressing modeled mass balance, ELAs are on average 15 m higher than satellite-observations of ELA. Modeled ELAs are within 100 m 90% of the time. All ELA estimates outside 100 m of satellite observations occur at Zillmer Glacier in 2017, and Conrad Glacier in 2018.

Comparing observed and modeled flux-bin b a , we find an average M E of +0.03 ± 0.07 m w.e. (Table 3) and an average M A E of 0.55 ± 0.18 m w.e. If we instead use three-year mosaics of ice velocity rather than annual mosaics, M E and M A E increase slightly to 0.06 and 0.60 m w.e., respectively. Using ITS_LIVE velocities increases average M A E to 0.85 m w.e. (Supplementary Figure 4). Using ice thickness from Clarke et al. (2013) increases error with a M E and M A E of −0.42 and 0.99 m w.e., respectively (Supplementary Figure 5). Resampled to 200 m, the resolution of Clarke et al. (2013) ice thickness, M A E increases from our other ice thickness sources by 0.18 m w.e. The visual distribution of b a residuals (observed—modeled) around zero, best illustrated by the kernel density estimation in Figure 10, is encouraging. The kernel density estimation represents the probability distribution as a non-parametric estimator of density. A slight negative bias in modeled mass balance (positive b a residuals) appears over the middle elevations, but the mean residual between 0.3 and 0.8 normalized elevation where the kernel plot visualizes this bias is −0.004 ± 0.45 m w.e., with a M A E of 0.63 m w.e. There is a slight positive bias over the lower elevations with an average residual of −0.20 m w.e. ± 0.55 m w.e. ANOVA tests between the binned b a observations and each of the three sets of modeled b a estimates across the three years for each glacier have an average p-value of 0.82, with a minimum p-value of 0.58 for the Conrad Glacier IPR-based b a .

Our estimates of height change due to firn compaction range from 0.01 to 0.89 m a⁻¹ for specific bins. Including our estimate of firn compaction for flux bins which contain retained accumulation, our balance gradients are 10.6% steeper than without firnification. The inclusion of firn compaction slightly reduces error for modeled d b / d z for OGGM and FAR19, but increases d b / d z error for IPR, and greatly increases d b + / d z error for all three sources. Including firn compaction creates a positive flux-bin b a bias as indicated by M E (−0.11 ± 0.07 m w.e). M A E increases by 0.04 m w.e. to 0.59 m w.e. with firn compaction. Column-averaged-density is typically > 850 kg m⁻³ in the accumulation zones, though most glaciers have a bin with a density around 815–840 kg m⁻³. Minimum density is typically high on the glaciers where accumulation is high, and thickness relatively low.

We estimate ice velocity uncertainty ( σ v ) to be ± 2.7 m a⁻¹ (Eq. 8). Average ice velocity uncertainty ( σ v ) ranges from 15% for Conrad Glacier, to 55% for Kokanee Glacier, and averages 31%. Average ice flux uncertainty ( σ Φ , Eq. 10), ranges from 19% for Conrad Glacier to 85% for Kokanee Glacier, and averages 41%. Average modeled σ b a (Eq. 11) is 0.57 m w.e.

We consider mass to be conserved if the sum of h (glaciological balance height change) and modeled V z falls within ± 1σ of the ALS d h / d t for a given bin during a given year (Berthier and Vincent, 2012). For 84% of bins, representing 86% of total glacier area, mass is conserved (Figure 11). The three ice thickness estimates all achieved similar mass conservation on average ( ± 3%). We applied our GPS surveys to correct for mass change between data collection (Pelto et al., 2019), but fresh snow is difficult to resolve (Belart et al., 2017). In some cases, we suspect that insufficient data are responsible for mass conservation failure with either: 1) inadequate accumulation zone glaciological samples; 2) mass change between field and ALS surveys; or 3) insufficient LiDAR coverage. We further discuss these issues in Supplementary Material 3. Using field data to assess mass conservation can be problematic given the relatively sparse coverage of in situ measurements for some bins (Figure 2), and the time elapsed between ALS surveys and field observations. Ideally, mass balance would be modeled over the spatial extent of the bins (e.g., Berthier and Vincent, 2012), yet we consider this beyond the scope of our study. The relatively large terror on modeled mass balance may achieve mass conservation for some bins which have large mass-balance residuals. We thus consider our mass conservation assessment a first-order approximation.

Our modeled gradients from all three ice thickness sources accord with observed gradients (Table 3). We chose single linear functions (Figure 8), and piecewise linear functions (Figure 9) to define balance gradients (Furbish and Andrews, 1984), however, b a estimates from our method can be used to directly describe the mass-balance profile (Kuhn, 1984), or fit with any desired function. The agreement among modeled balance gradients is interesting given the differences in cross-sectional area between the three sources (Figure 3). It is unsurprising that OGGM and IPR ice thicknesses produce greater cross-sectional area overall than does FAR19 since Pelto et al. (2020) found that observed ice thickness was underestimated by previous studies (Huss and Farinotti, 2012; Clarke et al., 2013; Farinotti et al., 2019). Similar balance gradients from the three ice thickness sources is likely due to consistent ratios of cross-sectional area between adjacent gates for ice flux estimates despite different bed shapes and cross-sectional areas.

Individual flux-bin b a do not perform as well as the gradients they produce. Average M E (Table 3) is relatively low and b a residuals cluster around zero over the glacier elevation ranges (Figure 10), suggesting that the modeled estimates perform well in net, but large M A E represents the variability of modeled b a observations. We do not find a systematic bias in modeled flux-bin b a relative to observations. Large ANOVA p-values between binned b a estimates from the different ice thickness implementations suggest that the differences between the means are not statistically significant. Errors in modeled mass balance (Supplementary Material 3) may stem from: 1) errors in ice flux; 2) sparse validation data ( b a ); 3) gaps in ALS coverage; 4) differences in acquisition dates between ALS and glaciological surveys; 5) ice density used; or 6) firn compaction.

Our mass balance gradients are relatively insensitive to segmentation of the glacier (the number and location of flux gates), however, the flux bin estimates are sensitive to the segmentation. Using more gates we find that the altitude profile is more variable. Primary causes for this effect appear to be ice velocity sampling (e.g., adjacent gates on areas of faster flow and slower flow, respectively) and the effect of width scaling on cross-sectional area.

Including firn compaction steepens our mass balance gradients via an increase in accumulation zone b a . Greater slope improves some b a gradients, overall improving the fit of d b / d z for OGGM and FAR19, but degrading that for IPR and of d b + / d z for all sources. IPR d b / d z was degraded due to the relatively larger mass balance estimated for some accumulation zone flux bins, where for a number of cross-sections, greater IPR thickness relative to OGGM and FAR19 (Figure 3 H–J, N, O, T) resulted in greater ice flux. Our rationale for excluding firn compaction is to demonstrate that this method can be used in the absence of glaciological observations, and because accurately determining firn compaction is challenging (e.g., Sold et al., 2015). Our method for estimating firn compaction requires knowledge of average accumulation from field observations and the spatial pattern of accumulation (Figure 4). This undermines our ability to estimate mass-balance gradients in the accumulation zones from remote data only. We also tested using the average annual compaction indicated by the firn model to estimate firn compaction, with similar results. The firn model compaction rates were larger than our compaction estimates, and steepened the gradients by 19%, increasing error. We further discuss firn compaction in Supplementary Material 1.

ITS_LIVE velocities systematically underestimate ice velocity (Figure 5). That our feature-tracking-derived ice velocities from 3 m resolution DEMs and optical satellite imagery are greater than those from ITS_LIVE is unsurprising, particularly for smaller glaciers with a large ratio of edge to interior pixels, yet the magnitude of the difference is surprising. Ice velocity and d b / d z (Supplementary Figure 4) both being 39% smaller relative to our velocities demonstrates how ice flux scales with ice velocity in a continuity approach, and the value of accurate velocity products. These differences hold up when sampled at ITS_LIVE resolution rather than our 25 m resolution (Supplementary Figure 3). The fit of d b / d z degrades for all glaciers with ITS_LIVE velocities (Supplementary Figure 4), however, the effect is smallest for the Conrad Glacier, which at 16.9 km² is much larger than the other glaciers (Table 1). ITS_LIVE multi-year mosaics performed worse than the annual ITS_LIVE mosaics. Our glaciers feature complex terrain and are irregular in dimension, producing radial flow that may be difficult for the resolution of ITS_LIVE. We suspect that for slightly larger glaciers, or glaciers with a simple geometry, velocity products like ITS_LIVE may be readily applicable. ITS_LIVE velocities worked well for a similar approach on the large (1,096 km²) Kaskawalush Glacier (Young et al., 2020).

Worse fit of d b / d z (Supplementary Figure 5) relative to the other ice thickness sources and increased error on b a estimates when using ice thickness data from Clarke et al. (2013), is likely in part due to underestimating ice thickness for these relatively small mountain glaciers (Pelto et al., 2020), but may also be due to lower resolution. Despite degrading the fit of d b / d z , Clarke et al. (2013) ice thickness still produced gradients within ± 2.0 S E of observations.

Challenges to calculate b a via flux gates include the accuracy, coverage, and uncertainty of velocity fields, ice thickness, and surface topography. The number and representativeness of field measurements is important if assessing the accuracy of modeled b a . Our results are sensitive to the quality of our ALS surveys. The worst model performance occurred at Zillmer Glacier in 2017 (Fig. 8N), with the ALS survey occurring on November 3, 2017, 74 days after the glaciological visit (Supplementary Material 3).

The relative influence of σ H and σ v on σ Φ (Eq. 9) varies by bin and glacier. σ v dominates uncertainty for gates with ice velocity below 10 m a⁻¹. Our σ H of 10.2% for the model estimates of H implies that σ H is smaller than σ v for most gates. For regional-scale studies, σ H is typically around 30% (Huss and Farinotti, 2012; Helfricht et al., 2019). Implementing our approach for glaciers without ice thickness observations would likely imply a greater σ H than reported here. In turn, percent σ v may need to be reduced relative to ours to achieve reasonable σ Φ . This implies either improved velocity products or faster moving ice. With increased σ H , smaller, slow-moving glaciers may have a prohibitively large σ Φ for accurate flux estimates. For faster moving ice, σ v is often < 5–10% (Berthier et al., 2003; Berthier and Vincent, 2012; Burgess et al., 2013; Dehecq et al., 2019) which would decrease σ Φ substantially. As demonstrated by agreement between observed and modeled d b / d z for the five glaciers, modeled ice thicknesses from FAR19 and OGGM are of sufficient quality to represent the subglacial topography for estimates of mass balance for these glaciers.

Our velocity fields from individual image pairs reveals tolerable inter-annual variability in ice velocity but often contained data gaps (Table 2). Variability among individual velocity fields is likely due to the quality of the DEMs from variability in ALS point density, noise and artifacts in the velocity fields and data coverage, as much as to real changes in ice velocity. We chose to produce a mosaic of velocity fields for each year separately to best represent potential inter-annual variability, though find that the mean difference between the annual and three-year mosaics were minor and within uncertainties. Using a single velocity mosaic rather than annual mosaics only slightly increases M E on modeled b a . For relatively small mountain glaciers, velocity changes over short periods (e.g., three years) will typically be within velocity uncertainty, and a single ice velocity can be used to represent the period, as done by Vincent et al. (2021). Theoretically, the change in velocity due to creep from 4.95 m of thinning over the tongue of Conrad Glacier from 2016 to 2018 would yield a 2.5% decrease in velocity, well within our reported uncertainty of ice velocity. This expected slowdown due to reduced driving stress via thinning is important over longer periods (Dehecq et al., 2019). Our correction of modeled and observed ice thickness for thinning during our three-year study had a negligible effect on ice flux and b a estimates. With the current rate of thinning over the ablation zones of our glaciers, we would expect thinning to produce a change in velocity which would impact our flux-derived b a estimates after 5–10 years. Kokanee, Nordic and Zillmer glaciers demonstrate that glaciers or areas of glaciers with mean ice velocity less than 10 m a⁻¹ are challenging to measure, but possible with high resolution DEMs or imagery.

The limitation of our validation dataset is the sparse nature of our glaciological observations for some bins. Calculations of b a from glaciological observations are often biased due to uncaptured spatial variability (O’Neel et al., 2019; Zemp et al., 2013). We sampled the ablation zones of Conrad and Illecillewaet glaciers more densely than the accumulation zones (Figure 4). Pelto et al. (2019) found that the three transverse cross-sections of stakes in the ablation zone of Conrad Glacier show minor variability in b a whereas measurements in the accumulation zone show greater variability. The Fluctuations of Glaciers database also shows a bias toward ablation measurements, with a nearly 4 to 1 ratio of ablation to accumulation zone measurements (WGMS, 2018). Illecillewaet Glacier provides an example of this spatial bias where only 1–2 measurements are made above 2,500 m, representing 75% of its area. Further discussion of mass balance is in Supplementary Material 3.

Farinotti et al. (2019) (FAR19) produced globally distributed ice thickness maps and we suggest that they can be used to represent ice thickness to estimate b a from ice flux. The success of OGGM and FAR19 relative to IPR demonstrates that modeled ice thickness can adequately represent cross-sectional ice area, or at least represent the ratio of cross-sectional area for adjacent flux-gates, for ice-flux derived b a estimates for some mountain glaciers. Ice velocity estimates are becoming increasingly widespread and accurate (Heid and Kääb, 2012; Burgess et al., 2013; Dehecq et al., 2019; Gardner et al., 2020), suggesting that velocity can be reasonably determined for many glaciers. While ITS_LIVE velocities were not successful here for most years, we emphasize that our glaciers have complex topography and geometries, and are relatively small, so suggest that ITS_LIVE velocities, or similar resolution products, may perform well for glaciers with simpler geometry and in particular, for larger, faster flowing glaciers (e.g., Young et al., 2020). When extending a mass-continuity approach over accumulation zones, we recommend estimating firn column depth and density to inform total ice column density, and find this to be more important than accounting for firn compaction. For either approach, firn column depth or firn compaction must be scaled by firn area, otherwise modeled mass balance will be biased. We acknowledge that our relative success may stem from the availability of high-resolution geodetic data, which may not be available for other studies. We suspect that lower resolution data (e.g., Brun et al., 2017; Dussaillant et al., 2019) could perform well in our framework, particularly if used for larger glaciers, or if the mass-balance profile is aggregated from tens or hundreds of glaciers. If high resolution data truly are a prerequisite for success with a continuity approach, this would still increase the number of glaciers for which mass balance gradients are estimated from tens or hundreds (Rea, 2009; WGMS, 2018) to thousands or tens of thousands given the growing availability of high-resolution data (e.g., Shean et al., 2016, 2020; Belart et al., 2017; Klug et al., 2018; Pelto et al., 2019). Annual remote sensing elevation data are becoming increasingly common (e.g., Klug et al., 2018; Pelto et al., 2019; Shean et al., 2020), though most such datasets are often still nearer decadal scale (e.g., Zemp et al., 2013; Magnússon et al., 2016). While we employed a flux-gate method, if well-distributed velocity estimates are available, a gridded method (e.g., Mcnabb et al., 2012) or a full-stokes ice flow model (e.g., Jarosch, 2008) could also be used.

Recent efforts to remotely estimate surface mass balance are limited but growing (Bisset et al., 2020; Gao et al., 2020; Young et al., 2020; Vincent et al., 2021). Some of these approaches resolve point mass-balance but require data inputs which render them more suitable for individual glaciers (Young et al., 2020), or specific areas of individual glaciers (Vincent et al., 2021). Vincent et al. (2021) derive point surface mass balances from vertical ice velocities and surface elevation change. Their method estimates point surface mass balance, and demonstrates the potential for expanding the limited number of point observations of surface mass balance available globally, but relies upon an assumption of constant vertical ice velocities and requires a dense network of in situ observations, which limits application of their approach in space and time. Young et al. (2020) use a fully distributed mass-balance model driven by downscaled and bias-corrected climate-reanalysis data, requiring spatially distributed glacio-meteorological data (e.g., reanalysis products, weather station timeseries, in situ accumulation and ablation measurements). Their approach is best suited to quantify the mass budget of individual glaciers, particularly large, regionally significant glaciers. Remote-sensing-based mass-continuity approaches, like the one employed here and by Bisset et al. (2020), may be readily scaled up, and should expand the number of glaciers for which the mass balance-altitude relation can be quantified, providing model validation (of balance gradients) in a framework like that of OGGM (Maussion et al., 2019).