Introduction
The crystal orientation fabric of glacier ice co-evolves with deformation and can locally enhance the rate of deformation by orders of magnitude (Shoji and Langway, Reference Shoji and Langway1985; Pimienta et al., Reference Pimienta, Duval and Lipenkov1987). How exactly fabric anisotropy evolves and affects the large-scale flow of ice masses has received a lot of attention in the literature, covering a broad range of topics. Questions examined include: How is the flow and age–depth relationship at ice divides affected by fabric anisotropy (Durand and others, Reference Durand2007; Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009; Pettit and others, Reference Pettit2011; Martín and Gudmundsson, Reference Martín and Gudmundsson2012)? How does fabric anisotropy affect the flow of grounded ice and ice shelves (Ma and others, Reference Ma2010)? How does fabric anisotropy affect the dynamics of mountain glaciers and ice streams compared to temperature variations (Hruby and others, Reference Hruby2020)? Do near-surface variations in fabric contain climatic information (Kennedy and others, Reference Kennedy, Pettit and Di Prinzio2013; Kennedy and Pettit, Reference Kennedy and Pettit2015)? Can fabric be used as a proxy with memory of past flow conditions (Thorsteinsson and others, Reference Thorsteinsson, Waddington and Fletcher2003; Wilson and Peternell, Reference Wilson and Peternell2011) for discovering e.g. palaeo ice streams (Lilien and others, Reference Lilien, Rathmann, Hvidberg and Dahl-Jensen2021; Llorens and others, Reference Llorens2021)? How might fabric affect inferred basal sliding and hence mass fluxes (Rathmann and Lilien, Reference Rathmann and Lilien2021)?
Addressing such questions requires being able to accurately model the large-scale flow of anisotropic ice. This demands, in turn, a bulk anisotropic rheology that can capture the relevant viscous effects of fabric anisotropy for a given application. Many anisotropic rheologies have been proposed (see e.g. Montagnat and others, Reference Montagnat2014) and often fall into one of two categories:
(1) The bulk rheology is taken to be the grain-averaged, anisotropic monocrystal rheology, for some appropriately defined average (e.g. Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009; Rathmann and others, Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen2021, and references therein).
(2) The grain c-axis orientation distribution function (ODF) is assumed to possess certain symmetries, thereby allowing a bulk rheology to be derived from e.g. plastic potential theory (e.g. Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005; Rathmann and Lilien, Reference Rathmann and Lilien2021, and references therein).
Both approaches have pros and cons; choosing one over the other is a trade off. The first approach places no constraints on what ODFs are permitted, whereas the second approach is, strictly speaking, only valid for a subset of possible ODF patterns (symmetries). Although this makes the first approach attractive, it depends on high-order fabric structure tensors (ODF moments) for flow exponents n > 1 – e.g. structure tensors through order eight for the popular flow exponent of n = 3 – which makes deriving the inverse rheology (needed for numerical modelling) too involved even for simple cases like assuming a homogeneous stress field over the polycrystal scale (Rathmann and others, Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen2021). While a linear-viscous monocrystal rheology (n = 1) or a nonlinear monocrystal fluidity that is independent of crystal orientation (Pettit and others, Reference Pettit, Thorsteinsson, Jacobson and Waddington2007) does not necessarily suffer from this problem, neglecting the observed directionally-dependent nonlinear fluidity of monocrystals (Duval and others, Reference Duval, Ashby and Anderman1983) can lead to bulk strain-rate enhancements being underestimated by at least an order of magnitude (Rathmann and others, Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen2021).
The second approach is widely adopted in ice-flow modelling but requires specifying the directional viscosities (due to fabric) in the directions of the ODF symmetry axes (elaborated on below). The problem is therefore twofold in the sense that the directional viscosities must be calculated using a separate model that depends on the orientation fabric (a so-called ‘micro-macro’ (μ-M) model, that may or may not function as an average over the grains). For example, the bulk anisotropic ‘GOLF’ rheology in the model Elmer/Ice (Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005) uses tabulated directional viscosities calculated from a separate viscoplastic self-consistent model.
The GOLF rheology in Elmer/Ice is currently the only operational anisotropic rheology of the second category above. The GOLF rheology assumes that the fabric is orthotropic; that is, the ODF has three planes of reflection symmetry with normals mi (i = 1, 2, 3) (Fig. 1a). The initial version of the rheology (Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006) considered a linear bulk viscosity, but has (e.g. Ma and others, Reference Ma2010; Licciulli, Reference Licciulli2018; Lilien and others, Reference Lilien, Rathmann, Hvidberg and Dahl-Jensen2021) been extended following Martín and others (Reference Martín, Gudmundsson, Pritchard and Gagliardini2009) (henceforth M09 approximation) and Pettit and others (Reference Pettit, Thorsteinsson, Jacobson and Waddington2007) (henceforth P07 approximation) to include a nonlinear viscosity or fluidity, respectively, taken to depend only on the second isotropic invariant of the strain-rate or stress tensor, respectively, by analogy to Glen's isotropic flow law.
The validity of assuming Glen's viscosity or fluidity is, however, untested, and does not conform with plastic potential theory where the fluidity (viscosity) should depend nonlinearly on the six orthotropic invariants of the stress (strain-rate) tensor. That is, when fabric is not isotropic but rather orthotropic, the second invariant in Glen's flow law decomposes into a weighted sum over six invariants, where the weights depend on the directional enhancement factors along the orthotropic fabric symmetry directions (elaborated on below). Indeed, this was pointed out by Ma and others (Reference Ma2010) who noted that no theoretical or experimental results are currently available to discriminate between the approximations.
In this letter, we provide such a comparison between the M09 approximation (assuming Glen's viscosity), the P07 approximation (assuming Glen's fluidity) and the unapproximated orthotropic rheology derived from plastic potential theory. In doing so, we find that the three orthotropic rheologies largely agree and are generally consistent with Dye 3 ice-core deformation tests. In addition, we give the forward and inverse analytical forms of the three orthotropic rheologies useful for future ice-flow modelling.
In the following, we start by introducing the orthotropic rheology, including the two approximations, followed by a validation and comparison of the three rheologies.
Notation
Throughout, inner, double inner and outer products, are denoted by ${\bf A}\cdot {\bf B} = \sum _{j}A_{ij} B_{jk}$, ${\bf A}\cdot \cdot \, {\bf B} = \sum _{i}\sum _{j} A_{ij} B_{ji}$ and ab = a ib j, respectively, where lower and upper case bold denote vectors and second-order tensors, respectively.
Orthotropic bulk rheology
The orthotropic bulk rheology can be constructed from plastic potential theory by demanding that the rheology be invariant under reflections with respect to mi. Objectivity implies that the forward rheology must depend on the six stress-tensor invariants of the symmetry (reflection) transformations, i.e. τ · · mimi and (τ · τ) · · mimi, or equivalently (Naumenko and Altenbach, Reference Naumenko and Altenbach2007)
Written compactly, the forward rheology is (see Naumenko and Altenbach (Reference Naumenko and Altenbach2007) for a derivation)
where the nonlinear fluidity is
Here, A is the flow-rate factor, λ i are six free material parameters, I i = I i(τ) is assumed implicit, and the index tuples are defined as
Notice that a factor of 1/2(n−1)/2 has been absorbed into A compared to the conventional definition of A in the literature.
The rheology (1)–(2) may be posed in a form that is more relevant to glaciology by expressing λi in terms of directional strain-rate enhancement factors (caused by the crystal orientation fabric), defined relative to Glen's isotropic law ${\dot {\bf \epsilon }}^{{\rm Glen}} = A( {\boldsymbol \tau }\cdot \cdot \, {\boldsymbol \tau }) ^{( n-1) /2} {\boldsymbol \tau }$. Specifically, we consider the shear and longitudinal strain-rate enhancements w.r.t. the three symmetry axes m1, m2, and m3 (Fig. 1b):
where $\hat {{\boldsymbol \tau }}( {\bf m}_{i},\; \, {\bf m}_{j})$ is an idealized compression (i = j) or shear (i ≠ j) stress-tensor function aligned with the fabric symmetry directions:
Calculating the six independent constraints from Eqn (5), gives
where E ij remain to be specified using a separate model that depends on the local orientation fabric (elaborated on below).
Inverse rheology
Posing Eqns (1)–(2) in its inverse form, ${\boldsymbol \tau} ({\dot {\epsilon}})$, amounts to solving a nonlinear matrix equation. Unlike other anisotropic rheologies, such as the transversely isotropic rheology (e.g. Rathmann and others, Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen2021), the orthotropic rheology does not tensorially depend on τ but on its projections τ · · mimj. Inverting the rheology is therefore possible by calculating the invariants $I_{l}( { \dot {\bf \epsilon }})$ (l = 1, 2, 3, 4, 5, 6) using Eqns (1)–(2) and solving for I l(τ), combined with the assumption that the inverse law should, too, depend tensorially only on ${\bf m}_{j_i}{\bf m}_{j_i}-{\bf m}_{k_i}{\bf m}_{k_i}$ and ${\bf m}_{j_i}{\bf m}_{k_i} + {\bf m}_{k_i}{\bf m}_{j_i}$. It follows from long but arithmetically straight-forward calculations that
where $I_{i} = I_{i}( { \dot {\bf \epsilon }})$ is assumed implicit, and the nonlinear viscosity is
For convenience, the shorthand γ is defined as
Indeed, we have numerically verified that the inverse satisfies ${ \dot {\bf \epsilon }}( {\boldsymbol \tau }( { \dot {\bf \epsilon }}_{{\rm 0}}) ) = { \dot {\bf \epsilon }}_{{\rm 0}}$ for various ${ \dot {\bf \epsilon }}_{{\rm 0}}$. Note that the inverse rheology was simplified by applying the identity m1m1 + m2m2 + m3m3 = I, and that $I_{j_i}-I_{k_i} = -3/2 { \dot {\bf \epsilon }}\cdot \cdot \, {\bf m}_{i}{\bf m}_{i}$.
Pettit and others (Reference Pettit, Thorsteinsson, Jacobson and Waddington2007) approximation
Pettit and others (Reference Pettit, Thorsteinsson, Jacobson and Waddington2007) proposed capturing the bulk flow nonlinearity in a linear-viscous anisotropic rheology by replacing the fluidity with that of Glen's law:
If Eqn (2) is replaced by (9), the bulk flow enhancements (5) imply that
The inverse rheology can be determined following the same procedure outlined above, yielding Eqn (6) but with the viscosity
where
Martín and others (Reference Martín, Gudmundsson, Pritchard and Gagliardini2009) approximation
Martín and others (Reference Martín, Gudmundsson, Pritchard and Gagliardini2009) proposed capturing the bulk flow nonlinearity in a linear-viscous anisotropic rheology by instead replacing the viscosity with that of Glen's law:
In order to derive the forward and inverse orthotropic rheology consistent with the viscosity (13), consider replacing the fluidity of the forward rheology (2) with
where λi′ ≠ λi. The corresponding inverse rheology, ${\boldsymbol \tau }( { \dot {\bf \epsilon }})$, is then given by Eqn (6) with Glen's viscosity (13) if
and
Relating λi to the enhancements E ij using Eqn (5), yields
and
where λi(E 11, E 22, E 33) can be inverted for numerically.
Comparison of rheologies
It is not immediately clear how the three orthotropic rheologies compare. Although the rheologies are constructed to reproduce the same directional enhancement factors when the principal stress and fabric directions are aligned by virtue of Eqn (5), the rheologies are not guaranteed to agree when the stress and fabric directions depart from these idealized ‘calibration states’ – except for n = 1 where the rheologies are identical by construction. To quantify the potential discrepancies, we consider two experiments: (1) numerically reproducing existing deformation tests made on Dye 3 ice-core samples, and (2) calculating strain-rates for synthetic ODFs that are modelled separately.
Dye 3 deformation tests
Shoji and Langway (Reference Shoji and Langway1985) investigated how strain rates depend on the misalignment between the axis of uniaxial compression (stress direction) and the preferred direction of a strong single maximum fabric in Dye 3 ice-core samples. Specifically, they considered the longitudinal strain-rate parallel to the stress direction, divided by that predicted by Glen's law for an isotropic fabric. In terms of the notation above, they calculated the strain-rate ratio
where $\dot {\epsilon }_{vv}$ is the experimentally determined quantity, v(θ) = [sin(θ), 0, cos(θ)] is the compressive axis, and θ is the angle between the compressive axis and the single-maximum direction. Only when θ = 0° are the three rheologies guaranteed (constructed) to agree (by virtue of Eqn (5)), and when θ = 90° by virtue of symmetry in the rheology. Considering intermediate angles allows, therefore, to make the differences between the rheologies clear, and to validate the rheologies against Dye 3 deformation tests.
Following Shoji and Langway (Reference Shoji and Langway1985), we take n = 3 and note that both A and the stress tensor magnitude cancel by the division in (19) (the same applies, in principle, to other isotropic contributions, such as impurities or strain-hardening, insofar as their effect can be captured in A). We furthermore assume that the orientation fabric is effectively unidirectional (all c-axes aligned with ${\hat {\boldsymbol z}}$) which is approximately true (Herron and others, Reference Herron, Langway and Brugger1985). The six bulk enhancement factors E ij (the remaining free parameters of the rheologies) are calculated using a grain-averaged, transversely isotropic monocrystal rheology that depends on the ODF, calibrated to reproduce directional enhancement factors observed from simple-shear deformation experiments (see Rathmann and Lilien (Reference Rathmann and Lilien2021) for details). Notice that the viscoplastic self-consistent approach adopted by Elmer/Ice (Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005) could equally well have been used to calculate E ij, and that E ij are constants independent of θ as far as the orthotropic rheologies are concerned.
Figure 2 shows the predicted relative strain-rates (lines) compared to those measured by Shoji and Langway (Reference Shoji and Langway1985) (markers). We find that the M09 approximation provides the best approximation to the unapproximated orthotropic rheology, but that all three rheologies compare relatively well with observations. However, the symmetry around θ = 45°, predicted by the orthotropic rheologies, is not supported by the deformation tests; that is, the deformation tests suggest that compression perpendicular to the single-maximum direction is slightly softer than compression parallel to it.
The fluidity always cancels in the division (19) for the P07 approximation, implying that the strain-rates predicted by the P07 approximation, relative to Glen's flow law (i.e. Fig. 2 and results below), are identical for all n. Hence, in the linear limit (n = 1) all three rheologies produce relative strain-rates identical to those shown here for the P07 approximation (recall the three rheologies are constructed to conform in the linear limit).
We mention in passing that the community appears lately to have become more receptive to evidence that alternative flow exponents, n ≃ 4, might be relevant in some circumstances (Bons and others, Reference Bons2018; Qi and Goldsby, Reference Qi and Goldsby2021; Millstein and others, Reference Millstein, Minchew and Pegler2022). The grey lines in Figure 2 show the corresponding behaviour for n = 4. Although the data provided by Shoji and Langway (Reference Shoji and Langway1985) is calculated assuming n = 3, and therefore cannot be used to discriminate between different n, we note that the difference between n = 3 and n = 4 is relatively small compared to n = 1 (dotted black line).
Synthetic fabrics
We also quantified the potential discrepancies between the three rheologies by calculating the strain rates predicted for a synthetic ice parcel with an evolving fabric. Specifically, we used our spectral fabric model (Rathmann and Lilien, Reference Rathmann and Lilien2021; Rathmann and others, Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen2021) to generate synthetic fabric (ODF) histories of an ice parcel subject to confined vertical compression or vertical simple shear, prescribed in terms of a constant strain-rate tensor. The two experiments consider exclusively the effect of lattice rotation (strain-induced rotation of c-axes) for ~3/4 of the total modelled parcel strains, whereas discontinuous dynamic recrystallization (DDRX) is assumed to dominate throughout the remaining deformation. We refer the reader to Rathmann and Lilien (Reference Rathmann and Lilien2021) for details on the model, but note that the (normalized) stress tensor – which determines the orientation of nucleated grains during DDRX – is assumed to be coaxial with the strain-rate tensor, and in this way no ice-flow modelling is involved. Examples of the simulated ODFs for different parcel strains, ε, are shown in Figures 3b–e and 4b–e for the two deformation experiments (red dots denote the fabric principal directions).
Given the synthetic ODFs at each degree of parcel strain, we calculate E ij in the same way as done above for the Dye 3 experiments.
Figures 3a and 4a show the longitudinal ($\dot {\epsilon }_{zz}$) and shear ($\dot {\epsilon }_{xz}$) strain rates predicted by the three rheologies, relative to Glen's law, for a constant (shared) stress-tensor that is coaxial with the strain-rate tensor used to model the parcel deformation. Note again that both the rate factor, A, and the stress-tensor magnitude cancel by considering strain rates relative to (divided by) Glen's law. We find that the largest difference between the M09 approximation and the unapproximated rheology is 5% (across both experiments), whereas for the P07 approximation the largest difference is 30%.
We point out that orthotropy might be a poor approximation for fabrics (ODFs) strongly affected by DDRX, such as seen in Figure 3e where a tetragonal symmetry is found. In this case, the double maximum pattern (two maxima in each hemisphere) introduce two additional symmetry axes along which the directional viscosities cannot be specified in an orthotropic rheology.
Discussion and conclusion
Overall, we find that both M09 and P07 approximations provide a reasonable nonlinear extension to the linear orthotropic rheology when compared to the unapproximated rheology derived from plastic potential theory. We suggest, therefore, that previous studies relying on an orthotropic bulk rheology with either the M09 or P07 approximation (see Introduction) are on relatively safe grounds even though not rigorously justified. Of course, carefully constructed deformation tests are needed to determine which version of the orthotropic rheology best agrees with the behaviour of real polycrystalline ice. For this purpose, deformation tests like those conducted on Dye 3 ice-core samples (Shoji and Langway, Reference Shoji and Langway1985) might be useful (i.e. varying the misalignment between the compressive stress direction and the preferred fabric direction). However, existing experimental results (e.g. Shoji and Langway, Reference Shoji and Langway1985, considered here) seem too ambiguous for discriminating between the rheologies (Fig. 2). In this light, the fact that the rheologies are found to produce relatively similar responses is a comforting result.
We emphasize, however, that the P07 approximation to the unapproximated orthotropic rheology is less accurate than the M09 approximation. This is particularly true for ice under compression with a stress direction that is misaligned with the preferred fabric direction (Fig. 2), possibly relevant in warm, highly stressed areas of ice sheets where DDRX is active (Figs. 3 and 4) (Duval and Castelnau, Reference Duval and Castelnau1995). As numerical ice-flow models develop support for DDRX in the future, our results provide an important caveat: the choice of a bulk flow nonlinearity in anisotropic rheologies (of the second kind described in the introduction) might affect modelled ice velocities for DDRX-induced fabrics. The two approximations might additionally lead to larger discrepancies in coupled ice-flow modelling as reported by Martín and Gudmundsson (Reference Martín and Gudmundsson2012); if small differences in viscosity lead to different deformation and fabric, minor differences could, potentially, be amplified to become significant.
We end by noting that both the forward and inverse analytical forms of all three rheologies were newly provided here in coordinate-independent form (i.e. not the fabric eigen basis) for easy use in future numerical ice-flow modelling.
Code availability
The model code (rheologies, fabric evolution and enhancement-factor calculations) is available at https://github.com/nicholasmr/specfab.
Acknowledgements
We wish to thank Editors Alan Rempel and Ralf Greve, and two anonymous reviewers for providing comments that greatly improved this paper. Notably, this includes more clearly demonstrating the differences between the three orthotropic rheologies by the numerical experiments presented here. The research leading to these results has received funding from the Villum Investigator grant nr. 16572 as part of the IceFlow project.