Reduced basal motion responsible for 50 years of declining ice velocities on Athabasca Glacier

Abstract

The time-evolution of glacier basal motion remains poorly constrained, despite its importance in understanding the response of glaciers to climate warming. Athabasca Glacier provides an ideal site for observing changes in basal motion over long timescales. Studies from the 1960s provide an in situ baseline dataset constraining ice deformation and basal motion. We use two complementary numerical flow models to investigate changes along a well-studied transverse profile and throughout a larger study area. A cross-sectional flow model allows us to calculate transverse englacial velocity fields to simulate modern and historical conditions. We subsequently use a 3-D numerical ice flow model, Icepack, to estimate changes in basal friction by inverting known surface velocities. Our results reproduce observed velocities well using standard values for flow parameters. They show that basal motion declined significantly (30–40%) and this constitutes the majority (50–80%) of the observed decrease in surface velocities. At the same time, basal resistive stress has remained nearly constant and now balances a much larger fraction of the driving stress. The decline in basal motion over multiple decades of climate warming could serve as a stabilizing feedback mechanism, slowing ice transport to lower elevations, and therefore moderating future mass loss rates.

1. Introduction

The response of glacial systems to warming climate conditions over the past several decades and into the future has numerous, significant societal implications. Ice mass loss from alpine glaciers and the Greenland and Antarctic ice sheets will control the rate of sea level rise over the course of the 21st century (Zemp and others, 2019; Rounce and others, 2023). The decline and loss of glaciers also has numerous local impacts, such as profound changes in water flow, sediment transport, ecosystem dynamics and slope stability (Clague and Shugar, 2023). Understanding current and estimating future rates of ice mass loss require constraints on several physical processes which control the flow dynamics of these systems. The viscous deformation of glacial ice has long been a topic of research and is generally assumed to be well described by the empirical Glen flow law relating stress and strain rate in polycrystalline ice (Nye, 1952; Hooke, 1981). Debris entrainment in the several meters of ice above the glacier bed can impact ice deformation rates in this transition zone (Rempel and others, 2023). The basal motion of a glacial system combines the effects of deformation within the subglacial till and sliding at the ice-bed interface (Weertman, 1957; Truffer and others, 2001). In general, basal drag decreases with higher slip velocity or increased water pressure for idealized hard-bedded glaciers with water-filled cavities (Helanow and others, 2021). However, when applied to observed bedrock topographies, process-based models show that this rate-weakening drag effect may become negligible (Helanow and others, 2021). This finding suggests that a universal slip law may exist which correctly predicts basal motion over both hard and soft-bedded glaciers (Zoet and others, 2020; Helanow and others, 2021). For many temperate glaciers, the processes contributing to basal motion constitute most of the total ice velocity (Raymond, 1971; Amundson and others, 2006; Vincent and Moreau, 2016). Therefore, constraining parameters which control basal sliding is critically important to understanding flow dynamics (Raymond, 1971; Amundson and others, 2006; Adhikari and Marshall, 2012). Basal water pressure is known to be a primary control of basal motion because water can decouple ice from the underlying bedrock and increase sliding speeds (Iken and Bindschadler, 1986). The difference between subglacial water pressure and ice overburden pressure determines the net effective water pressure of subglacial drainage networks (Röthlisberger, 1972; Schoof, 2010) and it is generally accepted that the effective pressure is a main control of basal velocity regardless of the nature of the bed (Iken and Bindschadler, 1986; Clarke, 1987).

Basal motion has the potential to serve as either a positive or negative feedback to ice mass loss in a warming climate. Some studies have suggested that enhanced meltwater generation from warmer temperatures could result in increased basal motion, resulting in a positive feedback mechanism that would greatly accelerate ice mass loss in Greenland (Parizek and Alley, 2004). Recent work suggests that basal motion slows following periods of enhanced melt, potentially due to the development of an efficient subglacial drainage network which increases the effective pressure (Truffer and others, 2005; Schoof, 2010; Tedstone and others, 2015; van de Wal and others, 2015; Stevens and others, 2016). Dense borehole observations documenting a complex variety of channelized, distributed and isolated portions of the subglacial drainage network at any given time complicate the conceptual model of a drainage network which gradually evolves over the course of a melt season to accommodate higher meltwater input (Rada and Schoof, 2018, 2023). Numerical modeling of subglacial drainage evolution in Greenland also shows that poorly connected regions of the bed as opposed to increased channelization may better explain late summer and fall changes in ice flow dynamics (Hoffman and others, 2016). Studies in Greenland have also shown that both summer melt and wintertime supraglacial lake drainages can cause ice accelerations due to enhanced sliding associated with these transient meltwater inputs (Maier and others, 2022, 2023). Other work in Greenland has documented a slowdown in ice velocity during a period of high melt from 2000 to 2012 followed by ice acceleration from 2013 to 2019 caused by a ~15% reduction in surface melt (Williams and others, 2020). These observations suggest basal motion is sensitive to changes in meltwater input on annual to decadal timescales. For other regions such as High Mountain Asia, long-term decadal changes in ice velocity are largely explained by changes in ice thickness and driving stress, suggesting these systems are less sensitive to changes in basal conditions (Dehecq and others, 2019). Other work has observed declining surface velocities on land-terminating mountain glaciers throughout the world due to negative climatic surface mass balance (Heid and Kääb, 2012). However, studies regarding the spatial and temporal variation of basal motion in response to water input are generally limited to relatively short timescales ranging from diurnal to decadal (Burgess and others, 2013; van de Wal and others, 2015). The direction and magnitude of basal motion change on multi-decadal to centennial timescales is critical to understanding whether a stabilizing or destabilizing feedback will occur because of climate warming. Enhanced basal motion would lead to higher rates of mass transport to lower elevations, where the ice can melt or calve more quickly. Conversely, reduced basal motion would retain ice at higher elevations, where it is subject to lower rates of melt. Sliding parameters derived from short-term observations of water pressure and basal velocity may not predict long-term sliding behavior well (Iken and Truffer, 1997). The lack of research observing and describing how basal motion evolves on decadal to centennial timescales provides the direct motivation for this work.

Athabasca Glacier, located in the Canadian Rockies, provides an ideal study site for investigating the long-term evolution of glacial basal motion (Fig. 1a). Past glaciological research on Athabasca Glacier provides the necessary baseline datasets to make our current study possible. In the 1960s, multiple borehole studies (Paterson and Savage, 1963; Raymond, 1971) measured both the rates of internal ice deformation and simultaneous surface velocities, allowing for historical estimates of basal motion. Raymond (1971) calculated deformation and basal motion over an annual period using inclinometer data from multiple boreholes along transects transverse to the ice flow direction. These historical datasets critical to our current study have been digitized and more recently republished (Armstrong and others, 2022).

Figure 1. (a) Athabasca Glacier model domains used for the historical (red) and modern (orange) time periods in the finite element model Icepack. The locations of the Paterson stake surface velocity measurements and the Raymond boreholes are shown. Background satellite image is from Worldview-2 imagery (© Maxar 2018). (b) Cross-sectional model domain of Athabasca Glacier transverse profile coincident with the Raymond (1971) borehole transect ‘Section A’ marked by the green line in panel (a). Ice thicknesses along this transect are constrained from boreholes drilled to the bed in 1966.

Here, we use two modeling frameworks to show how basal motion has evolved from the 1960s to the present day due to changes in ice geometry and meltwater production modifying the driving stress and basal friction fields. First, we implement a 2-D ice flow model to compute englacial velocities through a transverse cross section of the Athabasca Glacier valley. We run and tune the model using known values from the historical datasets and then perform contemporary model runs using a modern-day geometry and surface slope. Next, we perform a broader analysis using a 3-D numerical model, Icepack (Shapero and others, 2021), to gain insight as to how the basal friction, basal velocity and stress distribution have changed throughout the lower Athabasca Glacier over the past 50 years. Due to its simplicity, the cross-sectional model allows us to run the forward model many times and estimate the ice rheology. The 3-D finite element model allows for a broader, glacier-wide view of the parameters of interest. Both models result in similar conclusions and corroborate the findings from a much simpler model assessment in Armstrong and others (2022).

2. Methods

2.1. Study area

Athabasca Glacier (52.19° N, 117.26° W) is an outlet glacier of the Columbia Icefield located in Jasper National Park, Alberta, Canada. It is an ionic glacier for tourism in the Canadian Rockies with access via tour buses via an ice road seen in Figure 1a. The glacier flows from ~3400 m above sea level through three large icefalls. Athabasca Glacier has a total area of ~15 km² and a length of ~9 km (RGI 7.0 code: 12982) from its ice divide to the terminus (Ednie and others, 2017; Randolph Glacier Consortium, 2023). It has a large accumulation to total area ratio of ~73% (Ednie and others, 2017), but has experienced significant rates of terminus retreat (20 m a⁻¹ from 1919 to 1948) and surface elevation change (−0.91 m water equivalent (w.e.) a⁻¹ from 1919 to 2009) (Tennant and Menounos, 2013). Our study area consists of the lower ablation area of the glacier below the lowest elevation icefall and is ~3 km long and 1 km wide (Fig. 1a). We model ice flow through the same valley cross section studied extensively with borehole inclinometry measurements from July 1966 to August 1967 (Raymond, 1971; Armstrong and others, 2022). We use a standard coordinate system (x-axis down-glacier, y-axis transverse to flow and z-axis with ice depth) and produce a cross-sectional transect in a yz-plane, coincident with the locations of five historical boreholes or ‘Section A’ from Raymond (1971; Fig. 1b). We create two meshes based upon historical and modern ice thickness measurements at this valley cross section using Gmsh (Geuzaine and Remacle, 2009). The historical mesh seen in Figure 1b is based upon ice thickness observations from these boreholes and is interpolated between known point measurements. This results in several abrupt corners in the bed perimeter which are unlikely to be real, but the influence on ice flow will be minor and localized. Deformation data from full-depth boreholes 2, 3, 4 and 5 are used in our analysis.

We define two map-view model domains corresponding to a modern (orange) and historical (red) Athabasca Glacier extent (Fig. 1a). These geometries are uploaded into the numerical ice flow model, Icepack for use in a diagnostic, first-order hybrid model which resolves both sliding and deformation velocities (Shapero and others, 2021; Section 2.2.3). We have chosen the lateral extent of our model domains based upon where we have a high degree of certainty regarding ice thickness and bed topography from a recent radar survey of the glacier below the icefall (Armstrong and others, 2022). This constraint excludes several historical velocity point measurements on the debris-covered ice near the lateral moraines.

2.2. Cross-sectional modeling

We delineate an outline of the bed perimeter and produce a finite element mesh over this domain. To calculate the out-of-plane (i.e. down-glacier or x-direction) flow velocities (u) through this cross section, we simplify the Stokes equations by assuming the transverse and vertical velocities (v and w), and all out-of-plane velocity gradients (i.e. du/dx) are zero. This is strictly correct only for a top surface with no longitudinal gradients or cross-glacier slope and no in-plane components of basal motion, but well approximates the first-order controls on ice flow (Truffer and others, 2001; Amundson and others, 2006; Armstrong and others, 2016). This reduces the Stokes equations to the Poisson equation, relating the gravitational forcing to the ice viscosity and spatial velocity gradients (Eqn (1)). In this equation ice velocity (u) is a function of ice viscosity (η) and the force of gravity depends upon the ice density (ρ), gravitational acceleration (g), the ice thickness (H) and the surface slope (z _s). The 2-D Nabla operator ($\nabla$) indicates derivatives in the in-plane coordinates.

(1)$$\nabla \;\cdot ( ( {\eta ( u ) \nabla u} ) = \rho gH\nabla ( {z_{\rm s}} ) $$

Ice viscosity is calculated according to Eqn (2), using the rate factor (A), strain rates and the flow law exponent (n).

(2)$$\eta ( u ) = \displaystyle{1 \over 2}A^{{-}1/n}\left({\displaystyle{1 \over 4}{\left({\displaystyle{{\partial u} \over {\partial x}} + \displaystyle{{\partial u} \over {\partial y}}} \right)}^2} \right)^{( 1-n) /2n}$$

The velocity-dependent viscosity of the ice is then set to an initial value and the system of equations is iteratively solved (Piccard iteration) until convergence (Amundson and others, 2006; Armstrong and others, 2016). The boundary conditions used include a stress-free surface and prescribed basal velocities defined along 58 segments of the bed perimeter. For the historical model run, we prescribe the observed basal velocity profile as the boundary condition (Raymond, 1971). We next determine the values of the rate factor (A) and flow law exponent (n) which minimize the misfit between observed borehole and modeled ice deformation.

2.3. Ice rheology tuning

The ice rheology parameters which properly describe the viscous ice deformation according to Glen's flow law have been a subject of numerous studies on a large range of spatial scales (Glen, 1955; Goldsby and Kohlstedt, 2001; Millstein and others, 2022). From Glen's initial laboratory experiments to recent work on the topic, constraining the stress exponent and rate factor which describe ice flow behavior in this power–law relationship is a critical part of modeling ice flow (Behn and others, 2021).

Here we investigate the parameter space of the flow exponent, n, and rate factor, A, which enter the system of equations in the viscosity model we describe above. All model runs were solved using the same mesh with a historical ice thickness estimate and using a constant, laterally averaged surface slope value of 3.4° as reported at the time (Raymond, 1971). We perform a grid search, running the model with all possible combinations of these parameters over a plausible range (A = 2.3 × 10⁻²⁴ to 3.6 × 10⁻²⁴ Pa⁻ⁿ s⁻¹; n = 2.6–4.0). After the model converges for a given set of parameters, we quantify misfit between the modeled englacial velocities and the observed borehole velocities from Raymond (1971). We extract a ‘synthetic borehole’ from the model output at each historical borehole location to estimate the depth-dependent down-glacier velocity profile u(z). We then calculate root mean square error (RMSE) between each modeled deformation curve and the observed borehole deformation digitized from the Raymond study at 1 m intervals. Cumulative RMSE values for each parameter combination suggest the optimal values of n = 3 and A = 2.70 × 10⁻²⁴ Pa⁻³ s⁻¹. Each of the four boreholes’ misfit is weighted equally in the cumulative RMSE.

2.4. Modern day model targets

After determining the most plausible ice rheology parameters (A and n) using the known historical ice thickness, surface slope and englacial velocities, we proceeded to run the cross-sectional model to ascertain modern rates of basal motion. We used a modern ice thickness estimate and surface slope of the modern-day ice surface averaged across the valley width (~1100 m) (Armstrong and others, 2022). A modern borehole campaign to measure contemporary ice deformation rates on the Athabasca Glacier is ongoing, but the data have not yet been analyzed to constrain modern ice deformation. The lack of modern, observed ice deformation rates means that it is not possible to have an englacial velocity model target as was done with the historical, observed data in Section 2.3. Instead, we run the ice viscosity model with the new model domain under the assumption that flow law parameters have not changed.

Surface velocity measurements of the Athabasca Glacier were recorded from March to September 2020 using an array of high precision GPS stations, three of which are coincident with the Raymond borehole transect (Armstrong and others, 2022). Velocities calculated from this field campaign are likely biased high relative to the annual average because they spanned the winter/spring speed up event and fast summer motion but do not include winter data, when ice flow is typically slower. The data were processed and utilized to estimate an annual average surface velocity as described in Armstrong and others (2022). This contemporary surface velocity dataset for the Athabasca Glacier is used as our model target for annual surface velocity. We assume that the flow law parameters and spatial distribution of basal motion are unchanging in time and infer the modern values of basal velocity by iteratively rescaling the historical basal motion field to produce agreement with observed modern day surface velocity in the cross-sectional ice flow model.

2.5. Numerical modeling in icepack

While the cross-sectional model allowed us to quickly tune model parameters and estimate modern basal velocities along a single transverse profile, the numerical ice flow model Icepack allows us to investigate the evolution of basal motion across our entire study reach. We use Icepack to run 3-D model inversions of Athabasca Glacier to infer the basal friction field. Icepack is a finite element model built upon the partial differential equation solvers and functionality of Firedrake (Shapero and others, 2021). The model domains used for this study are shown in Figure 1b as the red and orange boundaries for the 1966 and 2018 glacier geometry derived from orthorectified imagery (Armstrong and others, 2022). These model domains each have four perimeter boundaries: the glacier terminus, up-glacier icefall and two sidewalls. We define the up-glacier edge and terminus as Dirichlet boundaries by specifying the ice surface velocity of the system. For the modern day, these velocity boundaries are prescribed from observed Landsat-derived velocities (Fig. S1; Armstrong and others, 2022). For the historical model runs, the velocity boundaries are 140% of these modern day values to approximate the observed point velocity measurements from the 1960s (Fig. S1). We define the sides of the domain as ‘sidewall boundaries’ or a mixed Dirichlet–Robin boundary (Shapero and others, 2021). These sidewall boundary conditions constrain the system by ensuring no normal ice flow and accounting for resistive stresses exerted by the valley walls using a lateral sliding law (Shapero and others, 2021). The input datasets seen in Supplementary Figures S1–S4 are interpolated onto the finite element mesh and used to initialize the velocity boundary conditions, surface elevation, ice thickness and bed elevation of the model domain for each time (Tennant and Menounos, 2013; Farinotti and others, 2019; Armstrong and others, 2022). We use the modern-day surface elevation (Tennant and Menounos, 2013; Menounos and others, 2019) and a modern ice thickness (Armstrong and others, 2022) to estimate the bed elevation (Fig. S2). Assuming this bed elevation is constant through time, we derive a historical ice thickness estimate for the entire model domain by subtracting the bed elevation from a historical estimate of surface elevation derived from 1966 aerial photos (Armstrong and others, 2022; Fig. S3). The input data are then extruded into a single layer 3-D mesh using high-degree basis functions and Legendre polynomials using functionality built into Icepack. We use the ‘HybridModel’ in Icepack which solves for both plug (‘shallow shelf’) and shear (‘shallow ice’) flow as would be expected for a valley glacier such as the Athabasca Glacier with components of both vertical shear and basal sliding (Shapero and others, 2021).

After setting up the model mesh, data initialization and boundary conditions, we then set up the friction inversion, which is iteratively solved using a conjugate gradient method with a Tikhonov-style regularization (Calvetti and others, 2000; Konovalov, 2012; Shapero and others, 2021). Icepack parameterizes bed friction as a Weertman-style law (Weertman, 1957) according to Eqn (3), where w is the exponent in the nonlinear friction law and basal shear stress, τ _b, is a function of the sliding velocity, u _b, and the friction, C (units: ${{{\rm MPa}\cdot {\rm y}{\rm r}^{( w/( {\rm w}-1) ) }} \over {{\rm m}^{( w/( {\rm w}-1) ) }}}$)

(3)$$\tau _{\rm b}( {u_{\rm b}, \;C} ) = -C\vert {u_{\rm b}} \vert ^{{1 \over {( {w-1} ) }}}u_{\rm b}$$

We introduce a new variable, θ, to ensure that the friction coefficient remains positive and physically meaningful during the inversion (Shapero and others, 2021). We define θ as the natural logarithm of the basal friction coefficient (C), which becomes the variable we seek to recover during the inversion using measured surface velocities for both the modern and historical time periods. The forward model then calculates a surface velocity field based upon the updated basal friction field until the solution converges by minimizing a misfit function which is the sum of both a loss functional (penalizing misfit between observed and modeled surface velocities) and a regularization functional (penalizing data overfitting through physically unrealistic ‘rough’ solutions). For the modern surface velocities, we implement Eqn (4) as the loss function which integrates the misfit (E) between observed (u _obs) and modeled (u) surface velocities weighted with an observational error of σ. We use Landsat-derived surface velocities interpolated onto the mesh of the model domain as the observed velocity dataset for our study area, u _obs. The integral is evaluated across the entire plan-view model domain, Ω, for Eqns (4) and (6).

(4)$$E( u ) = \displaystyle{1 \over 2}\mathop \smallint \limits_\Omega ^\;\left({\displaystyle{{u-u_{{\rm obs}}} \over \sigma }} \right)\;^2{\rm d}x$$

For the historical surface velocities, measurements are only available at a discrete number (k = 40) of points across the glacier's surface (Paterson, 1962) and thus the loss functional is implemented as a summation across these limited points as in Eqn (5).

(5)$$E( u ) = \mathop \sum \limits_k \displaystyle{{{\vert {u( {x_k} ) -u_{{\rm obs}}} \vert }^2} \over {2\sigma _k^2 }}^\;$$

The regularization functional (R) implemented is the same when inverting for both the modern and historical basal friction fields as seen in Eqn (6):

(6)$$R( \theta ) = \displaystyle{{\alpha ^2} \over 2}\mathop \smallint \limits_\Omega ^\;\vert {\nabla \theta } \vert \;^2$$

As with other ill-posed inverse problems, minimizing the total misfit becomes a balance of minimizing the data-model misfit without introducing unrealistic, large spatial gradients into the basal friction solution. The regularization term contains an L ² norm of the gradient of the friction coefficient which is needed for the inversion to be stable. We utilize an L-curve approach to determine the regularization factor, α, which optimizes the tradeoff between data fit and solution smoothness as has been used in similar inversions (Habermann and others, 2013). Decreasing α increases the importance of model-data misfit to the solution at the expense of increasing the model norm and decreasing the smoothness of the solution. For both the historical and modern glacier geometry, we conduct a series of friction inversions and plot the norm of the regularization function versus the model-data misfit for the converged solution (Fig. 2). We chose regularization factors of α = 100 and α = 200 for the modern and historical geometries, respectively. These values are located near the corner of the L-curve where data misfit becomes only marginally better with increasing solution roughness (arrows in Fig. 2). All results presented for basal friction, velocity fields and stress distribution use these values in the regularization function.

Figure 2. Tikhonov regularization L-curves for friction inversions conducted with modern ice parameters (a) and historical parameters (b). We find the tradeoff between model-data misfit and the roughness of the solution is optimized at approximately α = 100 and α = 200 (denoted by arrows) for the regularization function. Model misfit is higher when comparing model output to a Landsat-derived raster in modern inversions than between model output and historical point velocity measurements.

2.6. Force balance

We perform a 2-D force balance calculation to estimate the change in the gravitational driving stress (τ _dx) and the resistive terms of longitudinal stress gradient (F _lon), lateral drag (F _lat) and basal shear stress (τ _bx). We use a methodology which calculates the terms of the force balance along the flow direction (van der Veen and Whillans, 1989; O'Neel and others, 2005):

(7)$$\tau _{{\rm d}x} = \tau _{{\rm b}x} + \displaystyle{\partial \over {\partial x}}( {H{\bar{R}}_{xx}} ) + \displaystyle{\partial \over {\partial y}}( {H{\bar{R}}_{xy}} ) \;$$

We calculate τ _dx according to a standard expression for gravitational driving stress as seen in Eqn (1). The final two terms on the right-hand side of Eqn (7) respectively represent the longitudinal stress gradients (F _lon) and the lateral drag (F _lat). These are calculated using the ice thickness (H) and the $\bar{R}_{xx}$ and $\bar{R}_{xy}$ components of the depth-averaged membrane stresses. These stresses (Eqn 6 and 7 in O'Neel and others, 2005) can be directly evaluated in Icepack. We do not follow the sign convention used in O'Neel and define all stresses on the right-hand-side of Eqn (7) to be positive if they resist flow. We assume that the bridging stresses $\bar{R}_{zz}$ are negligible (O'Neel and others, 2005). We align the x-coordinate along flow and the y-coordinate transverse to flow by rotating the resistive stress components ($\bar{R}_{xx}$ and $\bar{R}_{xy}$) using appropriate tensor rotation. We calculate the τ _dx, F _lat and F _lon components of the force balance at a gridcell resolution of 400 m and determine τ _bx as the residual of these three terms for both the modern-day and historical Icepack velocity solutions. We perform these calculations for five overlapping 400 m by 400 m centerline grid cells with centers translated 200 m down glacier from the previous cell. The center of the upper-most gridcell corresponds to zero along the flow-aligned x-profile. We calculate the average τ _dx within a given gridcell. We determine the average F _lat along the sides of the cell along the valley margins and the average F _lon along the upstream and downstream sides of the cell using a finite difference approach.

3. Results

3.1. Ice rheology

Model-observation misfit was minimized (RMSE = 3.51 m a⁻¹) using rheology values of approximately n = 3.0 for the Glen flow law exponent and A = 2.70 × 10⁻²⁴ Pa⁻³ s⁻¹ for its pre-factor, although there is a considerable co-dependency of n and A (Fig. 3). We utilize these rheology values for all following analyses for both the cross-sectional and 3-D model for both the historical and modern periods.

Figure 3. Cumulative RMSE between modeled and observed englacial velocities calculated for each permutation of 15 different power law exponents (n) and 10 rate factor values (A). The contour plot shown in A vs n parameter space calculates the total misfit from four full-depth boreholes (Fig. 4). The optimal ice rheology values occur at the red X, with n = 3 and A = 2.70 × 10⁻²⁴ Pa⁻³ s⁻¹.

Figure 4. Modeled (blue) and observed (orange) depth-dependent velocity profiles, u(z), for four borehole sites from Raymond (1971). Modeled profiles were calculated using the ice rheology pair which minimized the cumulative misfit between observed and modeled ice deformation (n = 3, A = 2.70 × 10⁻²⁴ Pa⁻³ s⁻¹) across all boreholes (Fig. 3).

Model-observation misfit using these globally optimized parameters varies substantially between boreholes (Fig. 4). Each panel plots the observed englacial velocities from digitized borehole data (blue) and the modeled velocity which minimized the global RMS error between these velocity profiles (Fig. 4). Using the ice rheology parameters which minimize cumulative misfit does not necessarily produce the best fit for each borehole individually.

3.2. Valley cross-sectional model

After determining the most suitable ice rheology parameters, we use these values with a modern ice geometry. The modern-day mesh uses an updated ice thickness estimate and a steeper modern surface slope of approximately 3.8° across the valley profile (Armstrong and others, 2022). We run the model using this modern-day geometry at the same location studied by Raymond (Fig. 1a). Prescribed basal velocities ~30% lower than the historical value (Fig. 5b), or a decline of ~12 m a⁻¹ at the centerline, are required for the modeled surface velocity to match the observed GPS-derived surface velocities (Fig. 5a). Comparison of the modern GPS-derived surface velocities and historical surface velocities reported by Raymond (1971) show that the surface velocity has decreased by ~15 m a⁻¹ at the centerline (Armstrong and others, 2022). We therefore find the decline in basal motion constitutes ~80% of the total observed decline in surface velocities over this 50-year period. We also find that the transverse shape of the basal velocity magnitude (a landmark result of the Raymond (1971) study) remains similar for the modern time.

Figure 5. (a) Modeled historical (blue) and modern (orange) transverse surface velocity profiles with observed 1967 (pink circles) and observed 2020 GPS surface velocities (green circles). These results show that prescribing basal sliding at 70% of the historical values (30% decline) in a modern ice geometry allows the model to reproduce observed 2020 GPS surface velocities. (b) Two basal velocity profiles, 100% (blue) and 70% (orange) of their historical values as reported by (Raymond, 1971). X = 0 corresponds to the western margin.

3.3. Icepack plan view hybrid flow model

Athabasca Glacier shows a significant increase in the inverted basal friction coefficient (C) throughout the model domain (median increase = 127%; Fig. 6), which corroborates the prescribed decline in basal sliding observed in the cross-sectional flow model (Section 3.2). This increase in basal friction is evident throughout the study area and spans the full width of the glacier (Fig. 6). These results show that much of the modern glacier is underlain by a ‘stickier’ bed than in historical times. Due to the diffusive nature of ice flow (Balise and Raymond, 1985; Gudmundsson, 2003), the spatial pattern of surface velocity fields for the converged model solutions has remained similar (Fig. 7). The total decline in surface velocities (Fig. 7a) can be partitioned into contributions from a modest decline in ice deformation (Fig. 7b) and a larger decline in basal motion (Fig. 7c). The percent of total surface velocity change due to declining basal motion ranges from 50 to 80% throughout the model domain (median = 59%; Fig. 7d). These model results can be viewed along the longitudinal (Fig. 8) or transverse (Fig. S5) profile coincident with the cross section discussed in Section 3.2. The surface velocity (black), ice deformation velocity (blue) and basal velocity (red) have evolved over the past 50 years to a modern, slower ice flow regime for both the longitudinal and transverse profile. Along the longitudinal profile, total surface velocity change ranges from −15 to −25 m a⁻¹ near the terminus. Throughout the model domain, most of the decline in total velocity over the past 50 years is due to a significant decline in basal motion with a modest contribution from a decline in ice deformation (Figs 7c, 8c).

Figure 6. Basal friction coefficient, C, fields inferred from model inversions for a historical (a) and modern (b) parameterization of Athabasca Glacier. Basal friction has increased throughout most of the model domain (c). Red lines in panel (a) correspond to longitudinal and transverse profiles seen in Figures 8 and S7.

Figure 7. The change in total ice surface velocity (a), ice deformation velocity (b), basal velocity (c) and the percent of total change due to declining basal motion (d) between model inversions for a historical and modern parameterization of Athabasca Glacier.

Figure 8. Surface (black), basal (red) and ice deformation (blue) velocities along a longitudinal valley profile (seen in Fig. 7) for a historical (a) and modern (b) parameterization of Athabasca Glacier. Most of the total decline in ice velocity is attributable to a decline in basal velocities (c). X = 0 corresponds to the up-glacier end of the profile closest to the icefall (Fig. 6a).

Our study area experienced a modest decline in driving stress over the past 50 years (median: −0.047 MPa) with larger declines on the lower portion and minimal declines at higher elevations (Fig. 9a). This finding is similar to previous work (Armstrong and others, 2022) which utilized a simpler flow model and found minimal change in driving stress over the past 50 years due to the counterbalancing effects of ice thinning and surface slope steepening. The moderate decline in driving stress and the resulting decline in ice deformation (Figs 7b, 8c) are not sufficient to explain the total decline in surface velocities.

Figure 9. (a) Change in Athabasca Glacier driving stress from the 1960s to the present day. Driving stress has shown a modest decline with largest changes near the terminus. (b) Force balance calculation results at five 400 m by 400 m grid cells with x = 0 shown as the orange box in panel (a). Driving stress and lateral drag show the largest declines in magnitude. (c) Proportion of driving stress balanced by each resistive term in the force balance.

The gravitational driving stress is resisted by the drag along the bed (τ _bx), lateral drag (F _lat) and resistance to flow from longitudinal stress gradients (F _lon) as described in Eqn (7) (O'Neel and others, 2005). We find that the importance of F _lon in resisting down-slope motion has declined from minor significance in the historical era to a negligible role in the present (Fig. 9b, green). The observed decline in lateral drag for all five grid cells is consistent at ~0.05 MPa (~30% of τ _dx), due to declining centerline velocity and thinning ice. F _lat has declined the most of all resistive stresses to accommodate lower driving stresses both in absolute magnitude (Fig. 9b, blue) and relative proportion of τ _dx (Fig. 9c, blue). Basal shear stress remains similar in magnitude in the modern and historical model runs, but accounts for a greater proportion of the forces counterbalancing the driving stress (Fig. 9c, orange).

4. Discussion

4.1. Ice rheology

The simplicity of the cross-sectional ice flow model allows an exhaustive search over a plausible range of the ice flow law parameters A and n. We utilize the relative simplicity of this model as a strength which allows us to explore the parameter space in a computationally efficient manner. We have confidence in these results, because the assumptions of the model (no out-of-plane gradients and negligible transverse and vertical velocities) are supported by observations. For example, an order-of-magnitude analysis shows that in-plane gradients for the longitudinal velocities (~0.1 a⁻¹) are at least an order of magnitude higher than out-of-plane gradients (~0.01 a⁻¹). We explored this flow law parameter space fully to determine which values of the rate factor and flow law exponent minimize misfit between modeled outputs and the known ice deformation rates from four historical boreholes (Figs 3, 4). In aggregate, the total misfit is minimized using values of approximately n = 3 and A = 2.7 × 10⁻²⁴ Pa⁻³ s⁻¹. We note that that there is a plausible range of these parameters (n = [2.8, 3.2] and A = [2.6 × 10⁻²⁴, 2.8 × 10⁻²⁴ Pa⁻ⁿ s⁻¹]) which produce similarly good total misfit (Fig. 4). While these values do not always provide the best fit for individual boreholes, the ice rheology values which minimize the total misfit agree reasonably with previously published values for Athabasca Glacier (n = 3 and A = 3.17 × 10⁻²⁴ Pa⁻³ s⁻¹; Adhikari and Marshall, 2012) and the commonly used values of n = 3 and A = 2.4 × 10⁻²⁴ Pa⁻³ s⁻¹ for other temperate glaciers (Cuffey and Paterson, 2010). Raymond (1971) reports best fits to individual Athabasca Glacier boreholes using values of n = [3.6, 4.2] and A = [2.2 × 10⁻²⁷, 5.5 × 10⁻³¹] Pa⁻ⁿ s⁻¹. While these higher values can produce improved fits for individual boreholes, the global misfit approach yields results which are much closer to the ‘standard’ values for temperate ice (Cuffey and Paterson, 2010). We investigate different choices of ice rheology in our friction inversion in Icepack as well to determine the sensitivity of our results to changing the rate factor A. We determined the ice velocity components along the longitudinal and transverse profile using values of n = 3 and A = 3.2 × 10⁻²⁴ Pa⁻³ s⁻¹ which are outside of the minimized RMSE error contours from our prior rheology tuning (Fig. 3). We find that changing the ice rheology parameters does not significantly change ice deformation compared to the rheology values suggested from the cross-sectional model. We see that ice deformation velocity changes along both transects are on the order of 1 m a⁻¹ (Figs S6, S7) compared with the previous inversions (Figs 8, S5). Importantly, our primary conclusion that basal motion dominates the long-term decline in Athabasca Glacier ice velocities remains the same.

Other studies have demonstrated that numerous physical mechanisms will influence the flow law exponent describing viscous ice flow. Diffusion creep (n = 1), sliding along ice grain boundaries (n = 2) and dislocation creep in the crystal lattice (n = 4) all contribute to the viscous flow of ice (Goldsby and Kohlstedt, 2001; Behn and others, 2021; Millstein and others, 2022). The flow law exponent used for a given study area provides a single lumped estimate of the contributions from all these different processes. Recent work has demonstrated that a value of n = 4 greatly improved fits to the observed ice velocity fields across numerous fast-flowing, extensional ice shelves in Antarctica, rather than the long-held ‘standard’ value of n = 3 (Millstein and others, 2022). Increased nonlinearity (n = 4) applied to ice flow in northern Greenland results in greatly reduced estimates for basal sliding in these areas (Bons and others, 2018). For Athabasca Glacier, we find that the relatively slower, compressive flow regime and warmer ice yields improved fits with a lower value of the flow law exponent (n = 3). These differing results show the importance of calibrating the ice rheology parameters which govern the relationship between applied stress and ice strain rate for the given model domain of interest.

4.2. Evolution of basal motion on decadal to centennial timescales

The results of both our cross-sectional flow model and basal friction inversion in Icepack show that the basal velocities of Athabasca Glacier have declined substantially over the past five decades. A significant strength of the higher order model implemented in Icepack is the ability to gain insight throughout our entire study reach using more complete physics to estimate ice flow. Our primary conclusion from these varied modeling approaches remains the same: that basal motion has declined significantly and constitutes the majority (50–80%) of the total observed decline in surface velocity.

We note that the results of our basal inversion in Icepack contain boundary effects that are most evident within 100 m of the valley walls (Figs 7, S5). There may also be boundary effects at the upper (icefall side) and lower (terminus side) where we specify velocity (Dirichlet) boundary conditions due to artifacts in the inversion. This may also be due to the limited available point measurements of surface velocity (Fig. S4) which we use in the loss function while inverting for basal friction. These historic surface velocity point measurements tend to be clustered along a single longitudinal profile with far fewer observations laterally and immediately adjacent to the marginal boundaries of our model domain. We caution against interpreting model outputs near these boundaries as a realistic indication of change due to these edge effects.

The presented results for basal friction are one possible solution of the inversion. As such, it is possible that some small-scale spatial details are a result of overfitting, rather than real features (Fig. 6b). The striped nature of these features indicates that they are likely an artifact of the bed topography data, which were interpolated between transverse radar profiles. We will therefore not attempt to interpret these features. If we increase our choice of alpha in the regularization function, we can achiever smoother solutions with greater model-data misfit in the loss functional (Fig. 2). If we decrease our choice of alpha, we can reduce misfit at the expense of greatly increasing the model norm which introduces roughness and larger spatial gradients to the solution (Fig. 2). We have therefore chosen values for our regularization function which attempts to optimize this trade-off using an L-curve approach (Eqn (6)). Such a trade-off is a known limitation of inverting for basal conditions from surface velocities, but this method can still provide insight to basal conditions at length scales greater than the ice thickness (Gudmundsson and Raymond, 2008). While this inversion methodology gives us insight into the macro-scale changes in basal friction that have occurred on Athabasca Glacier, one must consider this trade-off before attempting to interpret any variations at small length scales.

Both our models broadly agree with the results of a 1-D (centerline) shape factor-modified shallow ice model including parameterized longitudinal stresses described in Armstrong and others (2022). Their study reported a total reduction in surface velocity of approximately −15 m a⁻¹ (~45%) throughout the study reach and estimated the decline in basal motion as the residual between the modeled ice deformation velocity and observed surface ice velocity. Their estimates conclude that decreasing basal motion accounts for a minimum of 46% and on average 91% of the total decline in observed velocities on Athabasca Glacier over the past 50 years. The present study employs higher fidelity physical models and produces results in agreement with these earlier findings, which strengthens the conclusion that declining basal motion is primarily responsible for the observed decrease in surface velocity. The models presented in this study make fewer simplifying assumptions and rely on fewer parameterizations of physical processes, thus providing more robust results. The previous work also estimated the basal friction along a single flowline for ice geometries from 1961 and 2020 using a Weertman-style sliding law. The authors reported relatively uniform increases (150–400%) in basal friction along this longitudinal profile. Our results from Icepack inversions utilizing the same sliding law show a 100–500% (median: 127%) increase in basal friction that is not limited to the centerline, but rather increases throughout the model domain (Fig. 6).

Our overall conclusion that declining basal velocities dominate the observed long-term slowing of Athabasca Glacier conflicts with recently published work on nearby Saskatchewan Glacier which suggests increased basal sliding from the 1950s to 2010s (Stevens and others, 2023). These contradictory results on Athabasca and Saskatchewan Glacier seem particularly surprising given the similarities of the systems: both share an ice divide as outlet glaciers of the Columbia Icefield, flow through similar elevation ranges and have terminus positions only ~7 km apart. Given their proximity and obvious macro-scale similarities, such diverging long-term ice dynamics would be an interesting and perhaps unexpected finding that warrants further study. The Saskatchewan Glacier estimates rely on modern-day surface velocity using relatively short GPS time series (2–18 days) collected in August 2017 and 2019 and historical ice surface velocity measurements and ice strain rate measurements from a single, shallow borehole (43 m deep of an estimated total ice thickness of 401 m; Meier, 1957) to estimate ice rheology. We suspect some seasonal biasing resulting from the short Saskatchewan velocity measurement time span may partly explain our conflicting results. However, if Saskatchewan Glacier has in fact experienced increased basal sliding while neighboring Athabasca Glacier has experienced a considerable decline, further work is warranted to understand the nuanced mechanisms of how these glaciers have undergone such a divergent evolution in basal motion over the past 50 years.

4.3. Mechanisms causing basal motion decline

We implement a Weertman-style sliding law in our Icepack analysis to invert for the friction coefficient. We choose this implementation given that more complex sliding laws (Schoof, 2005; Gagliardini and others, 2007) require estimates of subglacial water pressure and we currently lack this dataset for our study area. Given this limitation, our analysis is limited to numerical ice flow models which are not coupled to a basal hydrology model. We cannot resolve the physical processes impacting the basal hydrology and changing the net effective pressure, which are instead lumped into a catch-all friction coefficient. A physical explanation for the increased modern friction values could be the earlier onset of melt and later freeze-up dates resulting in longer melt seasons which allow the basal drainage network to mature from a distributed to channelized network earlier in the melt season and reduce winter velocities (Van Wychen and others, 2012; Burgess and others, 2013; Tedstone and others, 2015; Hoffman and others, 2018).

This long-term increase in basal friction runs counter to observations of basal motion increases following rapid, transient increases of meltwater flow which likely pressurized a subglacial drainage network (Bartholomaus and others, 2008). Such behavior was hypothesized to result in faster ice flow in a warming world (Parizek and Alley, 2004). However, these increases in sliding caused by enhanced meltwater input are often shown to be relatively short-lived as the subglacial conduit system rapidly evolves to accommodate the higher water flux, thus increasing basal drag. Sustained, high meltwater inputs therefore do not necessarily result in prolonged periods of enhanced sliding. For example, high melt seasons have been shown to be associated with lower velocities both on alpine glaciers (Truffer and others, 2005) and the Greenland Ice Sheet (van de Wal and others, 2015). Prior work on Findelengletscher (Switzerland) showed that a significant decline (from ~110 to 55 m a⁻¹) in surface velocities from 1982 to 1994 was caused by reduced basal sliding (Iken and Truffer, 1997). The slowing basal motion was attributed to increases in basal friction, with those authors positing that the interconnectedness of subglacial cavities potentially decreased over many melt seasons, causing a decline in basal sliding. A well-connected cavity system would result in water pressure changes impacting a greater proportion of the bed surface area, while a poorly connected cavity system isolates changes in water pressure thus reducing the sliding velocity. Observations of water pressure in moulins and ice velocity in Greenland show that local ice velocity is well correlated to water pressure, but out of phase with water pressure observations from 0.3 to 2 km away (Andrews and others, 2015). These observations suggested that a decline in ice velocities later in the melt season may be due to changes in the connectivity of unchannelized portions of the bed (Andrews and others, 2015), a result supported by the modeling work of Hoffman and others (2018). Enhanced surface meltwater production (50% increase) in the Greenland ablation area has not resulted in increased ice velocities, but rather a decline in basal velocities (Tedstone and others, 2015). Additional factors other than the magnitude of meltwater production can also impact basal friction. Increases in surface slope and hydraulic gradient in Greenland primarily control the significant changes in the subglacial drainage network which increase basal friction (Maier and others, 2022). Surface slope has steepened, and hydraulic gradient has likely increased on Athabasca Glacier which would similarly explain the increases in friction predicted by our modeling results (Armstrong and others, 2022). Recent observations on valley glaciers in Alaska and the Yukon show that while a distributed drainage network tends to evolve toward a more efficient system over a melt season, significant spatial heterogeneity likely complicates the impact on basal friction (Harper and others, 2005; Rada and Schoof, 2023). Some portions of the bed remain permanently isolated which suggests large gradients in hydraulic diffusivity. Increasing drainage efficiency causes the total area of disconnected regions to increase which may have a significant impact on the basal velocities of these systems glacier. While we cannot resolve the intricacies of subglacial drainage beneath Athabasca Glacier and how they have changed over time, either increasing drainage efficiency or greater prevalence of disconnected drainage could explain the increased basal friction and associated decline in basal motion.

Constraining the exact timing, magnitude and mechanics of such long-term changes to the basal drainage network is beyond the scope of the current study. While the overall decline in basal motion of Athabasca Glacier due to increased basal friction is clear, the precise mechanisms causing these changes remain hypotheses. However, since the physical basal environment is unlikely to have changed significantly, this is likely a consequence of changing basal hydrology. Regardless of the physical causes, our results show that the change in basal velocities experienced over the past 50 years on Athabasca Glacier have provided a stabilizing feedback mechanism which reduces ice transport to lower elevations.

4.4. Stress redistribution of Athabasca glacier

Our Icepack model results allow us to explore how the lateral drag, longitudinal stress gradient and basal shear stress distribution of Athabasca Glacier have changed over the past 50 years. We find that the lateral drag has declined most significantly throughout the model domain. The longitudinal stress gradients have become almost negligible in the total force balance for the glacier's modern geometry. While the magnitude of basal shear stress has remained similar, the proportion of driving stress resisted by basal shear stress has increased along the entire transect (Fig. 9c). These results are expected for a valley glacier geometry such as Athabasca Glacier which has experienced significant ice thinning and decreased ice deformation velocities. The thinner and slower centerline ice produces less lateral drag along the valley margins and therefore relies more upon basal shear to balance the driving stress. Overall, as the glacier geometry evolves into a shallower ice regime, ice throughout the model domain relies more upon basal shear than membrane stresses to resist driving stress compared to the past. The relatively small changes in basal stress might be somewhat surprising at first, but we hypothesize that this is to some degree a reflection of the strong nonlinearity in ice rheology that demands a critical stress to be reached for meaningful ice deformation to occur. The near-constant basal stress together with the slower basal motion is consistent with our finding of higher friction coefficients. Overall, this analysis yields insight as to how the force balance is likely to change on rapidly thinning mountain glaciers.

5. Conclusions

In this work, we show that changes in driving stress cannot explain the observed decrease in Athabasca Glacier's surface velocity over the last five decades, despite substantial thinning and steepening over that period. Our cross-sectional modeling, performed at a profile collocated with historical boreholes, shows a ~30% reduction in basal motion is required to accurately predict modern surface velocities which have slowed by ~15 m a⁻¹ since the 1960s. Slowing basal motion constitutes a majority (80% on the centerline) of the total decline in observed surface velocities at this cross-valley profile.

The glacier-wide basal friction inversion suggests that the decline in basal motion is not limited to any one part of Athabasca Glacier. The inversion results show basal friction has increased (median: +127%) throughout the model domain and basal motion has correspondingly decreased. Most of the total decline in surface velocities (up to 80% near the terminus) is attributable to declines in basal motion with modest contributions due to changes in internal deformation. As basal motion has declined throughout the domain, the relative role of basal shear stress in resisting the gravitational driving stress has increased. Data from a 2022–2023 field campaign to measure basal water pressures and borehole deformation on Athabasca Glacier will provide further insight into the dynamics of this system.

Our findings show that, despite substantial 20^th century thinning (0.91 m a⁻¹ w.e.; Armstrong and others, 2022), Athabasca Glacier has benefited from a stabilizing feedback mechanism as basal velocities have declined over the past 50 years. Ice is being transported to lower, warmer elevations more slowly than in the past which could reduce mass loss rates compared to a hypothetical scenario of increased or constant basal motion. If long-term declines in basal motion are widespread across other glacial systems, such a feedback mechanism would buffer ice mass loss rates in the coming decades and centuries. Constraining basal motion throughout mountain glaciers and ice sheets will continue to prove critically important for modeling glacial response to a warming climate.