A modified viscous flow law for natural glacier ice: Scaling from laboratories to ice sheets

Here, we present a physically derived model for the parameters $n$ and $A$ in Glen’s flow law. Our model builds on a series of laboratory experiments (14, 15, 38) that support a composite flow law of the form given in Eq. 1 and provide some estimates of the relevant parameters. While the interpretation and details of these studies are still debated in the literature (12, 39, 40), the framework for and broad takeaways from our model should be valid as they depend on long-standing principles in the rheology of polycrystalline materials and empirical evidence (41–43). We further discuss assumptions made in employing this composite flow law in Discussion of this paper. A similar method of estimating $n$ from a composite flow law, derived from Eq. 1, was presented by Behn et al. (32), in which they applied a different form of a grain size model to investigate depth-varying $n$ in idealized experiments and from ice core data, to establish a physical explanation for the canonical $n=3$ estimate. Here, we build on this work to explore estimates of $n$ for varying stresses and temperatures using an idealized framework, and across the AIS using satellite observations.

Results from the laboratory studies show that the typical stresses and temperature conditions in ice sheets reduce Eq. 1 to the sum of two mechanisms, dislocation creep and grain boundary sliding. This is because diffusion creep and basal slip-accommodating grain-boundary sliding are most relevant at lower stresses than are typically found in fast-flowing regions (surface velocity $\ge 30$ $\mathrm{m}$${\mathrm{a}}^{-1}$) of ice sheets (14, 15, 38). Therefore, the composite flow law (Eq. 1) for fast-flowing regions of ice sheets simplifies to

where ${\tilde{Q}}_i$ is the activation energy for mechanism $i$, ${\tilde{A}}_{i_0}$ is the prefactor for mechanism $i$, $R$ is the ideal gas constant, $m_{\mathrm{gbs}}=1.4$ is the grain size exponent (15), and, as previously discussed, dislocation creep has no dependence on grain size. These relationships between stress and strain rate have been shown to vary with temperature such that there is a sudden increase in deformation rate for both grain-boundary sliding and dislocation creep at a critical transition temperature $T^{\ast }$. This transition is hypothesized to be a result of the onset of premelting (9, 11, 44, 45). This behavior is presently parameterized by a step-change in activation energy and prefactor for each deformation mechanism at a transition temperature, such that for each deformation mechanism, there are two values of activation energy and two values of the prefactor: one for high-temperature creep and one for low-temperature creep. We denote values applicable to low-temperature creep with a superscript${}^{-}$ and values applicable to high-temperature creep with a superscript${}^{+}$. Based on the laboratory experiments by Goldsby and Kohlstedt (15), the transition temperatures are $T^{\ast }=255$ $\mathrm{K}$ for grain-boundary sliding and $T^{\ast }=258$ $\mathrm{K}$ for dislocation creep. In applying this model, we use values of activation energies and prefactors from ref. 15 (Table 1).

We apply this composite flow law to illuminate the spatial variability in Glen’s flow law parameters ($n$ and $A$; Eq. 2). In practice, this approach may be preferable to incorporating composite flow laws directly into ice-flow models since many of these models incorporate Glen’s flow law as presented in Eq. 2. Using instead the sum of power laws (Eq. 3) can complicate the internal workings of ice sheet models, not least by requiring ice flow models to become more complex by incorporating grain size models to calculate $d$. Importantly, the method presented here obviates the need for ice-flow modelers to explicitly solve for grain size, since the grain size model is embedded here within the model results for $n,A$.

While all deformation mechanisms are active at any given time in deforming ice, the relative contributions of each mechanism set the value of the stress exponent in Glen’s flow law $n$ (15, 38). Therefore, the regions where a single mechanism is the dominant contributor to ice deformation can be represented with the respective value of the stress exponent $n$: Where dislocation creep dominates, $n\to {\tilde{n}}_{\mathrm{dis}}=4$, and where grain-boundary sliding dominates, $n\to {\tilde{n}}_{\mathrm{gbs}}=1.8$. Where neither mechanism is the dominant contributor, we expect $n$ to take on an intermediate value of $1.8<n<4$. The value of the Glen flow-rate parameter $A$ would vary according to the value of $n$ and according to ice temperature and grain size (22).

To compute values of $n$, we assume that the balance of flow mechanisms (and thus the value of $n$) vary smoothly in stress-temperature space. This is a generally reasonable assumption given that the value of $n$ is unlikely to change rapidly with strain rate. Using the scalar form of Glen’s flow law (Eq. 2), we calculate $n$ as the slope

where ${\dot{ϵ}}_{\min }=\dot{ϵ}-\mathrm{\Delta }\dot{ϵ}$ and ${\dot{ϵ}}_{\max }=\dot{ϵ}+\mathrm{\Delta }\dot{ϵ}$. We compute ${\tau }_{\min }$ and ${\tau }_{\max }$ from Eq. 3. Here, we use $\mathrm{\Delta }\dot{ϵ}={10}^{-3}\dot{ϵ}$ to be the range of effective strain rate for which the flow mechanism may remain approximately constant. The choice of $\mathrm{\Delta }\dot{ϵ}$ does not significantly affect the results for all values $\mathrm{\Delta }\dot{ϵ}\le 0.1\dot{ϵ}$ (SI Appendix, Fig. S4).

In this model of $n$, the estimated value of $n$ can change spatially and temporally based on three state variables: the strain rate, ice temperature, and grain size. To calculate $n$ going forward, we employ a combination of observations and physical models to find these state variables. Effective strain rates $\dot{ϵ}$ are observable from the gradient of surface velocity fields, which are routinely measured from remote sensing observations. To represent the dependence of deformation rate on ice temperature and grain size, two variables that affect the relative contributions of the deformation mechanisms, we couple Eq. 1 to a steady-state thermomechanical model (46) and a steady-state grain size model (47).

The thermomechanical model represents the vertical advection and diffusion of heat and the dissipation of energy through heat due to viscous deformation. Previous work using this model has shown that regions of rapid deformation (lateral margins of ice streams, basal layers of slower-flowing ice) may contain very warm ice due to the rate of shear heating (46). The steady-state grain size model represents grain size changes due to recrystallization, mechanisms that alter the size, shape, and orientation of ice crystals. This model calculates grain size as the balance between grain growth due to normal grain growth and migration recrystallization and grain reduction due to rotation recrystallization. In regions of rapid deformation and warm ice (such as lateral margins of ice streams and basal layers of slower-flowing ice), migration recrystallization typically dominates, resulting in coarse-grains (47–49). The temperature dependence of this process is parameterized in a way similar to creep, where the rate of migration recrystallization increases significantly above a critical temperature (45, 50). The physical explanation for this critical temperature is the onset of premelt enabling rapid grain boundary migration, which is the same physical process likely responsible for the shift in creep rate at high temperatures. We take the critical temperature in the grain-size model to be $T^{\ast }=255$ $\mathrm{K}$, the same critical temperature found for grain-boundary sliding in ref. 15. This is because both processes are related to grain boundary behavior and, thus, are likely to respond to temperature in a similar way. Further details of these models are discussed in Materials and Methods and can be found in the respective publications.

The models of ice temperature and grain size make some critical assumptions, explored in detail in the respective publications and summarized here. First, these models assume that temperature and grain size are in steady state. Therefore, the results shown here can be interpreted as constraining the mechanisms by which ice flows under a steady or slowly varying forcing. For grain size, this is a good assumption (47), but we incur small errors in the temperature field due to the long thermal diffusion time of ice (46). We expect the latter errors to be small relative to other assumptions and model simplifications. Second, they assume that ice advection is vertical, and therefore horizontal (both lateral and along-flow) advection is negligible. We apply this assumption to simplify the advective component of the model. Previous work suggests that our results should be reasonably accurate in Antarctica despite this simplifying assumption (46). We leave adding horizontal advective components into this model framework to future work but note that the deformation maps (Fig. 2) are valid regardless of our assumptions about advection. Only the results we show for some grounded regions of Antarctica can be influenced by our decision to neglect horizontal advection, the implications of which have been discussed in detail in previous publications (46, 51, 52). Finally, the ice temperature and grain size models assume that strain rate is constant with depth, approximating a situation in which the glacier is slipping over soft sediment or the bed is well lubricated with liquid water. This assumption limits our results in Antarctica to fast-flowing glaciers (51, 53–55), though the broader conclusions for the dependence of $n$ and $A$ on stress and temperature should be valid for areas where much of the surface velocity is attributable to deformation within the ice column (i.e., vertical shearing). As far as deformation mechanisms are concerned, we can reasonably assume that vertical shearing is the same as lateral shearing for a given ice temperature.

Our coupled model allows us to constrain the mechanisms of ice deformation in natural glacier ice and to estimate the viscous properties of ice for the full range of temperatures and stresses found in glaciers and ice sheets on Earth. We then solve for stress from $\dot{ϵ},d,T$ using an iterative nonlinear equation solver on Eq. 3.

We estimate Glen’s flow law parameters, $n$ and $A$, for varying stresses and temperatures (and thus, strain rates), and organize the results into deformation maps (Fig. 2). As shown in these maps, dislocation creep ($n=4$) dominates when stresses are above 100 kPa, while grain boundary sliding ($n=1.8$) dominates at lower stresses ($<$10 kPa), as expected from previous studies (14, 15, 38, 56).

At intermediate stresses of order 10 to 100 $\mathrm{kPa}$, a range that describes many regions in fast-flowing areas of Antarctica, the dominant creep mechanism depends strongly on temperature. For both low and high temperatures, grain-boundary sliding ($n=1.8$) dominates below $\sim$50 $\mathrm{kPa}$. Between 50 and 100 $\mathrm{kPa}$, there is a region of intermediate values $1.8\le n\le 4$. Notably, the region where both dislocation creep and grain-boundary sliding provide controls on flow ($n\sim 3$) is a narrow region in stress-space for a given temperature, suggesting that in most flow regimes, one deformation mechanism tends to dominate over the other. This has implications for ice sheet modeling, as $n=3$, representing both deformation mechanisms acting as controls on flow, is the commonly used value for the stress exponent in current ice sheet models. Furthermore, the transition between grain-boundary sliding-dominated flow to dislocation creep-dominated flow occurs at a higher stress for lower temperatures ($T<255$ $\mathrm{K}$) than it does for high temperatures ($T>260$ $\mathrm{K}$). Thus, dislocation creep controls flow for a greater range of stresses at high temperatures than at lower temperatures.

For $255\le T\le 260$ $\mathrm{K}$, dislocation creep dominates ($n=4$) for stresses above $10$ $\mathrm{kPa}$. This “dip” in the dominance of dislocation creep occurs due to the shifts in activation energy for creep. At $T=255$ $\mathrm{K}$, there is a sharp increase in activation energy for grain-boundary sliding, thus promoting the effects of dislocation creep since the activation energy for dislocation creep remains low at $255$ $\mathrm{K}$. At $T=258$ $\mathrm{K}$, the activation energy for dislocation creep increases abruptly, transitioning the flow to grain-boundary sliding-dominated flow.

At stresses below 30 kPa and temperatures colder than $255$ $\mathrm{K}$, the estimated value of $n$ is anomalously large ($n\sim 2.5$). This is due to elevated grain sizes; at such low temperatures and stresses, normal grain growth, a recrystallization mechanism that operates in regions of low deformation, causes grain size to grow. We do not expect results in these stress-temperature conditions to be realistic nor to impact the primary conclusions of this study because at such low stresses and temperatures, other deformation mechanisms (such as basal slip or diffusion creep) that are not represented in our model may play important roles.

There are a number of uncertainties underlying this model that may affect these estimates of $n$, including the choice of model for grain size, the values of $n_{\mathrm{dis}}$ and $n_{\mathrm{gbs}}$, and the parameters within $A_{\mathrm{dis}}$ and $A_{\mathrm{gbs}}$ (the prefactor $A_{0_i}$ and the activation energy $Q_i$ for each mechanism $i=\mathrm{dis},\mathrm{gbs}$). The values of $n_{\mathrm{i}}$ determine the specific values for grain-boundary sliding-dominated creep and dislocation creep-dominated creep but otherwise do not significantly affect the structure of the deformation map. Collectively, the prefactors $A_{0_i}$, the activation energies $Q_i$, and the grain sizes determine the stress at which the transition from grain-boundary sliding-dominated flow to dislocation creep-dominated flow occurs. Smaller grain sizes produce grain-boundary sliding-dominated flow at higher stresses than that for larger grain sizes (SI Appendix, Fig. S3). This is due to the grain size dependence of grain-boundary sliding in Eq. 3, in which the rate of grain-boundary sliding is inversely proportional to grain size.

The value of $n$ responds strongly to the activation energies (SI Appendix, Fig. S2). The choice of activation energy values at low and high temperatures for grain-boundary sliding and dislocation creep largely control the stress at which the flow transitions from grain-boundary sliding-dominated flow to dislocation creep-dominated flow (and therefore the stress for which $n=3$ is most applicable for a given temperature). For example, decreasing the activation energy for dislocation creep and increasing the activation energy for grain-boundary sliding would reduce the stress at which this transition occurs. A number of laboratory experiments have estimated these activation energy values for low temperatures, such as refs. 10, 15, 17, and 57, 58, 59. Laboratory studies (e.g. refs. 9 and 15) find estimates of activation energy at high temperatures as well, though there have been fewer investigations of this. Further, some evidence suggests that activation energy continues to change as a function of temperature as ice approaches its melting point (8, 20, 45), while currently it is treated as constant above the transition temperature. These results point to a need for further investigation of activation energy particularly at high temperatures to reduce uncertainties in the deformation of ice and the value of $n$. This is further explored in SI Appendix.

Changing the prefactors $A_{0_i}$ within approximately an order of magnitude produces less change in the deformation mechanism map than changing activation energies (SI Appendix, Fig. S1). This is especially relevant because the prefactors $A_{0_i}$ are most likely to change due to microstructural evolution, such as fabric development, and such development produces at most an order of magnitude change in the prefactor (60–64). This would imply that physical processes that can be parameterized as a multiplicative enhancement factor on $A$ have less of an effect on estimates of $n$ than uncertainties in activation energies or grain size estimates.

This model can provide insight into the variations in observational studies from $n\approx 2$ to 4 as a manifestation of measurements being taken at various stresses and ice temperatures (Fig. 2A). We highlight the results of a selection of observational studies with boxes in Fig. 2A. Many of the studies concluding $n=$ 2 to 3 have been done in conditions that fall below the boundary between $n=1.8$ and 4, with stresses of 10 to 100 $\mathrm{kPa}$ and temperatures $<$255 $\mathrm{K}$ (e.g., ref. 25 in Byrd Station, Antarctica and Camp Century, Greenland and ref. 65 in Devon Island Ice Cap, Canada). Studies on Taylor Glacier, Antarctica (66), and Roosevelt Island, Antarctica (67), conclude that $n$ may vary between $n=3$ and 4. These two studies considered slightly higher stresses and a wider range of ice temperatures, falling between the boundary of $n=3$ and 4 in our deformation map. Finally, studies across the upper region of the GrIS (29) and on Antarctic ice shelves (30) estimate values of $n\ge 4$ at stresses between approximately 30 to 200 $\mathrm{kPa}$. Neither study provides an estimate of ice temperature, so we take temperature to be uncertain (Fig. 2A), making direct comparison more difficult. However, for the given stresses, our model estimates $n\ge 3$ for all temperatures and estimates $n\approx 4$ for most of the stress-temperature space, suggesting reasonable agreement between these studies and our model.

The most significant deviation between observational studies and our model is with the study of Ranganathan et al. (68), which finds an estimate of $n=4.2$ in Athabasca Glacier, a temperate alpine glacier. Our model estimates a lower value of $n\approx 3$ than is found in observations at the temperatures ($T\approx 273$ $\mathrm{K}$) and stresses of the study region. Some likely causes of this discrepancy, as discussed above, are uncertainty in the estimates for high-temperature activation energy for grain-boundary sliding and dislocation creep, which have only been constrained in a handful of studies, or uncertainties in the grain sizes in temperate glaciers. The latter is especially relevant since the processes governing grain size are not as well constrained in ice at the melting temperature. The model we apply here also does not allow for values of $n>4$, which may suggest some uncertainty in the stress exponents associated with dislocation creep and grain-boundary sliding.

Applying our estimates of $n$ to Glen’s flow law, we calculate the prefactor $A$ for isotropic ice viscosity (Fig. 2B). In regions where $n\approx 4$, we estimate $A\le {10}^{-28}$$\mathrm{P}{a}^{-n}$${\mathrm{s}}^{-1}$, while where $n\approx 2$, $A>{10}^{-22}$$\mathrm{P}{a}^{-n}$${\mathrm{s}}^{-1}$. Given the difference in the stress exponent, the increase in the magnitude of $A$ for decreasing $n$ is expected because as the stress exponent decreases, $A$ must increase to produce strain rates that are of a similar order of magnitude. With estimates of both $n$ and $A$, we can estimate dynamic viscosity as $\eta ={\tau }^{1-n}/(2A)$. Dynamic viscosity varies from ${10}^{12}\text{to}{10}^{16}$$\mathrm{Pa}·\mathrm{s}$ for the typical range of stresses in ice sheets and is primarily sensitive to stress over temperature. There are no significant nonlinearities in dynamic viscosity near the transition stress between grain-boundary sliding-dominated flow and dislocation creep-dominated flow. This is because the change in the magnitude of $A$ at that transition offsets this effect, producing a roughly smooth increase in dynamic viscosity with stress.

Using this method and with reasonable estimates of $n$, strain rate, and applied stress, we can estimate ice viscosity in natural ice sheets, providing insight into the magnitude of ice softening due to mechanisms such as shear heating and recrystallization.

Our model demonstrates how fundamental rheological parameters are affected by ice flow conditions, and it enables estimates of the dominant deformation mechanisms and relevant viscous parameters across the AIS. This is possible because ice in Antarctica should be relatively dry; ultimately, we will be able to apply the model to wetter ice in Greenland once we better understand how intergranular liquid water content influences the balance of creep mechanisms. Here, we present estimates in the AIS with specific focus on Pine Island Glacier, Byrd Glacier, Bindschadler and MacAyeal Ice Streams, and Amery Ice Shelf, all of which are well-observed, fast-flowing areas that represent a range of dynamical characteristics. Computing $n,A$ requires observations of effective strain rate, ice thickness, and surface mass balance. Effective strain rates are derived from Landsat 7 and 8 velocity fields (3) using methods described in ref. 46. Ice thickness is calculated from basal topography from BedMachine (69) and surface elevation from the Reference Elevation Model of Antarctica (70). Surface mass balance, averaged over 1979 to 2019, is estimated from the Regional Atmospheric Climate Model (RACMO) (71). The estimates presented here are depth-averaged.

We estimate $n\ge 3.6$ in all fast-flowing areas of the AIS (Figs. 3 and 4). Notably, we estimate very few regions with $n=3$, the canonical value used in most ice sheet models. Within ice streams, the value of $n$ varies between $n\approx 3.6$ and $n\approx 4$ (Fig. 4). In general, we find $n\approx 4$ in regions of rapid deformation, such as the margins of glaciers, and we find $3.6\le n<4$ in slower-deforming regions, such as in the centerline of the ice shelves. In Pine Island Glacier, our model finds $n=4$ across the entire glacier, due to the magnitude of strain rates and velocity in Pine Island Glacier. Similarly in Byrd Glacier, we estimate $n=4$ for all of the grounded regions except the inflow boundary, and we estimate $n=4$ around the grounding line on the ice shelf. The value of $n$ begins to decay with distance from the grounding line on the ice shelf. Bindschadler and MacAyeal Ice Streams show more variation. On the ice shelf of Bindschadler and MacAyeal Ice Streams, we estimate $n\approx 3.7$, while we estimate $n=4$ in the grounded part of the ice streams and in the shear margins on the ice shelf. Finally, we estimate $n=4$ across Amery Ice Shelf, with a region of slightly lower values in the center of the ice shelf.

We find slightly lower estimates in the center of the Ross Ice Shelf than is found by Millstein et al. (30), which applied observations of surface strain rates and ice thickness to estimate the stress exponent $n$. Millstein et al. (30) estimated $n\ge 4$ across the ice shelf, with estimates ranging from $n=3.9\text{to}4.5$. In the same areas, we find $n\approx 3.75$. As discussed in the previous section, this could be due to parametric uncertainties, such as activation energy. In particular, underestimating the activation energy for grain boundary sliding and/or overestimating the activation energy for dislocation creep would produce lower estimates of $n$. Another likely cause is the choice of temperature and grain-size models. In particular, the ice temperature model is a one-dimensional model that neglects effects of along-flow and lateral advection of ice and assumes that ice temperature is in steady state. Similarly, the grain-size model is also one-dimensional and assumes steady-state. While in general, these are good assumptions for fast-flowing areas, these assumptions could produce small deviations from observations in the model estimates. A useful direction for future work would be to incorporate more complex temperature and grain-size models, which would also allow us to apply these estimates to slower-flowing regions (checkered regions of Fig. 3).

Our results here are also not inconsistent with the results of Behn et al. (32). They find that $n=3$ by coupling the composite flow law used in this study (15) to a grain-size model derived from ref. 74. We similarly find $n=3$ to be valid for a range of stresses and temperatures (Fig. 2A), though we estimate larger $n$ values for faster deforming regions. A key distinction between their estimates and the ones presented in this study is in the grain size model. We apply a grain-size model that accounts for different recrystallization mechanisms (47), producing coarser grain estimates than that of Behn et al. (32). This results in larger estimates of $n$. However, the framework and underlying physics explored in this study are consistent with that of Behn et al. (32), and their study similarly finds the importance of understanding grain size evolution in order to represent ice deformation in models.

We also present estimates of the flow-rate parameter $A$ across the AIS and in specific regions of the ice sheet (Fig. 3). The variations in $A$ are primarily the result of variations in $n$, since $A$ varies orders of magnitude in response to changes in $n$. For example, in Bindschadler and MacAyeal Ice Streams, we estimate $A\approx {10}^{-31}$$\mathrm{P}{a}^{-n}$${\mathrm{s}}^{-1}$ in the margins of the grounded region, while $A$ is approximately two orders of magnitude larger in the center of the ice shelf. This is likely because $n$ is larger in the margins of the grounded region.

To visualize regions of ice softening, we estimate dynamic viscosity $\eta$, which is a function of both $A$ and $n$, and thus changes in dynamic viscosity are a result of physical processes softening ice rather than calibration due to different parameter values. In Antarctic ice streams, lateral shear is primarily localized in the lateral margins, and thus in our model, the margins of ice streams are warmer and expected to have larger grain sizes than the center regions of the ice stream (47). Both of these processes affect estimates of ice viscosity. Here, we see generally that $\eta$ is lower in these rapidly deforming regions of the ice sheet. This is especially evident in the margins of Pine Island Glacier, where dynamic viscosity is approximately 1.5 orders of magnitude lower than it is in the centerline. This result is in line with a number of modeling and computational studies suggesting that ice is warmer and softer in shear margins (52, 75–79). We also expect strong (likely close to single maximum) fabric in lateral shear margins (60, 62, 63), which can decrease dynamic viscosity by an order of magnitude.

How we model ice deformation has significant effects on glacier dynamics, projections of glacier behavior, ice-sheet stability, and sea-level rise (Fig. 1). Our model provides physically informed estimates of the fundamental parameters underlying our representation of viscous ice flow: the stress exponent $n$ and prefactor $A$ in Glen’s flow law. The practical implications of our model are 1) a physical basis for scaling from laboratory ice to natural glacier ice, 2) the unification of ice deformation that captures and contextualizes the range of existing estimates of the stress exponent $n$, and 3) the establishment of a framework for estimating the values of $A$ and $n$ in Glen’s flow law (Eq. 2) based on first principles, laboratory experiments, and observations. This modeling framework can be readily applied to existing ice-flow models while respecting the various coupled physical processes, such as internal heating due to deformation and evolving grain sizes, as a way of improving our parameterization of ice deformation and as part of a broader community effort to make more reliable projections of sea-level rise.

Our results further depend on some key assumptions about ice flow. Importantly, they depend on the assumption that dislocation creep and grain-boundary sliding are both active in ice sheets, that the contributions of these two mechanisms operate independently such that their contributions can be summed according to the composite flow law (Eq. 1), and that they are the dominant two mechanisms controlling ice flow in natural conditions. Some studies have suggested that the behavior identified by Goldsby and Kohlstedt (15, 38) as grain-boundary sliding may in fact be descriptive of other processes, such as the accommodation of basal slip by grain-boundary migration, which acts as a recovery mechanism (40, 80). While more work needs to be done to determine the physical mechanism behind the $n=1.8$ regime, the analysis here should be robust to the specific mechanisms that produce the stress dependencies, since this work is interested primarily in determining under which conditions each regime is most applicable. However, uncertainties would be reduced by further investigation into which process is dominant and the incorporation of a physical understanding of that process into the model.

Relatedly, these results are further dependent upon the construction of the composite flow law itself, as presented in ref. 15. Multiterm flow laws have been presented in a number of studies to describe changes in flow mechanisms at varying stresses. These flow laws are all similar in that they are generally constructed such that a term with a lower stress exponent dominates at lower stresses and a term with a higher stress exponent dominates at higher stresses, in agreement with observational and laboratory data (13, 81–83). Many of these multiterm flow laws are empirical or phenomenological, setting up terms that do not explicitly describe physical mechanisms but rather are used to describe a changing nonlinearity of flow with stress (56, 82, 84, 85). Many of these multiterm flow laws have either two or three terms. Those with two terms generally have a linear term and a term with $n=3$, as in ref. 56. Those with three terms typically have terms with $n=$ 1, 3, and 5, as in refs. 82 and 85. While here we use the composite flow law presented in ref. 15, for future work, other multiterm flow laws could be used in this framework presented in this study to evaluate the range of $n$ estimates for varying flow laws.

Other processes not considered here include the effect of anisotropy (fabric). Fabric may increase the rate of ice deformation up to an order of magnitude, due to the development of a crystallographic preferred orientation from rotation caused by basal slip (60, 62, 86–88). Many modeling frameworks suggest that fabric can be parameterized by a strain rate enhancement factor of $\sim$2 to 10 (60, 61, 63, 89). As demonstrated in SI Appendix, Fig. S1, effects of processes, such as fabric, that can be represented as multiplicative enhancement factors are likely to affect dynamic viscosity through enhancement to the flow-rate parameter $A$ but are only likely to affect estimates of $n$ if fabric affects dislocation creep and grain-boundary sliding differently. Therefore, understanding the effects of processes such as fabric on different deformation mechanisms individually is an important step for future work.

Because we focus on the stresses and temperatures found in existing glaciers and ice sheets, this work does not account for other creep mechanisms such as diffusion creep and basal-slip accommodated grain-boundary sliding (Eq. 3). While Goldsby and Kohlstedt (15) suggest that mechanisms like diffusion creep and basal slip are unlikely to be dominant in rapidly deforming ice, these regimes should be considered when extending this work to slower-flowing regions of ice sheets, such as flow near ice divides.

While we neglect some of the complexity of ice flow, we believe this study is a step toward understanding the controls of ice flow and identifying the parameters in ice flow models that ought to be further constrained. Additional complexities in the deformation of ice could be incorporated into this framework by altering the composite flow law used.

We provide a framework for modeling the variability of the viscous flow parameters with realistic stress and temperature conditions, which provides a step toward incorporating increased physical complexity into ice flow models with minimal increases in the complexity of the models themselves. Interested readers will find the results shown in the deformation maps (Fig. 2) in an hdf5 file in SI Appendix, which can be used as a simple look-up table in ice-flow models and those used for sea-level rise projections. Such a framework should bolster the physical underpinnings of models used to study glacier dynamics and project sea level and to understand more broadly how glaciers and ice sheets evolve with climate.

To estimate the viscous stress exponent $n$ and prefactor $A$, we couple the composite flow law presented in ref. 15 with a model for ice temperature (46) and a model for steady-state grain size (47). The composite flow law explicitly represents multiple ice deformation mechanisms and is dependent upon the applied stress, the ice temperature, and the grain size, along with kinetic parameters (such as the activation energy and prefactor in the Arrhenius relation). Below, we discuss the models and relevant calculations.

We thank two anonymous reviewers and the scientific editor for their time and valuable feedback. We benefited greatly from discussions with Joanna Millstein, David Goldsby, Mark Behn, Colin R. Meyer, and Jerome Neufeld. We thank Marianne Haseloff for providing the code for the ice sheet stability model, and for providing the information for running the model. M.R. was supported by the School of Science Service Fellowship, the Martin Fellowship, National Science Foundation Directorate for Geosciences and the Natural Environment Research Council (NSFGEO-NERC) award 1853918 and the Nippon Electric Company (NEC) Corporation Fund for Research in Computers and Communications. B.M. acknowledges funding from NSFGEO-NERC award 1853918, NSF-NERC award 1739031, and the NEC Corporation Fund for Research in Computers and Communications.