REASSESSING THE FLOW LAW OF GLACIER ICE USING SATELLITE OBSERVATIONS

Accurate representation of the viscous flow of ice is fundamental to understanding glacier dynamics and projecting sea-level rise. Ice viscosity is often described by a simple but largely untested and uncalibrated constitutive relation, Glen’s Flow Law, wherein the rate of deformation is proportional to stress raised to the power n. The value n = 3 is commonly prescribed in ice-flow models, though observations and experiments support a range of values across stresses and temperatures found on Earth. Here, we leverage recent remotely-sensed observations of Antarctic ice shelves to show that Glen’s Flow Law approximates the viscous flow of ice with n = 4.1± 0.4 in fast-flowing areas. The viscosity and flow rate of ice are therefore more sensitive to changes in stress than most ice-flow models allow. Mass loss from ice sheets presents both the greatest potential contribution to future sea-level rise and the largest source of uncertainty in such estimates [1, 2]. In Antarctica, mass loss occurs principally through fast-flowing glaciers that flow into floating ice shelves, which provide resistive buttressing stresses that impede the seaward flow of ice and stabilize marine grounding zones [3, 4]. The rate at which glaciers flow is controlled by the shear-thinning viscous deformation of ice [5]. The most commonly adopted constitutive relation, known as Glen’s Flow Law, is often employed to quantify the deformation of glacier ice by relating the rate of deformation, hereafter called strain rate, to the deviatoric stress [6]. Glen’s Flow Law is most simply expressed as ̇e = Aτ n e (1) where ̇e is the effective strain rate, τe the effective deviatoric stress, n the stress exponent, and A the rate factor or flow law coefficient. Variation in the parameter A can be used to encapsulate the effects of temperature, grain size, grain orientation (fabric), impurities, and interstitial water content [7]. Glen’s Flow Law is routinely implemented in large-scale ice-flow models with the prescribed value n = 3 assumed to be constant in space and time [8, 9]. However, multiple mechanisms influence the viscous deformation of ice, each with a suggested value of n: dislocation creep (n = 4), grain-boundary sliding (n ≈ 2, with slight variance dictated by the direction of motion of dislocations), and diffusion creep (n = 1) all accommodate creep at the individual grain level and, in aggregate, describe the flow of glacier ice [10]. These mechanisms are not treated independently in Glen’s Flow Law (Eq. 1). Rather, it serves as a lumped parameterization representing the combined effect of all mechanisms. Generalized forms of the flow law have been proposed to account for multiple creep mechanisms, fabric, and grain size, but these have not been widely tested, calibrated, nor implemented [9–11]. The simplicity of Glen’s Flow Law has proven useful and, subject to suitable calibration under different conditions, has the potential to provide a reasonably accurate general description of the flow of glacier ice. Glen’s Flow Law (Eq. 1) with n = 3 shows broad consistency with sparse observations of natural ice flows such as borehole deformation measurements and ice flow velocities, as well as laboratory experiments on polycrystalline ice aggregates under PREPRINT REASSESSING THE FLOW LAW conditions relevant for ice sheets [6, 10, 12–20]. However, nearly 70 years after its introduction, the need remains to test and rigorously calibrate the parameters n and A in the natural environment. Figure 1: Visual summary of our methodology. The schematic shows how we begin with publicly available satellite observations of surface velocity vector ui and ice thickness H . Using the strain rate tensor, ̇ij , we calculate the effective strain rate ̇e = √ ̇ij ̇ij/2 and along-flow strain rate ̇xx. In our areas of interest, where ̇xx ≈ ̇e, we estimate the effective deviatoric stress τe = √ τijτij/2 ≈ τxx using the force balance detailed in the supplement, which gives τxx ∝ H (Eq. 2). The values of ̇e and τe are then correlated through a flow law, indicated by the horizontal dashed arrow labeled with Glen’s Flow Law (Eq. 1). To help address this long-standing problem and benchmark a flow law that can be used in ice-flow models, we test Glen’s Flow Law across wide areas of Antarctic ice shelves, the floating extensions of the ice sheet, using satellite observations. To do so, we require independent estimates of strain rates and (deviatoric) stresses (Eq. 1). The schematic in Figure 1 graphically describes the methodology, showing how we begin with independent observations of surface velocities and ice thicknesses, apply these to evaluate strain-rate ̇e and stress τe, and then conduct a regression analysis to infer the parameters in Glen’s Flow Law. We focus on ice shelves because the underlying ocean provides negligible shear resistance to ice flow, allowing for two important simplifications in our analysis. First, we can neglect drag at the base of the ice and thus consider a stress regime that is simpler for our purposes than would be expected for grounded ice, where basal drag presents a further unknown that must be constrained. Second, the lack of drag at the base means that strain rates are approximately constant with depth. For this reason, the horizontal strain rates we calculate from observations of the surface velocity fields approximate the strain rates at all depths. Ice shelves cover areas that are large compared with the sub-kilometer resolution of observations, providing ample opportunities to comprehensively observe broad regions of flow undergoing relatively simple one-dimensional deformation. As a result, we can focus on regions that are close to being in pure extension, where the ice spreads under its own weight in one direction and the governing equations of flow reduce to a simple two-term balance, detailed further in this report. This basic premise has been employed for decades to study the rheology of glacier ice [17, 19, 21] but has not been systematically applied on continental scales before now.

Mass loss from ice sheets presents both the greatest potential contribution to future sea-level rise and the largest source of uncertainty in such estimates [1,2]. In Antarctica, mass loss occurs principally through fast-flowing glaciers that flow into floating ice shelves, which provide resistive buttressing stresses that impede the seaward flow of ice and stabilize marine grounding zones [3,4]. The rate at which glaciers flow is controlled by the shear-thinning viscous deformation of ice [5]. The most commonly adopted constitutive relation, known as Glen's Flow Law, is often employed to quantify the deformation of glacier ice by relating the rate of deformation, hereafter called strain rate, to the deviatoric stress [6]. Glen's Flow Law is most simply expressed as˙ e = Aτ n e (1) where˙ e is the effective strain rate, τ e the effective deviatoric stress, n the stress exponent, and A the rate factor or flow law coefficient. Variation in the parameter A can be used to encapsulate the effects of temperature, grain size, grain orientation (fabric), impurities, and interstitial water content [7].
Glen's Flow Law is routinely implemented in large-scale ice-flow models with the prescribed value n = 3 assumed to be constant in space and time [8,9]. However, multiple mechanisms influence the viscous deformation of ice, each with a suggested value of n: dislocation creep (n = 4), grain-boundary sliding (n ≈ 2, with slight variance dictated by the direction of motion of dislocations), and diffusion creep (n = 1) all accommodate creep at the individual grain level and, in aggregate, describe the flow of glacier ice [10]. These mechanisms are not treated independently in Glen's Flow Law (Eq. 1). Rather, it serves as a lumped parameterization representing the combined effect of all mechanisms. Generalized forms of the flow law have been proposed to account for multiple creep mechanisms, fabric, and grain size, but these have not been widely tested, calibrated, nor implemented [9][10][11].
The simplicity of Glen's Flow Law has proven useful and, subject to suitable calibration under different conditions, has the potential to provide a reasonably accurate general description of the flow of glacier ice. Glen's Flow Law (Eq. 1) with n = 3 shows broad consistency with sparse observations of natural ice flows such as borehole deformation measurements and ice flow velocities, as well as laboratory experiments on polycrystalline ice aggregates under conditions relevant for ice sheets [6,10,[12][13][14][15][16][17][18][19][20]. However, nearly 70 years after its introduction, the need remains to test and rigorously calibrate the parameters n and A in the natural environment. Figure 1: Visual summary of our methodology. The schematic shows how we begin with publicly available satellite observations of surface velocity vector u i and ice thickness H. Using the strain rate tensor,˙ ij , we calculate the effective strain rate˙ e = ˙ ij˙ ij /2 and along-flow strain rate˙ xx . In our areas of interest, where˙ xx ≈˙ e , we estimate the effective deviatoric stress τ e = τ ij τ ij /2 ≈ τ xx using the force balance detailed in the supplement, which gives τ xx ∝ H (Eq. 2). The values of˙ e and τ e are then correlated through a flow law, indicated by the horizontal dashed arrow labeled with Glen's Flow Law (Eq. 1).
To help address this long-standing problem and benchmark a flow law that can be used in ice-flow models, we test Glen's Flow Law across wide areas of Antarctic ice shelves, the floating extensions of the ice sheet, using satellite observations. To do so, we require independent estimates of strain rates and (deviatoric) stresses (Eq. 1). The schematic in Figure 1 graphically describes the methodology, showing how we begin with independent observations of surface velocities and ice thicknesses, apply these to evaluate strain-rate˙ e and stress τ e , and then conduct a regression analysis to infer the parameters in Glen's Flow Law.
We focus on ice shelves because the underlying ocean provides negligible shear resistance to ice flow, allowing for two important simplifications in our analysis. First, we can neglect drag at the base of the ice and thus consider a stress regime that is simpler for our purposes than would be expected for grounded ice, where basal drag presents a further unknown that must be constrained. Second, the lack of drag at the base means that strain rates are approximately constant with depth. For this reason, the horizontal strain rates we calculate from observations of the surface velocity fields approximate the strain rates at all depths.
Ice shelves cover areas that are large compared with the sub-kilometer resolution of observations, providing ample opportunities to comprehensively observe broad regions of flow undergoing relatively simple one-dimensional deformation. As a result, we can focus on regions that are close to being in pure extension, where the ice spreads under its own weight in one direction and the governing equations of flow reduce to a simple two-term balance, detailed further in this report. This basic premise has been employed for decades to study the rheology of glacier ice [17,19,21] but has not been systematically applied on continental scales before now.
We use measurements of ice thickness provided through the BedMachine project [22], and surface velocity data from the NASA Inter-mission Time Series of Land Ice Velocity and Elevation (ITS_LIVE) project [23]. The surface velocity data, which encompass most of the Antarctic Ice Sheet at a grid spacing of 120 m×120 m, are derived from Landsat 4, 5, 7, and 8 imagery using the auto-RIFT feature tracking processing chain, providing reliable constraints on the two horizontal components of ice velocity [23]. We use these to calculate the horizontal strain rates˙ ij (for i, j = x, y the two horizontal coordinates) across all Antarctic ice shelves, as defined by 2˙ ij = (∂u i /∂x j + ∂u j /∂x i ), where u i represents the horizontal components of the ice velocity vector and x i the horizontal coordinates. To calculate the components of the velocity gradient, we apply a two-dimensional Savitzky-Golay filter with a polynomial order of one and square window of 3720 m (31 pixels) [24]. More detail on the strain rate calculations are given in the supplement.
After deriving strain rates from the surface velocity fields, we determine regions flowing in approximately pure extension, with a view to simplifying the force balance governing the local ice flow. The two-dimensional strain rate tensor˙ ij has three unique components (the off-diagonal terms are equal by definition) and a scalar invariant representing the effective horizontal strain rate˙ = ˙ ij˙ ij /2, where summation is implied for repeated indices. Note that the effective strain rate˙ e in Eq. 1 follows the same definition as for˙ but is applied to the three-dimensional strain-rate tensor. We focus on areas of the ice shelves that are solely confined by seaward pressure in the along-flow, or x, direction, and analyze areas in which the along-flow component of the strain rate tensor˙ xx is much larger than both lateral normal and shear strain rates (˙ xx ˙ yy ,˙ xy ). We combine these into the single criterion˙ xx ≈ √ 2˙ , corresponding to areas of the ice shelves where longitudinal extension is dominant. The more specific criterion˙ xx >˙ is used to define large, spatially coherent regions where the extensional component of strain dominates the flow ( Fig. 2 and Supplementary  Figures S1 and S2). Approximately 25% of the total surface area of all Antarctic ice shelves satisfies this criterion. In these areas it follows from the incompressibility of ice and the absence of drag at the base of ice shelves that˙ xx ≈˙ e , the three-dimensional effective strain rate in Eq. 1.
To estimate the effective deviatoric stress from remote sensing observations, we utilize a well-established reduced form of the Stokes equations that govern the viscous flow of glacier ice. Over the ice shelves, where negligible shear stress applies at both the upper (atmosphere) and lower (ocean) surfaces of the ice, we can adopt the depth-integrated form of the Stokes equations commonly referred to as the Shallow-Shelf Approximation (SSA), which contains only body forces and the horizontal gradients of the stress tensor elements. Based on the conditions described above, we can further reduce the SSA equations to a simple expression relating the along-flow deviatoric stress τ xx to local ice thickness H as: where g = g(1 − ρ/ρ w ) is the reduced gravity, representing the balance between the resistive longitudinal stress and the driving buoyancy force (see supplementary materials for full derivation). Here, we take ρ = 910 kg/m 3 as the mass density of glacier ice and ρ w = 1026 kg/m 3 as the mass density of seawater. Where the criteria for predominantly extensional flow is met (˙ xx ≈˙ e ), we expect τ xx ≈ τ e . Thus, the criteria we apply to the strain rate fields to identify areas in primarily extensional flow allows us to calculate effective stress τ e (Eq. 1) from observations of ice thickness and independently of the surface velocity fields used to calculate˙ e . And our focus on a single flow regime allows us to avoid the complexities that arise from viscous anisotropy in ice, which would require a non-scalar form of A to represent deformation in multiple directions.
Linear regressions fitted to the values for log(˙ e ) and log(τ e ) constrain n through the slope and A in the y−intercept, divulging values of the flow law parameters across viable regions of Antarctic ice shelves. To determine 95% confidence intervals on the regression of strain rate on stress, we implement a non-parametric bootstrap which allows us to estimate constraints on the determined value of n without making assumptions on the underlying structure of the distribution [25]. Our analysis encompasses regions of both large ice shelves, such as those shown in Figure 2, and smaller ice shelves that line the continent. We focus first on highlighted areas on the Ross and Filchner-Ronne Ice Shelves in Figure  2, which we extracted from areas along flow lines, with probable consistency between values of temperature, grain size, and fabric, and therefore A and n.
The log-log plots between strain rate and deviatoric stress shown in Figure 2 exhibit linear trends that are consistent with a power-law relation. These results provide strong evidence that, for a suitable choice of n, Glen's Flow Law is a valid approximation of the viscous flow of Antarctic ice shelves and, as discussed later, likely other dynamic regions of Antarctica. Critically, we find n ≈ 4 in the fast-flowing, extensional regions of Antarctic ice shelves. This result is consistent with early evidence for a higher value of the flow law exponent [6,19,21], and demonstrates that this higher value is applicable to natural ice flow at the continental scale.
The results of our full analysis covering all viable regions of Antarctic ice shelves is shown in Figure 3, which includes regions of both large and small ice shelves (mapped in Supplementary Figures S1 and S2). The normalized kernel density estimates of the bootstrapped values of the flow law exponent (Figure 3) indicate that n = 4.1 ± 0.4 in extensional regions of Antarctic ice shelves. Figure 3 shows the confidence with which our estimate stands across  geographic areas of different sizes and representing a range of stresses. Large areas extracted for analysis, > 1000 km 2 , have less spread in the error estimation and are centered closer to n = 4.1, whereas smaller areas exhibit a greater spread in the distribution. This is likely because the broader ranges of stresses and greater number of observations in the larger ice shelves provide relatively accurate inferred trends across the data. Notably, geographic regions from West Antarctica have higher values of n than regions sampled from East Antarctica. This observation can be plausibly attributed to higher sub-ice-shelf melt rates in West Antarctic ice shelves, where the bulk of ice is created on the ice shelf by compaction of snow as opposed to being inherited from the grounded glacier [26]. Additionally, there is a possible grain size dependence wherein warmer conditions would contribute to larger grains [27]. In such regions, larger grains, strain, and values of the stress exponent validate a hypothesis that ice deformation is facilitated primarily by dislocation creep [10,28]. Our results highlight further spatial variability in the precise values of the flow law exponent and rate factor across different ice shelves, and even different regions within single ice shelves (see Figure 3). We reserve for future work detailed analysis and modeling of these variations.
Values of the flow law rate factor, A, in our results are smaller than standard tabulated values for n = 3 [29], which span 10 −35 to 10 −27 Pa −n s −1 (see Supplement Figure S3). With the higher (integer) value n = 4, it follows that smaller values of the rate factor A relative to those for n = 3 are required to accommodate the same strain rate at a reference stress. Here, we do not give newly calibrated values for A because proper constraints on the physical properties of the ice, like temperature and grain size, are not currently available in these areas and require work that is beyond the scope of this study. Rather, we note that the smaller values found here validate our method for deriving Glen's Flow Law and we recommend that future efforts using a value n ≈ 4 utilize standard tabulated sources for A [29] and scale these values accordingly for the new value of n.
The result that n ≈ 4 challenges the long-held practice of assuming the flow law exponent is n = 3 everywhere, and at all times, in large-scale ice-sheet flow models. While our observations focus on regions that make up about a quarter of the extent of the ice shelves and experience stresses of order 100 kPa (Supplement Figure S4), complementary laboratory work showing that n = 4 is suitable at higher stresses [10] supports extending our conclusion that n ≈ 4 to other dynamic regions in Antarctica. Additionally, our conclusion complements a growing body of work advocating for the use of n > 3 in other areas of the cryosphere [13,14]. Taken together, this work calls for a broader community effort to quantify the uncertainties in the flow law parameters and the consequences of these uncertainties on models of glacier dynamics. A higher value of n increases the sensitivity of viscosity to changes in stress but the impact of n = 4 on large ice-flow models used for projections of sea-level rise and ice-sheet evolution remains unclear as few sensitivity analyses have been conducted [9] and n is not a parameter explored in current ensemble-model analyses [1,2]. The value n = 4 has the potential to increase the sensitivity of ice-sheet mass loss to ongoing climate change considerably relative to n = 3 due to the stronger dependence of flow rates to changes in resistive stresses.
By applying continental-scale satellite observations to standard models in glacier dynamics, we have validated Glen's Flow Law, a constitutive relationship that helps form the foundation of modern glaciology, and calibrated the stress exponent in Antarctic ice shelves. This work serves as a pathway towards a standard calibration framework for the community using publicly available remote sensing data. Our conclusion that n ≈ 4 across much of Antarctica's ice shelves is a step towards reassessing the governing equations of ice flow in the satellite age, and reveals an increased sensitivity of flow rates to applied stresses relative to the commonly used n = 3. As a consequence, future sea-level rise is likely more sensitive to climatic forcings than present models using common assumptions of the flow law allow.