Multi-decadal basal slip enhancement at Saskatchewan Glacier, Canadian Rocky Mountains

Abstract

Glacier motion responds dynamically to changing meltwater inputs, but the multi-decadal response of basal sliding to climate remains poorly constrained due to its sensitivity across multiple timescales. Observational records of glacier motion provide critical benchmarks to decode processes influencing glacier dynamics, but multi-decadal records that precede satellite observation and modern warming are rare. Here we present a record of motion in the ablation zone of Saskatchewan Glacier that spans seven decades. We combine in situ and remote-sensing observations to inform a first-order glacier flow model used to estimate the relative contributions of sliding and internal deformation on dynamics. We find a significant increase in basal sliding rates between melt-seasons in the 1950s and those in the 1990s and 2010s and explore three process-based explanations for this anomalous behavior: (i) the glacier surface steepened over seven decades, maintaining flow-driving stresses despite sustained thinning; (ii) the formation of a proglacial lake after 1955 may support elevated basal water pressures; and (iii) subglacial topography may cause dynamic responses specific to Saskatchewan Glacier. Although further constraints are necessary to ascertain which processes are of greatest importance for Saskatchewan Glacier's dynamic evolution, this record provides a benchmark for studies of multi-decadal glacier dynamics.

1. Introduction

Alpine glaciers are the focus of the longest running observational records of glacier flow (Seligman, 1949; Clarke, 1987) and respond rapidly to changes in climate (Marcott and others, 2019; Marshall, 2021). Thus, they provide empirical constraints on decadal to centennial dynamic responses that remain poorly understood in other parts of the cryosphere (review in Davison and others, 2019, see their Fig. 9). On these timescales and longer, models of ice dynamics rely upon parameterized forms of processes that operate on short timescales (seconds to years), extrapolating physical relationships from observational records to much longer timescales (millennia or more) (Fowler, 2011; Larour and others, 2012). Such treatments can result in significant differences between models that explicitly include short-period processes (Helanow and others, 2021) and those that use common parameterizations (Nowicki and others, 2016; de Fleurian and others, 2018). Long-running observational records of glacier dynamics pre-dating c1984 (the launch of Landsat-5, e.g., Dehecq and others, 2019) are sparse globally. Several of these sites have remained the focus of continuous observation across several decades (Fischer and others, 2015; O'Neel and others, 2019) while others have only recently been re-established (Thomson and Copland, 2017a; Armstrong and others, 2020). These provide a growing geographic distribution of long-running datasets that may help identify multi-decadal dynamics missed in common parameterizations of glacier flow.

Glacier flow is facilitated by sliding at the ice-bed interface and internal deformation of the ice. Internal deformation is driven by deviatoric stresses acting on glacier ice, which is influenced by the surface and bed geometries of a glacier (Nye, 1952; Hooke, 2005; Adhikari and Marshall, 2012; Armstrong and others, 2016; Young and others, 2019). Rates of deformation in response to these deviatoric stresses are modulated by the viscosity of the ice, which is dictated primarily by the prevailing deformation processes, thermal structure of the glacier and concentration of interstitial water (Glen, 1955; Cuffey and Patterson, 2010; Adams and others, 2021). Basal sliding is controlled by areally integrated traction at the ice-bed interface, which is related to the distribution and magnitude of water pressures along the ice-bed interface (Flowers, 2015, and references therein) and the composition and structure of the bed itself (Lliboutry, 1968; Truffer and others, 2001; Zoet and Iverson, 2020; Helanow and others, 2021).

Surface runoff exerts a first-order control on subglacial water flow in most glacier systems, reaching the ice-bed interface through moulins and englacial conduits (Fountain and Walder, 1998; Gulley and others, 2009; Andrews and others, 2014). Consequently, different hydrodynamic regimes within a given year are commonly associated with the different phases of the surface melting cycle: early melt-season, peak melting, late melt-season and the accumulation-season (which we refer to as winter, hereafter) (Mair and others, 2002; Rada and Schoof, 2018; Davison and others, 2019). The early melt-season is characterized by accelerating and highly variable sliding velocities. This ‘spring speedup’ arises from inefficient subglacial drainage configurations (i.e., cavities and films) that rapidly pressurize in response to rising runoff supply in a sequence of ‘spring events’ that each last days to weeks (Iken, 1981; Mair and others, 2003). As the melt-season progresses, efficient drainage features form (i.e., channels), which decrease subglacial water pressures and sliding velocities, often to below annual-average rates. By the late melt-season, hydrodynamics are largely controlled by these efficient drainage system configurations (Flowers, 2015). The transition between these two regimes typically occurs near the time of peak melting and the warmest annual air temperatures (Iken and Bindschadler, 1986; Mair and others, 2002; Vincent and Moreau, 2016). As runoff supply wanes in the late melt-season, efficient drainage elements constrict, returning the subglacial drainage system to an inefficient configuration during the winter. Sliding speeds often recover toward annual-average velocities during this regime (Burgess and others, 2013; Jiskoot and others, 2017). This succession of hydrodynamic regimes can be skewed by a number of factors within a given year related to changing climatic and hydrologic conditions. In particular, inter-annual modulation of short-period runoff supply variability may impact sliding-facilitated flow rates on longer timescales via dynamics that are poorly represented in current numerical models (Schoof, 2010; de Fleurian and others, 2018).

Warming climate generally hastens the onset of the melt-season and delays the transition to winter conditions, both of which increase the number of melt days in a given year. The timing of peak melting can also shift, thereby changing the relative contributions of early and late melt-season hydrodynamics to overall glacier motion. Although rising runoff supply increases the prevalence of basal sliding on sub-annual timescales (Iken and Bindschadler, 1986; Zwally and others, 2002; Tuckett and others, 2019), increases in runoff supply in late summer can decrease sliding speeds, as observed at temperate (Vincent and Moreau, 2016), polythermal (Thomson and Copland, 2017b) and polar glaciers (Nienow and others, 2017; Williams and others, 2020). This multi-annual reduction in sliding arises from earlier formation and greater longevity of efficient subglacial drainage configurations, which lengthen the annual contributions of late melt-season hydrodynamics to annual displacement. Similarly, sustained multi-annual ice-mass losses result in glacier thinning, which correlates with decreasing ice-surface velocities for many glaciers (Dehecq and others, 2019). In contrast to this canonical view, several observational records indicate that rising multi-annual runoff supply correlates with increasing annual sliding velocities in glacier ablation zones, which is attributed to a variety of hydrologic changes. These changes include increasing intensity/frequency of spring events (Gabbud and others, 2016), more supraglacial lake drainage events (Thomson and Copland, 2017a), formation of proglacial lakes (Tsutaki and others, 2011; Pronk and others, 2021), increasing runoff supply variability (Gabbud and others, 2016; Lane and Nienow, 2019), changes in subglacial water routing and storage (Harper and others, 2007; Bartholomaus and others, 2011; Schoof and others, 2014; Rada and Schoof, 2018) or combinations thereof.

In this study, we present a newly compiled record of flow behaviors at Saskatchewan Glacier, Alberta, Canada, that spans seven decades. We characterize modern flow behaviors using new on-ice observations acquired in August 2017 and August 2019. These are combined with the initial characterization of Saskatchewan Glacier from the 1950s (Meier and others, 1954; Meier, 1957) and a number of intervening studies of ice-surface elevations and velocities. Using a first-order glacier flow model, we estimate rates of internal deformation and basal sliding velocities at two sites in its ablation zone. We consider the relative timing of our results with respect to annual melt cycles, weather, climate normals and evolving glacier geometry to motivate process-based interpretations of Saskatchewan Glacier's dynamic evolution across seven decades. Finally, we compare this system to other long-studied glaciers.

2. Field site

Saskatchewan Glacier (near 52.150°N, 117.174°W, Fig. 1) is a temperate outlet glacier of the Columbia Icefield, flowing eastward into Banff National Park of Canada (NPC) near the Banff/Jasper NPC border. In addition to ice-field contributions, the glacier currently receives ice mass from one tributary glacier and from three additional tributaries prior to the 1950s. Saskatchewan Glacier incises massively bedded carbonates of the Cathedral Formation (Ford, 1983), forming decimetric to decametric step-shaped bedforms that are exposed in the valley walls and tributary glacier forefields. These bear abundant ice-rock contact features suggesting Saskatchewan Glacier's bed is relatively till free (Fig. 2 in Iverson, 1991, location noted in our Fig. 1a as I’91). The till apron observed in the forefield likely forms in the few hundred meters up-ice from the glacier terminus where conditions are favorable for till deposition (Hart, 2006; Eyles and others, 2015), similar to other hard-bedded glaciers (MacGregor and others, 2005, and references therein) including the neighboring Castleguard Glaciers (Zoet and Iverson, 2016; Woodard and others, 2021, location Z&I’16 in Fig. 1a). In 2019, Saskatchewan Glacier was ~8 km long and 1–2 km wide following a kilometer of retreat between 1919 and 1985 and another kilometer of retreat in the past 34 years. The 1955 footprint is shown in Figure 1a for reference (Tennant and Menounos, 2013; Intsiful and Ambinakudige, 2021). This retreat is driven by sustained mass loss, which has accelerated in the past two decades (Demuth and Horne, 2017; Ednie and others, 2017; Kinnard and others, 2021). Rapid retreat left a series of recessional moraines, small proglacial lakes, braided stream complexes and a till apron behind the 1854 terminus (Eyles and others, 2015). The low-lying region behind a moraine near the 1955 terminus (Fig. 1a) coalesced into a proglacial lake that now fills the southeastern half of the forefield (compare our Figs 1c and 2 in Eyles and others, 2015). During summers 2017 and 2019, we observed that the majority of surface runoff (rain and surface melt) flowed into moulins and emerged from a discharge stream along the northwest side of the valley (Fig. 11) or discharged directly into the proglacial lake (Fig. 1c).

Fig. 1. Overview of the Saskatchewan Glacier study region. (a) Overview of the eastern half of the Columbia Icefield with glaciers labeled and the 1955 (black outline, from Tennant and Menounos, 2013) and 2019 (blue outline, also in maps b–c) extents of the Saskatchewan Glacier catchment are outlined. Unnamed tributary glaciers (T1–T3) and the Castleguard Glaciers (numbers I–IV) are labeled (referenced in text). Our study areas in the upper sector (blue box), lower sector (orange box) and transects (red lines) are marked, as are locations of glacier forefields used in comparison to our study sites (red box and blue dots; see text). (b) Upper sector seismic deployment from 2017 and sites of relevant surveys in the 1950s (Meier et al., 1954; Meier, 1957, see legend). (c) Lower sector geophysical deployment and hydrologic features from 2019 and sites from the 1950s. (d) Regional overview map with the location of Saskatchewan Glacier (star) relative to the climate reference station (CRS) at Golden, BC and the GNSS reference station at Priddis, AB (diamonds), major population centers (circles) and province borders. Basemap imagery: Orthorectified 4-band PlanetScope scene, 20 August 2019. Courtesy of Planet.com.

3. Data and methods

In this study, we characterized the modern and past dynamics of Saskatchewan Glacier using new on-ice observations co-located with past on-ice and remote-sensing studies (Meier et al., 1954; Meier, 1957; Mattar and others, 1998; Danielson and Gesch, 2011; Tennant and Menounos, 2013; and others, 2018). During August 2017 we collected GPS and passive seismic observations in the upper sector (Fig. 1b), and in August 2019 we collected GNSS, ice-surface ablation, meteorologic and seismic observations in the lower sector and glacier forefield (Fig. 1c). We supplemented these data with regional weather, climate and GNSS reference station data (Fig. 1d, see data and materials).

3.1. Past studies and regional climate

3.1.1. Meier and others: 1952–1954

During summers of 1952–1954, Meier et al. (1954); Meier (1957) conducted the first study of the structure and motion of Saskatchewan Glacier focusing on the upper sector (Fig. 1b) but also surveying on-ice sites farther up-glacier, within the lower sector, and beyond the modern terminus (Fig. 1c) (Meier et al., 1954; Meier, 1957). They determined ice-thicknesses and bed inclination at five sections of the glacier using explosive-source active seismic surveys in 1952 and proposed the subglacial topography of SG consisted of a series of glacier overdeepenings (see profiles in Fig. 1a). From 1952 to 1954, they collected extensive surface displacement and ablation measurements to characterize ice-surface velocities with one-day repeat times. During July/August 1952 and August/September 1953, they conducted sub-diurnal repeat measurements of surface flow in the Castleguard Sector, reporting 50–100% variability in (sub)diurnal surface velocities relative to mean annual velocity (Figs 13 and 15 in Meier, 1957) that correlated with air-temperatures. A sizable rainfall event between 29 and 31 August produced an acceleration of the ice-surface followed by an abrupt deceleration over the subsequent 4–12 h, likely resulting from a temporary cessation of slip as rainfall waned on the 31st (Fig. 15 in Meier, 1957). During 1952, they installed a 43 m deep borehole in the Castleguard Sector that was resurveyed in 1954, providing a strain-rate/depth relationship for the glacier necessary to estimate internal deformation velocities for the system (Fig. 32 in Meier, 1957, and in our supplementary Fig. S5). Meier et al. (1954) and Meier (1957) concluded that annual displacements of Saskatchewan Glacier were almost entirely facilitated by internal deformation between 1952 and 1954, estimating a slip velocity of less than 5 m a^-1in the upper sector (stake 6-6; Fig. 1b) and approximately zero in the lower sector (stake 14-5; Fig. 1c).

3.1.2. Supporting studies: 1948–2016

Tennant and Menounos (2013) derived digital elevation models (DEMs) of Saskatchewan Glacier from orthophotos and spaceborne imagery spanning 1948–2009 with an average repeat interval of 9 years, bounding the studies of Meier et al. (1954) and Meier (1957), and they provided source data for mass-balance analyses by Demuth and Horne (2017). Ednie and others (2017) supplemented these results with in situ mass-balance surveys, and Kinnard and others (2021) used these as model inputs to assess the relative impacts of various climate parameters on Saskatchewan Glacier's mass balance, demonstrating that mass losses accelerated in the 1990s or 2000s from values in the 1980s. Ice-surface velocities at Saskatchewan Glacier were measured using remote-sensing and in situ methods in 1995/1996 (Mattar and others, 1998) and again in 2010/2011 using remote sensing (Van Wychen and others, 2018), providing crucial ice-surface velocity estimates between the time of this study and those in the 1950s (Meier et al., 1954; Meier, 1957). Although remote-sensing data products such as GoLIVE (Fahnestock and others, 2015; Scambos and others, 2016) provide a substantial increase in the temporal coverage of glacier surface velocities globally, they do not currently provide sufficiently accurate estimates of surface velocities for narrow alpine glaciers. Ice-surface velocity estimates, ice-surface DEMs and their respective data sources are summarized in the supplementary material (Tables S1 and S2).

3.1.3. Regional climate: 1940–2020

Vincent and others (2015, 2020) produced homogenized monthly temperature records for climate reference stations (CRS) across Canada for the last century, correcting for changes in sensors and station construction, among other factors. We use records from the CRS located ~100 km south of our field site in Golden, British Columbia (Golden, BC; Fig. 1d) to represent regional temperature trends and anomalies over the past seven decades shown in Figure 2a. Annual melt-season (April–September) and winter (October–March) mean air temperature anomalies between 1940 and 2020 are summarized in Figure 2b, which were standardized to the 1981–2010 climate normal for Golden CRS (see data and materials). Selection of the melt-season/winter transitions was based upon a review of satellite imagery of Saskatchewan Glacier from the past two decades using PlanetExplorer (see data and materials). A 30-year rolling-window average of the annual mean temperatures is shown as an approximation of earlier climate (Fig. 2b). The timing of ice-surface velocity observations is highlighted for reference, and peak annual temperatures are noted in Figure 2. Between the 1950s and the present climate has warmed by 1–1.5 °C at Golden CRS with the most recent warming initiating in the 1970s. These trends are consistent with regional analyses presented in Vincent and others (2015). We note that observations in the 1990s and 2010s coincided with cooling or below normal air temperatures (Fig. 2b) and that peak annual temperatures for 2011 and 2019 were not realized until August, compared to the general trend of peak annual temperatures in July (Fig. 2a).

Fig. 2. Air temperature and climate observations and anomalies from the Golden, BC climate reference station (elevation 785 m a.s.l.). (a) Monthly mean air temperature values from 1940–2020 (black) compared to the 1981–2010 climate normal (white dashed line). Temperatures from years that coincide with ice-surface velocity observations highlighted (Vincent and others, 2020). Peak annual temperatures are marked (colored circles) and noted in the legend. The approximate timing of the winter and melt-season at Saskatchewan Glacier are shaded and labeled. (b) Annual, winter, melt-season and climate mean temperature anomalies based on the 1981–2010 climate normal for the CRS at Golden. Vertical lines show the dates of ice-surface velocity measurements and are shaded to match the legend in (a).

3.2. Estimating internal deformation and sliding velocities

Glacier motion is the result of internal deformation and basal sliding and can be written as a sum of internal (V _int) and sliding (V _slip) velocities:

(1)$$V_{{\rm surf}} = V_{{\rm int}} + V_{{\rm slip}}$$

Modeling and observational studies indicate that this summation requires a number of simplifying assumptions and/or favorable site characteristics including, but not limited to, sufficiently long observational timescales, absence of longitudinal stresses and spatially homogeneous velocity values (Hooke, 2005; Armstrong and others, 2016, and references therein). We elected to use this simple formulation in this study and acknowledge that estimates of sliding velocities using this equation are first-order approximations. Internal deformation rates in polycrystalline ice are characterized by a nonlinear viscous rheology, which has the general form (Nye, 1952; Glen, 1955):

(2)$$\dot{\epsilon}_{ij} = \left({\tau_{ij}\over B}\right)^{n}$$

with the strain rate tensor $\dot {\epsilon }_{ij}$ related to the deviatoric stress tensor τ_ij by an effective viscosity B and flow-law exponent n. We make the simplifying assumptions that deviatoric stresses are largely concentrated on planes parallel to the bed and oriented in the along-flow direction (i.e., τ_xz ≠ 0, else τ_ij = 0). This turns Eqn (1) into a scalar equation, which we solve for $\dot {\epsilon }_{xz}$. In this simplified case, the driving stress at a given depth below the ice-surface (τ_xz) is calculated as:

(3)$$\tau_{xz} = S_{{\rm f}} \rho g \zeta {\rm sin}( \alpha) $$

with the depth below the glacier surface ζ (units m b.g.s.), valley shape-factor S _f, ice-surface slope α and glacier ice density ρ_ice (= 910 kg m^-3). The shape-factor accounts for lateral drag from valley walls in its original formulation by Nye (1952), which we show here in the formulation consistent with Eqn (3) (after Hooke, 2005):

(4)$$S_{{\rm f}} = {A\over H_{{\rm i}} P}$$

with cross-sectional area A, ice-bed contact perimeter P and centerline ice thickness H _i. Thus estimates of S _f range from 0.5 (a half-pipe) to 1.0 (a sheet). In the case of S _f < 1.0, use of Eqn (1) becomes further restricted to the glacier centerline. In this case, the internal deformation velocity at the glacier centerline can be estimated by inserting Eqn (2) into Eqn (3) and integrating across the height of the ice-column (ζ ∈ [0,  H _i], Hooke 2005):

(5)$$V_{{\rm int}} = {2\over 1 + n}\left({S_{{\rm f}}\rho g{\rm sin}( \alpha) \over B}\right)^{n} H_{i}^{n + 1}$$

Numerical modeling and observation-based analyses demonstrate that valley shape influences lateral drag in two ways: through the classic hydraulic radius correction factor (Nye, 1952) and an additional slip ratio correction factor that accounts for coupling between sliding and internal deformation (Truffer and others, 2001; Adhikari and Marshall, 2012). Increased sliding rates reduce coupling stresses at the ice-bed interface, decreasing internal deformation rates, resulting in Eqn (5) over-estimating values of V _int and Eqn (1) underestimating V _slip (Adhikari and Marshall, 2012). Nonetheless, we used the first-order model described above to estimate internal deformation and sliding rates of Saskatchewan Glacier, while bearing in mind these simplifying assumptions and their impacts on the accuracy of calculated values.

3.2.1. GPS acquisition and processing

In August 2017, we collected ice-surface position observations using built-in GPS antennae of six data-loggers (DiGOS/OmniRecs DataCube³) at seismometer sites in the upper sector (Fig. 1c). These data-loggers recorded positions every 8–9 s, providing 10⁵ position estimates at each station during their 2–5 day deployment but did not report observational uncertainties, which are typically 6–12 m for individual estimates. We omitted position outliers greater than six standard deviations from raw position means and then re-calculated single estimates of position means and standard deviations from the entire 2–5 day records. These estimates had average vertical uncertainties of 1.2 m (not exceeding 5 m) and horizontal uncertainties of 1.0 m (not exceeding 3 m).

In August 2019, we collected ice-surface position observations with a set of three Emlid RS+ single-band (L1) real-time kinematic (RTK) GNSS antennae in the lower sector. We maintained a GNSS reference station consisting of one antenna mounted on a tripod in the glacier forefield (BASE, Fig. 1c), and we used other antennae as rover units mounted on surveyor poles. We surveyed 72 campaign sites in the lower sector between 1 and 19August, marked with 1 m long bamboo wands installed in half meter deep boreholes that also served as ablation stakes. From 6 and 19 August, we installed and continuously operated a rover antenna at site ROV1 (Fig. 1c) by inserting its survey rod mount into a 0.8 m deep borehole. To prevent melt-out, we drilled new boreholes every 2–4 days, shifting the site by a few meters. We recorded all GNSS data at a 5 Hz sampling rate in RINEX format and campaign measurements were collected with a 1 min acquisition period.

We post-processed raw GNSS reference station (BASE, Fig. 1c) data using a static carrier-phase double-differencing routine implemented in the open-source RTKLIB package (see data and materials). To minimize effects of atmospheric delay and multi-path signals, we omitted data from satellites below 20° during all post-processing. For our initial processing step, we used reference data from the Canadian Active Control System (CACS) continuously operating GNSS station PRDS located ~240 km away in Priddis, Alberta (Priddis, AB; Fig. 1d). We then calculated the absolute position of the base station (BASE; Fig. 1c) as the mean of 1.19 × 10⁷ processed positions. Using data from BASE as reference, we determined the kinematic GNSS positions of the rover antennae via kinematic carrier-phase double-differencing in the RTKLIB package, correcting for antenna-heights in the process. Our processing configuration options are documented in the repository and additional details of the double-differencing routine can be found in the RTKLIB manual Appendix E7. Due to the short acquisition times of our survey, we only used high-precision fixed-integer-ambiguity resolution positions for velocity calculations, which had average horizontal uncertainties of 5 mm (not exceeding 16 mm) and average vertical uncertainties of 13 mm (not exceeding 45 mm). Using only ”fixed” status data gave rise to some data gaps in addition to two longer data-gaps from temporary power failures at BASE.

3.2.2. Ice-surface velocities

To estimate ice-surface velocities in the lower sector, we applied and expanded upon established methods (Hoffman and others, 2011; Mejia and others, 2021). We first rotated post-processed data and their covariance matrices into an along-flow (N44°E)/across-flow (N46°W) basis. Next, we corrected continuous data at ROV1 for displacements due to re-installation of ROV1 and abrupt, sustained shifts of the antenna due to small melt-out effects by estimating the average ice-surface velocity during re-installations and removing residual displacements (see supplementary material). We then removed impulsive signals from wind-gusts using an iterative, rolling-window standard score (z-score) detection algorithm akin to those used for earthquake detection to automatically identify and remove these record segments that exceeded a score of three (Withers and others, 1998). Removed segments rarely exceeded 1 min. Continuous displacement processing steps are illustrated in the supplementary material (Fig. S1).

We then estimated average ice-surface positions and ice-surface velocities with a centered moving-window application of a generalized weighted least squares (WLS) linear fitting to displacement data in the along-flow and across-flow directions, using inverse data variances as weights. Similar to Mejia and others (2021), we used a 4 h sampling window with a 1 s step size to focus on (sub)diurnal velocity trends in the continuous record. We restricted parameter estimation to windows with fewer than 25% null values (gaps) to suppress edge-effects from large gaps and the ends of our record. By using the rolling-window inversion scheme on post-processed observations we simultaneously inverted for window-average positions, window-average velocitiesand model covariance matrices for both parameters. We then calculated long-term-average velocities using the same WLS scheme from campaign GPS measurements at survey sites with at least three observations and from the entire observational record at ROV1 to provide a comparison to continuous velocity estimates at ROV1. Horizontal velocity uncertainties estimated with the 4 h rolling window WLS did not exceed 0.4 m a^-1 (1.5 cm d^-1, see Fig. S1e in the supplementary material).

3.2.3. Ice-surface profiles and feature extraction: 1948–2019

We extracted smoothed ice-surface elevation profiles and slopes along transects A–A’, B–B’ and C–C’ (Fig. 1a) from ice-surface DEMs (Tennant and Menounos, 2013) and long-term-average ice-surface positions in 2017 and 2019. We do not use data from Danielson and Gesch (2011) for ice-surface elevations due to the large uncertainties (see Table S2), but we do use their data for additional constraints on elevations of recently deglaciated bedrock surfaces. Raw data profiles were generated by projecting data points within 100 m of A–A’ onto the transect using the QGIS TerrainProfile plugin (QGIS Association, 2021; Jurgiel and others, 2021). We applied the same moving-window WLS linear fitting scheme as ice-surface velocities to estimate smoothed surface elevation profiles, slopes and parameter covariance matrices, using a 600 m long window with a 10 m step size and requiring a minimum of three data-points per window to estimate parameters. We used data-source-specific uncertainties from Tennant and Menounos (2013) (see Table S2 in the supplementary material) as the standard deviation of extracted DEM values and assumed a 2 m vertical uncertainty for their reference year: 1986. Similar to A–A’, we extracted elevation data within 100 m of B–B’ and C–C’ and projected them into the cross-section planes to estimate across-flow transect surface elevation profiles. Unlike A–A’, we used third-order WLS polynomial fittings to projected data for each transect, creating surface models and associated covariance matrices. We used these polynomial representations to inform a Bayesian analysis of internal deformation velocities detailed in a later section and in the supplementary material.

3.2.4. Valley geometry

We estimated ice-thicknesses in both sectors using the horizontal to vertical spectral ratio (HV) technique on ambient seismic recordings at geophone sites and velocity estimates from our active-source survey in 2019 (Figs 1b–c). We combined these with ice-surface elevation data from published DEMs (Danielson and Gesch, 2011; Tennant and Menounos, 2013) and our GPS observations. The HV technique estimates harmonic frequencies of compliant materials from ambient vibrations. The fundamental frequency (f ₀) in sufficiently simple glacier settings is representative of the local ice-thickness (see Preiswerk and others, 2018a, for details) and is related to ice thickness (H) by the shear-velocity of glacier ice (V _S) (Bard and Bouchon, 1980a, 1980b):

(6)$$H = {V_{{\rm S}}\over 4 f_0}$$

We provide a detailed description of the seismic deployments in 2017 and 2019 and the HV workflow in the supplementary material (Fig. S3), which follows methods from prior applications of HV to recover glacier thicknesses (Picotti and others, 2017a; Preiswerk and others, 2018a; Yan and others, 2018a; Kohler and others, 2019; Stevens and James, 2022). Similarly, details on active-seismic data analysis to estimate V _S are included in the supplementary material (Fig. S4). Using seismic estimates of ice-thicknesses and ice-surface elevation data from 2017/2019, we created a DEM for the ice-bed interface in both sectors.

Glacier erosion over the course of seven decades is likely small (less than 1 m, e.g., Hallet and others, 1996) compared to uncertainties from our seismic analyses (10¹ m), so we used the same valley-shape model to estimate parameters in Eqn (4) and vary only the ice-surface elevation model to approximate valley-fill geometries. Following methods from Hilley and others (2020), we approximated the valley-shape in each sector by fitting a third-order polynomial WLS to ice-bed interface elevations from 2017/2019 and recently exposed bedrock elevations estimated from DEM data. We used the combination of the ice-surface and ice-bed interface polynomial models to estimate parameters in Eqn (4) and interpolate ice-thickness and ice-surface elevation data in the lower sector using a radial basis function interpolation scheme to create a DEM of the ice-bed interface.

3.2.5. Internal deformation velocities

We estimated an effective viscosity (B) for Saskatchewan Glacier by inverting borehole strain-rate data from Meier (1957) and depth-dependent driving stresses calculated using Eqn (2) as the kernel for a generalized WLS inversion. We used reported observation depths as true vertical depths and the surface slope reported at the borehole site in 1952 (α = 4.41°; p. 59 in Meier, 1957) to calculate driving stresses with Eqn (3). We found that using the reported S _f from 1952 (S _f = 0.62; Meier, 1957) produced unrealistically low values of B compared to those expected for temperate ice (Cuffey and Patterson, 2010), so we assumed S _f = 1 in our calculation of driving stresses (i.e., negating effects of lateral drag) based on how shallow the borehole was (ζ_max = 43 m) compared to the total thickness of the ice in 1952 (H _i = 401 m, Meier et al., 1954). Meier (1957) reported strain-rates at nine levels in the borehole so we conducted a leave-one-out-cross-validation analysis (LOOCV) to assess sensitivity of estimated B values to individual data points. Finally, we calculated internal deformation velocities in both sectors for each available DEM. We quantified V _int uncertainties by propagating observational uncertainties through Eqns (2)–(5) with Monte Carlo Markov Chain (MCMC) simulations. Simulations for transect/DEM pairs consisted of 1000 realizations of perturbed parameter values informed by parameter covariance matrices to estimate the posterior distribution of V _int based on our estimates of ice-surface geometries, bed geometries, surface slopes and ice effective viscosity but other limitations associated with a potential reduction in coupling resulting from slip in the shape factor assessment were not included. Further details on the MCMC simulations are included in the supplementary material.

3.2.6. Sliding velocities

Ice-surface DEMs and ice-surface velocity observations aligned only in 2019 (see Tables 1 and 2), requiring interpolation of V _int estimates to estimate sliding velocities via Eqn (1). We used a linear interpolation scheme to estimate the interpolated V _int values at the mid-point date of associated ice-surface velocity observations (see Table S1) and analytically propagate uncertainties, assuming no covariance. In actuality, internal deformation and sliding velocities covary, but we lacked sufficient data to quantify this relationship (e.g., Adhikari and Marshall, 2012).We tested the statistical significance of sliding velocity estimates using a two-sample z-test:

(7)$$z-score = {\bar{V}_{{\rm surf}} - \bar{V}_{{\rm int}}\over \sqrt{\sigma_{V_{{\rm surf}}}^{2} + \sigma_{V_{{\rm int}}}^{2}}}$$

from the mean values of ice-surface ($\bar {V}_{{\rm surf}}$) and internal deformation ($\bar {V}_{{\rm int}}$) velocities and their variances $\sigma _{V_{{\rm surf}}}^{2}$ and $\sigma _{V_{{\rm int}}}^{2}$, respectively. We used this statistic to test a null hypothesis of no-sliding at the p = 0.97 (z-score = 2) confidence bound, rejecting the null hypothesis for z-scores higher than two; that is, an estimated sliding velocity using Eqn (1) was statistically significant if its z-score exceeded two (Townsend, 2002).

Table 1. Estimates of geometric parameters and associated driving stresses estimated from the upper and lower sector transects from the 1948 DEM

These are the reference values for percent differences shown in Figures 6a–b.

Table 2. Summary of estimated internal deformation and sliding velocities in the upper and lower sector of Saskatchewan Glacier, and associated z-scores (see text)

Dates correspond to acquisition times of surface velocity estimates in Table S1.

3.3. Runoff supply modeling

We extracted insights on inputs to the subglacial hydrologic system in the lower sector during August 2019 by modeling runoff supply (the sum of rates of surface-melting and rainfall). To estimate a time-series of runoff supply to the lower sector during August 2019, we used the enhanced temperature index (ETI) model of Pellicciotti and others (2005) in its incoming-radiation formulation (Bouamri and others, 2018, their Eqn 8):

(8)$$\bar{M} = \left\{\matrix{ f_{\rm T} \bar{T} + f_{{\rm I}} \bar{I},\; \hfill & { \bar{T} > T^{\ast }}\hfill \cr 0,\; \hfill & {\bar{T} \le T^{\ast }} \hfill }\right.$$

This model estimates the average melt production rate $\bar {M}$ for a period of time using observations of average air temperature ($\bar {T}$) and incoming shortwave radiation ($\bar {I}$), empirically derived coefficients f _T (temperature coefficient) and f _I (incoming shortwave radiation coefficient) and a threshold temperature T* (= 0°C, generally). We collected temperature, incoming shortwave radiation and rainfall observations every 30 min with an Ambient Weather WS-2000 automated weather station (AWS) installed in the forefield (Fig. 1a, elevation 1783 m a.s.l.) ~200 m lower than sites in the lower sector. To estimate melt rates, we collected repeat ice-surface ablation measurements at all survey locations with measurement accuracy of 5 mm and calculated pairwise estimates of $\bar {M}$ and synchronous estimates of $\bar {T}$ and $\bar {I}$ at each site, producing two to five estimates per station with sampling periods of 2–19 days for a total of 433 data points. We then used these estimates of $\bar {M}$, $\bar {T}$, $\bar {I}$ to invert for f _T and f _I with a generalized least-squares using Eqn (8) as the kernel. We propagated data uncertainties through this inversion by applying an 10 000 sample MCMC simulation in which data were randomly perturbed by their uncertainties and model parameters were estimated from the resultant posterior distribution. We then used Eqn (8) to conduct the forward problem, calculating hourly estimates of melt rate (M) from hourly mean AWS observations (T and I) and model parameter estimates (f _T and f _I). Finally, we summed melt rate estimates with synchronous rainfall rates to calculate a time-series of runoff supply for the lower sector during August 2019.

4. Results

We begin by presenting observations of surface location solutions, velocities, ablation, local weather and glacier geometry in the ablation zone of Saskatchewan Glacier from our surveys in 2017 and 2019. Next, we present extracted geometric, flow and rheologic parameters from previous studies to estimate internal deformation velocities across seven decades. Finally, we present a statistical assessment of sliding velocities based on surface velocity observations. Measurement and parameter uncertainties are presented as a single standard deviation unless otherwise stated.

4.1. Observations: 2017/2019

We show AWS data used for our ETI model in Figure 3a, which display typical diurnal cycles, and rainfall rates in Figure 3b, which are sporadic and rarely exceed 1 mm w.e. h⁻¹. Our inversion for ETI model coefficients resulted in estimates of f _T = 5.65 ± 0.11 mm w.e. d⁻¹ °C⁻¹ and f_I = 0.0085 ± 1.4 × 10⁻⁵$f_{\rm I}$ mm w.e. m² d⁻¹ W⁻¹. These values are comparable to values derived in studies of other alpine systems: f _T ∈ [3,  6]  mm w.e. d⁻¹ °C⁻¹, f _T ∈ [0,  0.12] mm w.e. m² d⁻¹ W⁻¹ (Pellicciotti and others, 2005; Bouamri and others, 2018). Our surface runoff rate model is shown in Figure 3b along summed surface runoff rates (rainfall and melting). We calculated daily cumulative runoff to serve as a proxy for relative subglacial water flux, setting the start of the hydrologic day at 4:00 local time when modeled M reached minimum daily values. There was a power failure for the AWS on 7 August resulting in a data gap during which time a thunderstorm occurred and we account for this gap in our daily cumulative runoff estimates by using the average M from available values. Runoff estimates are consistent with our observations of surface ablation and rainfall played a minor role in total surface runoff rate trends during our August 2019. Daily cumulative runoff estimates show three multi-day cycles with multi-day plateaus around 65 mm w.e. d⁻¹ separated by minima on the 3rd, the 12th and the 17th when rates dropped below 50 mm w.e. d⁻¹.

Fig. 3. Continuous observations and models of meteorologic data, surface runoff and ice-surface velocities in the lower sector in August 2019. (a) Air temperatures (T) and incoming shortwave radiation flux (I) measurements from AWS. (b) Runoff supply rates from melting (maroon), rainfall (blue), their sum (black) and their uncertainties (translucent envelopes). Daily cumulative melting rates are overlain (bars) with days with partial data coverage noted. (c) Ice-surface velocities (V _surf), their 24 h rolling-average ($\bar {V}_{{\rm surf}}$), calculated internal deformation velocities (V _int) and inferred sliding velocities ($\hat {V}_{{\rm slip}}$, red shading) at ROV1 are presented in m a⁻¹ and cm d⁻¹. Uncertainties of V _int are visible, but those for V _surf are not due to their small values (see text).

The average ice-surface velocity at ROV1 calculated from its entire record (3.88 × 10⁶ samples) was 63.1 ± 0.001 m a⁻¹, which is moderately higher than values from campaign GNSS measurements at adjacent sites with comparable acquisition periods (7 or more days), but within standard errors ([37–53] ± [1–15]m a⁻¹, see Fig. S2 in the supplementary material). This indicates that the stitching procedure to correct instrument moves was reasonably successful.

Continuous ice-surface velocities at ROV1 are shown in Figure 3c with annual (m a⁻¹) and daily (cm d⁻¹) rate units for ease of reference. We included these uncertainties in Figure 3c, but they barely exceed the width of the plotted line. Additionally, we show the 24 h rolling-average of the surface velocity record ($\bar {V}_{{\rm surf}}$) for comparison with daily cumulative runoff trends in Figure 3b. We observed prominent diurnal cycles in the surface velocity record, with maximum values regularly exceeding 100 m a⁻¹ and minimum values falling below 20 m a⁻¹. Daily mean velocities ranged from 50–70 m a⁻¹, which is consistent with the overall average velocity at ROV1 and neighboring campaign velocity measurements (Fig. S2). Velocity variability was higher before 9 August compared to subsequent observations. This is also reflected in the daily mean velocities. We interpret these observations as periods of more-responsive (before the 9th) and less-responsive (after the 9th) basal hydrodynamic conditions to comparable rates of runoff supply to the bed (compare Figs 3b–c).

4.2. Glacier geometry: 1948–2019

The elevation profiles we extracted from DEMs and our GPS/GNSS observations for along- and across-flow transects are shown in Figure 4. We display the smoothed surface elevation profiles estimated simultaneously with surface slopes on transect A–A’ (Fig. 4a) and elevation profiles in across-flow transects approximated as first-order polynomials are overlain on transects B–B’ and C–C’ (Figs 4b–c). Surface elevations we identified as recently exposed bedrock from data in Tennant and Menounos (2013) and Danielson and Gesch (2011) are shown in Figures 4b–c, which we used with ice-bed interface elevations HV analyses (see supplementary material) to model the valley elevation profiles using a 3rd order polynomial fit (Figs 4b–c, maroon line). These valley elevation profiles are comparable to those estimated by Meier (1957) and HV estimates of the ice-bed interface elevations are not significantly different from the active-source seismic estimates of Meier et al. (1954) (see Figs 4a–c). Ice thickness uncertainty estimates were larger for the upper sector as a consequence of lower f ₀ values associated with thicker ice (see Eqn (6)).

Fig. 4. Ice-surface, bed and exposed bedrock elevation data (points) and modeled surfaces (lines) for the (a) along-flow (A–A’), (b) upper sector (B–B’) and (c) lower sector (C–C’) cross-sections. Ice-surface elevation data and models are colored by year and transect intersection locations are marked in red. Data and model fits for the valley profile are shown in maroon and data from (Meier, 1957) are denoted as M57. Plot notations are communal and documented in the legend. DEM years correspond to Table S2.

Figure 5a shows our model of subglacial topography in the lower sector using the entire seismic array. We observed several topographic features that are consistent with the scale and geometry of those observed in neighboring forefields that incise the Cathedral Formation (see Figs 1a and 5b). Specifically, we observe a central depression bounded up-flow by a steep slope and down-flow by flow-perpendicular ridge that we interpreted as an overdeepening bounded by a headwall and a bedrock riegel (Patton and others, 2016). There is a second depression on the down-flow end of the seismic array that we also interpret as the upper end of a subsequent overdeepening. The riegel has an along-flow notch with an approximate relief of 6 m, and it is situated directly down flow from two active moulins observed during 2019. These factors make this a likely position for a subglacial channel (Dow and others, 2011; Banwell and others, 2016). These interpretations are illustrated in Figure 5. The 2017 survey lacked sufficient array density to map subglacial topographic features in the upper sector.

Fig. 5. Interpreted subglacial topographic features in the Cathedral Formation from (a) HV and GPS surveying in the lower sector of Saskatchewan Glacier and (b) satellite imagery of the Castleguard I forefield (locations in Fig. 1a, frame colors identical). Topographic feature interpretations are labeled in red (OD = overdeepening) and hydrologic features are marked in blue (see legend and text). Locations of geophones, ROV1 and transects (red lines) are shown for reference. Basemap imagery in (a) is the same as Figure 1 and imagery in (b) was sourced from GoogleEarth Pro.

Estimates of the ice-filled valley cross-sectional area (A), ice-bed interface contact perimeter (P), centerline ice-thickness (H _i), valley shape factor (S _f), surface slopes along A–A’ (α) and driving stresses (τ_xz) for each ice-surface DEM in both sectors are summarized in Figure 6 as percent differences relative to estimates from 1948 (see Table 1). We found that the S _f has not appreciably changed, with an increase of 3% in the upper sector and 10% in the lower sector between 1948 and 2017/2019. The ice-surface has monotonically lowered in both sectors (Figs 4a and 6) thinning the glacier by 27 and 52% in the upper and lower sectors, respectively. The rate of thinning has roughly doubled in the lower sector in the last decade compared to the preceding six decades (Fig. 6b). Ice-surface slopes progressively steepened, with a 31% increase in the upper sector and a 63% increase in the lower sector during this record. Despite these changes to surface geometries, driving stresses remained relatively constant, declining <5% in the upper sector and ~18% in the lower sector. This indicates that effects of thinning on driving stresses were largely compensated by the effects of steepening ice-surfaces.

Fig. 6. Percent changes in flow cross-section geometric parameters and driving stress for the (a) upper sector and (b) lower sector relative to modern estimates (see Table 1). Scaled standard deviations are shown as dotted envelopes.

4.3. Flow dynamics: 1948–2019

Our inversion of borehole strain-rate data returned estimates of B = 0.19 ± 0.01 MPa a^1/3, with a slightly lower estimate of B (0.18 ± 0.01 MPa a^1/3) when individual data below 15 m b.g.s. were omitted (see Fig. S5 in the supplementary material). Although our estimate of B is derived from a shallow borehole where strain-rates are low, the inversion produced comparable estimates to average values of B for temperate ice (0.21, MPa a^1/3 Cuffey and Patterson, 2010) and the values recovered from full ice-thickness borehole deformation observations at the neighboring Athabasca Glacier (location in Fig. 1a) using similar analyses as those presented here (B = 0.15 MPa a^1/3, Raymond, 1971) and accounting for higher-order flow mechanics (B = 0.23 MPa a^1/3, Adhikari and Marshall, 2012). We therefore take the ensemble average from the LOOCV analysis (B = 0.19 ± 0.02 MPa a^1/3) as a reasonable parameter estimate and use this value to calculate internal deformation velocities shown in Figure 7. We then interpolated estimates of internal deformation velocities ($\hat {V}_{{\rm int}}$) to coincide with ice-surface velocity measurements (summarized in Table 2) along with estimates of sliding velocities ($\hat {V}_{{\rm slip}}$) and the two-sample z-scores based on the distributions of V _surf and $\hat {V}_{{\rm int}}$.

Fig. 7. Ice-surface velocities, internal deformation velocities (purple) and their uncertainties between 1948 and 2019 in (a) the upper sector (transect B–B’) and (b) the lower sector (transect C–C’). Surface velocities are colored by their sampling time (see legend). Uncertainties are shown at one (box) and two (whisker) standard deviations, and mean values are marked by thin white lines. Ice-surface velocity measurements with uncertainties smaller than 1 m a⁻¹ are marked with squares rather than box-plots (see Table S1).

We find that surface velocities are statistically indistinguishable from internal deformation velocities in the lower sector during the 1950s (Table 2) and estimate that internal deformation accounts for the majority of flow in the upper sector during this time. These results are consistent with the internal-deformation dominated mode of flow characterized in Meier's analyses (Meier, 1957). Sliding velocity estimates in the 1990s and 2019 are statistically significant, returning z-scores greater than two in most cases. During the 1990s and 2010s our mean estimates suggest that sliding accounted for roughly half the observed surface velocities in both sectors, and 60–80% of observed surface motions in the lower sector during the 2010s (Table 2). Estimates for the 2011 sliding velocity had larger uncertainties arising from surface velocity uncertainties in (Van Wychen and others, 2018, see Table S1) and increased surface roughness from the 2009 DEM (Fig. 3a and Table S2) that increased surface slope uncertainties (Table 1 and Fig. 6), and by extension, internal deformation velocities (Fig. 7 and Table 2). Consequently, these estimates of sliding velocity failed to refute the null hypothesis at the z-score = 2 confidence bound in both sectors. However, we note that the average values for this estimate are comparable to estimates from the 1990s and 2019 that successfully rejected the no-sliding null hypothesis.

5. Discussion

Observations of Saskatchewan Glacier spanning decades indicate that glacial motion in the late melt-season has undergone a transition from a regime dominated by internal deformation to one that prominently features basal sliding. Detailed velocity estimates from 2019 support the link between water inputs to the subglacial drainage system and basal sliding, indicating subglacial hydrology likely plays a role in this transition. The transition occurred without a significant change in driving stresses – with effects of ice thinning counterbalanced by ice-surface steepening. Whereas velocity observations in the winter months are fewer in number, they also indicate an increased relevance of sliding between the 1950s and 2000s. The trend of increasing sliding velocities in response to warming climate deviates in some ways from observations of glaciers elsewhere (Vincent and Moreau, 2016; Davison and others, 2019). The development of a proglacial lake (Fig. 1c) and the subglacial geomorphology of Saskatchewan Glacier may offer explanations for this anomalous behavior. Lakes tend to support elevated basal water pressures near glacier termini, and their formation at a glacial toe can cause a sustained enhancement of glacier sliding near the terminus (e.g., Tsutaki and others, 2011; Pronk and others, 2021). Additionally, the decimetric to decametric bedrock steps that form from subglacial erosion of the Cathedral Formation may promote more widespread and more stable conditions for the subglacial cavity network formation compared to other locations (see Helanow and others, 2021; Woodard and others, 2021). This could enhance the connectivity of these classically inefficient drainage features compared to other glacier settings, increasing the area over which water pressures can rapidly respond to changes in surface runoff to the bed (Murray and Clarke, 1995; Schoof and others, 2014; Rada and Schoof, 2018). The shape of basal bedrock morphology specifically affects the drag response provided by the glacier bed (Schoof, 2005; Zoet and Iverson, 2015, 2016; Zoet and others, 2021). Though many glacier beds display a variety of basal morphological features (Woodard and others, 2021) that could cause the glacier to encounter a wide range of adverse slopes upon cavitation, the Cathedral Formation erodes into an end-member, stepped-bed morphology with a relatively limited range of obstacle shapes (see Woodard and others, 2021), which may make it distinctively sensitive to changes in the subglacial effective pressure.

5.1. Lower sector dynamics: August 2019

Our estimate of the internal deformation velocity for 2019 indicates that the vast majority of observed surface motion over the study period is the result of sliding in the lower sector with a high degree of statistical confidence (red filled area in Fig. 3c; values in Table 2). Though this interpretation hinges upon a number of simplifying assumptions detailed in our methods section, potential correction factors to the flow-law model used here to better replicate higher-order modeling results (Adhikari and Marshall, 2012) would decrease our estimates of V _int further, particularly when sliding velocities are high and ice-bed coupling is weak. Furthermore, the surface velocity minimum on 11 August (6.1 ± 0.3 m a⁻¹) is not significantly different (z-score = 1.44) from our estimate of V _int (8.9 ± 1.7 m a⁻¹, Table 2). We interpret this as a short-lived cessation of sliding similar to the event observed by Meier (1957) in 1953. Following the framework of Adhikari and Marshall (2012), this could be a period of high ice-bed coupling and higher V _int values along the glacier centerline that may represent a maximum value for V _int during August 2019. We conclude that our estimates of sliding velocities are reasonable first-order approximations and that glacier motion during August 2019 was dominated by sliding.

We observed coherent diurnal and multi-day cycles between ROV1 velocities and runoff supply rates during August 2019 indicating a linkage between variations in subglacial hydrology and sliding velocities. Similar coherence of diurnal cycles in runoff and surface velocities is observed near subglacial channels where water pressures quickly respond to changes in runoff inputs (Murray and Clarke, 1995; Mair and others, 2003). ROV1 is located down-flow of two active moulins and a notch in the subglacial topography, which would tend to favor channel formation (Figs 1c and 4, see Fountain and Walder, 1998; Banwell and others, 2016). Sensitivity of flow velocities near ROV1 to runoff inputs to the bed was likely further compounded by the presence of a bedrock riegel directly up-flow of the site, which also tends to increase the sensitivity of basal water pressures to changes in subglacial water fluxes (Hooke and Pohjola, 1994; Dow and others, 2014). Thus, changes in the daily average and variability of surface velocities likely reflect changes in subglacial water fluxes and drainage system architecture that modulate basal water pressures.

Diurnal cycles are present in both time-series but there is an abrupt decline in the amplitude of daily velocity variations around 9 August that initiates a multi-day decline in mean daily velocities (Fig. 3c). This diurnal amplitude transition occurred during a period of elevated, but relatively stable, runoff supply. We interpret this as growth of the subglacial drainage system, resulting in lower basal water pressures on the 9th compared to the 8th, and thus slower sliding velocities. The multi-day decline in mean daily velocities coincides with a fall in daily runoff supply through 13 August (Figs 3b–c) and continued to covary through the rest of our observational record. However, the amplitude of diurnal velocity variations did not recover to pre-9 August values. We interpret this as a relative expansion of the subglacial drainage system that suppressed pressurization of the bed in late August compared to early August, similar to hydrodynamics characterized in other alpine settings (Gimbert and others, 2016; Nanni and others, 2020).

August temperatures are largely past peak annual temperatures for this region, which generally occur in July. However, 2019 was a notably cooler melt-season (Figs 2a–b), which may have extended the early melt-season compared to adjacent years and caused our observations to coincide with the transition between early and late melt-season dynamics. This is consistent with our observations of declining diurnal velocity variability following 9 August (Fig. 3c). It may be the case that the transition on 9 August represents a shift from early to late melt-season hydrodynamics. In this case, average velocities before the 9th may represent near-maximum velocities for 2019, while those after the 9th may more closely represent the annual average velocity if Saskatchewan Glacier generally follows the classic hydrodynamic behaviors of other glacier systems dominated by runoff supply (Davison and others, 2019, and references therein). Sub-seasonal variability in runoff supply and dynamic adjustment of the subglacial drainage system can also trigger similar patterns (Gimbert and others, 2016), but without further context for the 2019 melt-season and constraints on subglacial drainage system architecture, these process interpretations can only be circumstantially supported.

5.2. Dynamics across seven decades

The suite of ice-surface velocity measurements in this study rarely exceeded a sampling period of one month (Table 1), necessitating their interpretation within the context of relative climate, weather and melt-season dynamics (e.g., Mair and others, 2002; and others, 2018, and references therein). The majority of melt-season observations were collected after annual peak temperatures (compare Table 1 and Fig. 3). Although these melt-season measurements provide an imperfect estimate of annual flow behaviors, the sparser winter surface velocity measurements (Table 1) can help provide further context for annual velocities, which may be as little as 10% below annual values depending on the relative intensity of the prior melt-season (Van Wychen and others, 2018, and references therein). For example, Shackelton Glacier – an comparably scaled outlet glacier of the neighboring Clemenceau Icefield – has annual centerline velocities comparable to winter velocities for Saskatchewan Glacier in the 2010s (compare Fig. 3 in Jiskoot and others, 2017, with our Fig. 7).

Surface velocity observations in the 1950s were statistically indistinguishable from estimates of internal deformation velocities in the lower sector (z-scores less than two, Table 2) and of significance in the upper sector (z-scores greater than two, Table 2) in agreement with Meier's conclusions that the flow of Saskatchewan Glacier was dominated by internal deformation (see Supplement 7 in Meier, 1957). In the 1990s–2010s, differences between surface velocity observations and internal deformation velocity estimates become statistically significant except for the 2011 winter observations in both sectors and the November 1995 observation in the lower sector (Table 2). Diminished z-scores arose from larger uncertainties in both surface velocities and driving stresses associated with changes in observational techniques from land-based surveying to remote sensing (compare Figs 3a and 7, quantified in Table S2). Overall our observations suggest a multi-decadal enhancement of basal sliding between successive melt-seasons, and perhaps between successive winters. While a number of hydrologic factors may contribute to this enhancement, the evolution of the glacier's surface geometry has also likely made the system more sensitive to changes in the hydrologic system.

5.3. Interpreting dynamics

Saskatchewan Glacier's modern-day dynamics provide insights into the evolution of its dynamics over the last seven decades. Driving stresses at both study sites remained relatively constant between 1948 and 2019 due to steepening surface slopes that kept pace with declining ice-thickness (see Figs 6a–b). Thinning ice reduces maximum overburden stresses on the ice-bed interface, which has multiple impacts on hydrodynamics at the bed. First, decreased ice overburden pressures are more easily counterbalanced by increases in basal water pressure, thereby requiring smaller variations in water supply to the bed to meaningfully influence sliding velocities (Iken and Bindschadler, 1986). Second, smaller overburden stresses decrease the creep-closure rate of subglacial cavities and channels, supporting the longevity of well-connected sections of the subglacial drainage system (Schoof and others, 2014; Rada and Schoof, 2018). Under ongoing warming (Fig. 2a), these combined effects of thinning on subglacial drainage system architecture may play a meaningful role in the observed enhancement of sliding at Saskatchewan Glacier. Further, the observations in the 1990s and 2010s were collected during/following colder winters, which are shown to promote faster sliding velocities during subsequent melt-seasons (Sole and others, 2013; Williams and others, 2020).

Hydrodynamics associated with the formation of the proglacial lake may also help to explain enhanced sliding velocities during melt-seasons in the 1990s and 2010s compared to the 1950s, similar to the enhancement of flow characterized at Rhonegletscher, Switzerland (Tsutaki and others, 2011) and in the Himalayan region (Pronk and others, 2021). Hydrologic and flow monitoring of Haut Glacier d'Arolla also shows an enhancement of sliding coincident with modern warming (Gabbud and others, 2016) and increased variability in daily subglacial discharge (Lane and Nienow, 2019). Increasing subglacial water flux variability on diurnal timescales is shown to enhance sliding compared to steady-state conditions (Schoof, 2010), which may help explain Haut Glacier d'Arolla's multi-decadal acceleration. Similar processes may be affecting Saskatchewan Glacier.

In addition to hydrologic facilitation of enhanced sliding, the geometry of Saskatchewan Glacier's bed may also contribute to the enhanced sliding we infer. Recent numerical analyses of sliding and subglacial cavity dynamics over realistic beds based on the neighboring Castleguard IV glacier forefield (Fig. 1 in Zoet and Iverson, 2016, approximate location in Fig. 1a) indicates that step-shaped bedforms promote widespread cavity formation compared to other bed morphologies (Helanow and others, 2021). The scale of these steps and cavities are well correlated (Lliboutry, 1968; Cohen and others, 2006; Andrews and others, 2014; Zoet and Iverson, 2016). The Castleguard Glaciers II–IV incise the Pika/Elodin Formations, which have decimetric to metric bedding that erodes into similarly scaled bedrock steps (Ford, 1983; Woodard and others, 2021). In comparison, the decimetric to decametric bedding and bedforms produced by glacier erosion of the Cathedral Formation should result in larger subglacial cavities than those modeled on representations of the Castleguard Glacier IV forefield (Helanow and others, 2021). Larger cavities increase the likelihood of developing a well-connected distributed drainage system that can rapidly respond to changes in surface runoff supply to the bed and potentially pressurize wider regions of the bed compared to systems with smaller bedrock obstacles (Harper and others, 2007). However, larger drainage networks also tend to have greater storage capacity, requiring greater fluxes of water to trigger pressure changes (Bartholomaus and others, 2008, 2011). Glacier thinning (and thus reduction of overburden pressure) would counteract some of this effect, requiring a smaller absolute basal water pressure increase to appreciably change the stress state of the ice-bed interface. Coupled with the inter-annual acceleration of mass loss at Saskatchewan Glacier (Kinnard and others, 2021), the associated increase in meltwater availability would tend to increase sliding-facilitated displacement, particularly if meltwater supply variability on short timescales increases as well, similar to the observed behaviors at Haut Glacier d'Arolla (Gabbud and others, 2016; Lane and Nienow, 2019).

Our first-order estimates of processes facilitating the motion of Saskatchewan Glacier provide a new multi-decadal record of glacier flow. They demonstrate a significant enhancement of basal sliding during the melt-season and perhaps in the winter months between the 1950s and the 1990s/2010s. Several factors may explain this result, whether alone or in concert with each other. However, further constraints are necessary to better ascertain the processes responsible for enhanced sliding at Saskatchewan Glacier. This work provides a somewhat anomalous case-study compared to most other glacier systems with long-running dynamics records. Further geophysical, remote-sensing and numerical investigations are warranted to better constrain potential ice-dynamic behaviors not well-characterized in commonly used models of glacier flow.

6. Conclusions

We revisited an early, comprehensive study of glacier dynamics at Saskatchewan Glacier and connected it to modern observations to create a new seven-decade-long history of flow for this hard-bedded, temperate, valley glacier. Our high-resolution records of surface velocities and surface runoff indicated a tight linkage between variations in water input to the glacier bed and sliding velocities in its ablation zone. We produced a higher resolution map of the ice-bed interface in the ablation zone, resolving a series of glacier overdeepenings that provided supporting evidence for this interpretation. We then used ice-surface DEMs spanning the 1940s–2010s and borehole deformation data from the 1950s to parameterize a first-order glacier flow model and estimate rates of internal deformation across this seven-decade record. We compared these with co-located estimates of ice-surface velocities across seven decades and found a statistically significant enhancement of sliding-facilitated motion across with sliding accounting for roughly one-third of observed surface velocities in the 1950s while sliding accounted for the majority of observed surface velocities in the 1990s/2010s. We found our results from the 1950s to be in good agreement with the results in (Meier, 1957), further supporting our interpretations.

We hypothesize that the tight linkage between runoff supply and sliding velocities characterized in 2019 is related to the multi-decadal enhancement of sliding-facilitated motion. Several hydrologic processes observed in other glacier settings provide possible explanations – such as the effects of thinning ice overburden on subglacial drainage system architectures or the formation of a proglacial lake. The presence of massive, step-shaped subglacial obstacles also likely contributes to dynamics and architecture of the subglacial drainage system at Saskatchewan Glacier, which is sensitive to variations in water inputs to the bed. Though many constraints are still necessary to ascertain the processes driving enhanced sliding at Saskatchewan Glacier, this multi-decadal record provides a benchmark for future studies analyzing the effects of changing climate and glacier geometry on overall glacier dynamics across similar timescales.