## Introduction

The viscosity of polycrystalline ice depends on physical properties such as the orientation distribution of individual grains (orientation fabric) (Shoji and Langway, Reference Shoji and Langway1985, Reference Shoji and Langway1988), grain sizes (Goldsby and Kohlstedt, Reference Goldsby and Kohlstedt1997; Kuiper and others, Reference Kuiper2020*a*, Reference Kuiper2020*b*; Fan and others, Reference Fan2020), and hardening effects due to crystal lattice defects (Faria and others (Reference Faria, Weikusat and Azuma2014) and references therein). The orientation fabric is particularly important because ice monocrystals deform by dislocation slip (shear) along crystallographic planes (Fig. 1 left), where basal plane shear (planes with normal **c**) is up to 100–1000 times softer than along any other plane (Weertman, Reference Weertman, Whalley, Jones and Gold1973; Duval and others, Reference Duval, Ashby and Anderman1983). As a consequence, anisotropic orientation fabrics can introduce local directional hardening of a polycrystal (Thorsteinsson (Reference Thorsteinsson2001, Reference Thorsteinsson2002) among others), which in turn can affect the deformation of ice masses at large scales (Thorsteinsson and others, Reference Thorsteinsson, Waddington and Fletcher2003; Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006; Pettit and others, Reference Pettit, Thorsteinsson, Jacobson and Waddington2007; Martín and Gudmundsson, Reference Martín and Gudmundsson2012; Hruby and others, Reference Hruby2020).

General monocrystal rheologies for ice, such as those proposed by Castelnau and others (Reference Castelnau, Duval, Montagnat and Brenner2008) and Suquet and others (Reference Suquet2012), should be able to account for the variation in shear resistance that occurs between different slip systems, including the observed non-linear power law behaviour (Kamb, Reference Kamb1961; Duval and others, Reference Duval, Ashby and Anderman1983). A popular simplified rheology is the transversely isotropic rheology of Johnson (Reference Johnson1977) (see Supplementary A and the notation section):

which approximates monocrystals as axisymmetric with symmetry axis **c** (Meyssonnier and Philip, Reference Meyssonnier and Philip1996). Here ${{\bf \tau}'} \,\hbox{and}\,\,{{ \dot {\bf \epsilon}}}'$ are the microscopic deviatoric stress and strain-rate tensors, *A*′ is an isotropic rate factor, and *n*′ is the power-law exponent. The enhancement factors ${{E}_{cc}'}$ and ${{E}_{ca}'}$ are the strain-rate enhancements for compression/extension along **c** and shear parallel to basal planes (Fig. 1):

where **a** denotes any direction transverse to **c**, and the indices *c* and *a* indicate the tensorial components in the direction of **c** and **a**, respectively.

By averaging the rheology (1)–(2) over a polycrystal (defined below), it is possible to construct a bulk (polycrystalline) rheology that provides an anisotropic extension to Glen–Nye's isotropic flow law (Svendsen and Hutter, Reference Svendsen and Hutter1996), or to infer bulk directional enhancement factors, the latter a prerequisite for bulk rheologies such as the orthotropic flow law proposed by Gillet-Chaulet and others (Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005). Past work has, however, focused on the linear (*n*′ = 1) monocrystal rheology (Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006; Durand and others, Reference Durand2007) or neglected the orientation-dependent contributions to the non-linear fluidity ${{\eta^{-1}}}$ (Pettit and others, Reference Pettit, Thorsteinsson, Jacobson and Waddington2007; Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009; Ma and others, Reference Ma2010; Martín and Gudmundsson, Reference Martín and Gudmundsson2012; Hruby and others, Reference Hruby2020) by setting

For strongly developed fabrics, bulk directional enhancement factors calculated using (1) with (5) (Gödert and Hutter, Reference Gödert and Hutter1998; Gagliardini and Meyssonnier, Reference Gagliardini and Meyssonnier1999) can meanwhile lead to smaller enhancements compared to basal slip (Schmid) models that include an orientation-dependent fluidity (Thorsteinsson, Reference Thorsteinsson2001). Bulk rheologies that rely on (5) as the grain fluidity might therefore produce different flow regimes compared to (2).

It is not clear over which regimes of orientation fabric strength (5) provides a good approximation to (2) (in the grain-averaged rheology), or what fabric information the grain-averaged rheology requires if (2) is to be used. In this letter, we address these questions by calculating bulk directional enhancement factors for single-maximum and girdle fabrics using both the full (2) and approximate (5) fluidities. In doing so, we find that the full fluidity (2) requires specifying more detailed information about the orientation distribution of grains (higher-order structure tensors), which leads us to propose a new non-parametric spectral fabric model instead of directly modelling the structure tensors.

### Notation

Throughout, primes shall be used to denote rheological parameters, stresses and strain rates of monocrystals, as opposed to non-primed variables used to denote bulk rheological parameters, stresses and strain rates of polycrystals. Boldface symbols denote vectors or tensors, the order of the latter being implicit by the context. Let **a** and **b** be vectors, and **A** and **B** be second-order tensors. The following products are then used: ${\bf a}{\bf \cdot} {\bf b} = a_{i}b_{i}$, **ab** = *a* _{i}*b* _{j} (the outer, dyadic product), ${\bf A}{\bf \cdot} {\bf a} = A_{ij}a_{j}$, ${\bf A}{\bf \cdot} {\bf \cdot} {\bf a}{\bf b} = A_{ij}a_{j}b_{i}$, ${\bf A}{\bf \cdot} {\bf B} = A_{ij}B_{jk}$, and ${\bf A}{\bf \cdot} {\bf \cdot} {\bf B} = A_{ij}B_{ji} = {\rm tr}[ {\bf A}{\bf \cdot} {\bf B}\rpar$, where summation over repeated indices is implied. The identity matrix is denoted by **I**, and the superscript ${\ssf T}$ denotes the matrix transpose.

## Averaged rheology

Following Azuma (Reference Azuma1994) (among others), we define the bulk strain rate, ${ {\dot {\bf \epsilon }}}$, as the average strain rate of grains in a polycrystal. Moreover, we adopt the Sachs hypothesis (Sachs, Reference Sachs1929) that the microscopic stress field, ${{\bf \tau}'}$, can be assumed homogeneous over the polycrystal and equal to the bulk stress at the polycrystal scale, ${{\bf \tau}}$. Taken together, the bulk strain rate is therefore

By specifying the individual grain stresses in this way, grain strain compatibility (aligned shared grain boundaries) is not guaranteed, and mechanical interactions between grains are not directly accounted for.

The Sachs hypothesis can broadly be regarded as an end-member case along a spectrum of homogenization schemes. At the opposite end, the Taylor hypothesis (Taylor, Reference Taylor1938) assumes a homogeneous strain rate over the polycrystal. The Taylor hypothesis is, however, not well adapted for strongly anisotropic fabrics: five independent slip systems are required to accommodate for an arbitrary bulk strain rate (${{\dot {\bf \epsilon }}}$ possesses five independent components), but basal planes provide only two (Castelnau and others, Reference Castelnau, Duval, Lebensohn and Canova1996). Although the realized stress and strain-rate fields inside polycrystals are somewhere between the two end-member cases, observations and experiments (Azuma and Higashi, Reference Azuma and Higashi1985; Azuma, Reference Azuma1995) suggest that assuming a homogeneous stress field is a better approximation of the two. For brevity, we refer the reader to Thorsteinsson (Reference Thorsteinsson2001), Castelnau and others (Reference Castelnau, Duval, Lebensohn and Canova1996) or Montagnat and others (Reference Montagnat2014) for details and more sophisticated homogenization schemes.

Given a prescribed stress that all grains are subject to, equation (6) reduces to an average over all grain orientations. Provided with a (unnormalized) grain c-axis distribution ${n(\theta\comma\; \phi)}$ defined over orientation space *S* ^{2}, the average strain rate of grains becomes

where dΩ = sin*θ*d*θ*d${\phi}$ is the infinitesimal solid angle, and

is the total number of grains. Unless stated otherwise, ${\langle{\cdot}\rangle\,{ = }\,\langle{\cdot}\rangle_{n(\theta\comma \phi)}}$ is hereafter assumed implicit.

Averaging the rheology (1)–(2) using (7), the bulk strain rate (6) is

where **c**^{k} is the *k*-th outer (dyadic) product of **c** with itself. In the linear case, *n*′ = 1, the averaged terms reduce to ${\langle {\eta^{-1}}{\bf c}^{k} \rangle}\,{ = }\,{\langle A' {\bf c}^k \rangle}\,{ = }\,A'{\langle{\bf c}^k \rangle}$, and determining $\langle {{\dot {\bf \epsilon }'}}\rangle$ therefore requires constructing the second- and fourth-order structure tensors, 〈**c**^{2}〉 and 〈**c**^{4}〉 (Advani and Tucker, Reference Advani and Tucker1987). Likewise, for the orientation-independent fluidity (5), the averaged terms reduce to $\langle {{\eta^{-1}}{\bf c}^{k}}\rangle = {A}'[ {\bf \tau }'{\bf \cdot} {\bf \cdot} \, {\bf \tau }'\rpar ^{[ {n}'-1\rpar /2}\langle {{\bf c}^{k}}\rangle$, and hence $\langle {{\dot {{\bf \epsilon}}'}}\rangle$ depends on the same structure tensors.

As noted, the rheology of monocrystals has been experimentally found to follow a power law with an exponent of 2 to 3 (Duval and others, Reference Duval, Ashby and Anderman1983). In the non-linear case of *n*′ = 3, the averaged terms in (9) become (for any integer *k*)

and determining $\langle {{\dot {{\bf \epsilon}}'}}\rangle$ therefore requires constructing the structure tensors 〈**c**^{2}〉, 〈**c**^{4}〉, 〈**c**^{6}〉 and 〈**c**^{8}〉. That is, adopting the orientation-dependent fluidity (2) with *n*′ = 3 over the orientation-independent fluidity (5) requires additionally specifying the higher-order structure tensors 〈**c**^{6}〉 and 〈**c**^{8}〉.

We mention in passing that while considering the Taylor hypothesis, too, would provide a more complete treatment of the problem, averaging the inverse rheology, $\langle {{\bf \tau }'[ {{\dot {\bf \epsilon }}}\rpar }\rangle$ (see Supplementary A), for *n*′ > 1 involves inverse fractional integrands that are hard to rewrite in terms of even-ordered structure tensors in analogy to (10).

With increasing order, structure tensors quantify the structure of $n(\theta\comma\, \phi)$ at an increasingly finer scale on *S* ^{2}. It is, however, not clear over which regimes of fabric strength that (5) is a good approximation to (2). In what follows, we address this question by calculating bulk enhancement factors both assuming (5) and while relaxing this assumption.

## Bulk directional enhancement factors

The rheology of an anisotropic polycrystal can be characterized by its bulk enhancement factors for compression/extension and shear w.r.t. the fabric principal directions. Given $\langle {{\dot {\bf \epsilon }'}[ {\bf \tau }\rpar }\rangle$, Thorsteinsson (Reference Thorsteinsson2001) suggested defining the bulk enhancement of the **v**–**w** component of ${{\dot {\bf \epsilon }}}$ as (**v** and **w** being arbitrary vectors)

where $\langle {{{\dot {\bf \epsilon }}}'[ {\bf \tau }\rpar }\rangle_{{\rm const.}}$ is the average strain-rate tensor of an isotropic polycrystal ${(n(\theta\comma\; \phi) = {\rm const.})}$.

It is important to note that *E* _{vw} does not depend on *A*′ or the stress magnitude because they are orientation independent and therefore cancel by virtue of the division in (11). Likewise, if $\eta^{-1}$ is approximated by (5) (orientation independent), then *E* _{vw} does not depend on $\eta^{-1}$ either. As a consequence, *E* _{vw} derived from the linear rheology (*n*′ = 1) is identical to that derived from any non-linear rheology (*n*′ > 1) with (5). It therefore suffices to compare the bulk enhancements resulting from the linear (*n*′ = 1) and non-linear (*n*′ = 3) realizations of the rheology (1)–(2).

For the sake of simplicity, we restrict ourselves to consider axisymmetric fabrics, which includes single-maximum and girdle fabrics, and denote the fabric symmetry axis **m** (Fig. 2). Complementary to **m**, let **t** denote any transverse direction to **m**, and let **p** and **q** denote the ± 45^{°} directions to **m**: ${\bf p} = [ {\bf m} + {\bf t}\rpar /\sqrt {2}$ and ${\bf q} = [ {\bf m}-{\bf t}\rpar /\sqrt {2}$. Mechanical deformation tests have found that shearing parallel to the basal planes is about ten times easier for a polycrystal with aligned c-axes compared to isotropic ice (Pimienta and others, Reference Pimienta, Duval, Lipenkov, Waddington and Walder1987), and that shearing of polycrystals with a preferred-direction fabric can be up to 10^{3}–10^{4} times softer along the plane with normal **m** (i.e. **m**–**t** shear) compared to the hardest direction, oriented at ≃ 45^{°} to **m** (i.e. **p**–**q** shear) (Shoji and Langway, Reference Shoji and Langway1985, Reference Shoji and Langway1988). Translated into enhancement factors, the two observations suggest that *E* _{mt} ≃ 10 and ${E}_{mt}/{E}_{pq} \simeq$ 10^{3}–10^{4}. We shall therefore be interested in determining *E* _{mt} and *E* _{pq} for axisymmetric polycrystals, and, in addition, the compression/extensional enhancement factor *E* _{mm}. By definition, *E* _{mm} and *E* _{mt} are the bulk enhancements of a polycrystal when subject to bulk pure- and simple-shear stresses w.r.t. **m** and **t**:

where ${\tau_0}$ is an arbitrary stress magnitude. Likewise,

A summary of the enhancement factors and defined directions is provided in Figure 2.

At this stage, an important caveat regarding the bulk strain rate (9) is warranted: as pointed out by Meyssonnier and Philip (Reference Meyssonnier and Philip1996), the mechanical environment of grains in a polycrystal is very different from that of isolated grains (monocrystals) for which (1)–(2) is assumed. Adopting values of ${{E}_{cc}'}$ and ${{E}_{ca}'}$ determined from mechanical tests on monocrystals can not necessarily be expected to produce the observed bulk enhancement factors; grain interactions can result in non-preferred slip systems being activated in proportions that are inconsistent with the values ${{E}_{cc}'}$ and ${{E}_{ca}'}$ derived from monocrystal mechanical tests. We shall therefore suppose that ${{E}_{cc}'}$ and ${{E}_{ca}'}$ are unknown. For a fair comparison between the linear and non-linear rheology (and hence (5) and (2)), ${{E}_{cc}'}$ and ${{E}_{ca}'}$ must be selected for both rheologies such that *E* _{mt} and *E* _{mt}/*E* _{pq} match observations in the end-member case of a maximally strong fabric (aligned c-axes). Only after tuning against this end-member case can the bulk directional enhancement factors that result from linear and non-linear grain rheologies be compared for intermediate fabric strengths and the approximation (5) assessed.

### Unidirectional fabric

In a maximally strong fabric, c-axes are aligned with the symmetry axis, **c** = **m**, and hence ${n}{[\theta \comma\; \phi \rpar } = \delta [ \hat{\bf r}-{\bf m}\rpar$, where $\delta [ \hat{\bf r}\rpar$ is the Dirac delta function and $\hat{\bf r} = [\sin(\theta)\cos(\phi)\comma\; \, \sin(\theta)\sin(\phi)\comma\; \, \cos(\theta)]$ is the radial unit vector. For this distribution, the numerator in (11) is easily determined from (9)–(10) by noting that $\langle {{\bf c}^{k}}\rangle _{\delta [ \hat{\bf r}-{\bf m}\rpar } = {\bf m}^{k}$. Calculating the denominator in (11) is more tedious but nonetheless possible by rewriting **c**^{k} in terms of spherical harmonics, $Y_{l}^{m}{[ \theta \comma\; \phi \rpar }$, and using the contraction rule (writing the product of any two spherical harmonics in terms of a new spherical harmonics series). In order to stay focused on the physics, we refer the reader to Supplementary B for details.

In Figure 3, *E* _{mt} (blue-filled contours), *E* _{mt}/*E* _{pq} (solid contours) and *E* _{mm} (dashed contours) are shown as a function of ${{E}_{cc}'}$ and ${{E}_{ca}'}$ in the case of the linear (Fig. 3a) and non-linear (Fig. 3b) grain rheology. With increasing ${{E}_{ca}'}$ and decreasing ${{E}_{cc}'}$, the soft-to-hard shear-enhancement ratio, *E* _{mt}/*E* _{pq}, is generally found to increase, *E* _{mm} is found to decrease, while *E* _{mt} monotonically approaches a maximum value (dark blue). Overall, the orientation-dependent fluidity (*n*′ = 3), compared to the orientation-independent case (*n*′ = 1, or assuming (5)), yields stronger enhancement factors (increased sensitivity to fabric strength). While it is possible to carefully select ${{E}_{cc}'}$ and ${{E}_{ca}'}$ for both *n*′ = 1 and *n*′ = 3 such that *E* _{mt}/*E* _{pq} ≃ 10^{3}–10^{4} according to Shoji and Langway (Reference Shoji and Langway1985, Reference Shoji and Langway1988), the value of *E* _{mt} is limited to *E* _{mt} ≤ 2.5 for *n*′ = 1 and *E* _{mt} ≤ 4.375 for *n*′ = 3, the latter value being approximately a factor of 2 smaller than experimentally determined values (Pimienta and others, Reference Pimienta, Duval, Lipenkov, Waddington and Walder1987).

### Evolving single-maximum and girdle fabrics

Putting aside the upper bounds on *E* _{mt} (treated in the discussion), the question remains whether (5) can provide a good approximation to fluidity (2) (in the grain-averaged rheology) if they are both tuned to reproduce *E* _{mt}/*E* _{pq} ≃ 10^{3}–10^{4} for a unidirectional fabric. For this reason, let us consider a single-maximum and girdle fabric which evolve from isotropy, a situation that is common in polar ice masses.

During vertical pure shear, grains in a polycrystal (ice parcel) tend to rotate towards the compressive axis (${\hat{\bf z}}$) and away from the extensional axis (Azuma and Higashi, Reference Azuma and Higashi1985; van der Veen and Whillans, Reference van der Veen and Whillans1994), the former coinciding with the symmetry axis, **m**. The preferred direction (single maximum) therefore strengthens with increasing parcel strain (decreasing ice-parcel height) in the absence of other recrystallization processes (Alley, Reference Alley1992; Thorsteinsson, Reference Thorsteinsson2002). This corresponds to the situation found at ice-sheet domes (Haefeli, Reference Haefeli1963; Nye, Reference Nye1963; Alley and others, Reference Alley1995) where the vertical strain-rate, $t_e^{-1}$, is taken to be approximately uniform, i.e. ${\bf \nabla }{\bf u} = t_e^{-1} {\rm diag}[ 1/2\comma\; \, 1/2\comma\; \, -1\rpar$. For girdle fabrics, c-axes are distributed along the preferred plane with normal **m**, which is the result of extension along **m** and compression in the transverse (preferred) plane in the absence of other recrystallization processes (Alley, Reference Alley1992). The corresponding velocity gradient is the time-reversed of what forms a preferred-direction fabric, i.e. *t* _{e} → −*t* _{e}. The girdle therefore strengthens with increasing ice-parcel height.

In the absence of other microstructural processes, grain rotation can be modelled as a (conservative) convection process on *S* ^{2} involving ${n(\theta\comma\; \phi)}$ (Gödert and Hutter, Reference Gödert and Hutter1998):

where $\dot {{\bf c}}{[ \theta \comma\; \phi \rpar }$ is the c-axis velocity field on *S* ^{2}, and the divergence operator acts on *S* ^{2}. Let us consider the simplest case where grains rotate in response to the macroscopic stretching, ${{\dot {\bf \epsilon }}} = [ {\bf \nabla }{\bf u} + [ {\bf \nabla }{\bf u}\rpar {}^{\mkern -1.5mu{\ssf T}}\rpar /2$, and spin, ${\bf \omega } = [ {\bf \nabla }{\bf u}-[ {\bf \nabla }{\bf u}\rpar {}^{\mkern -1.5mu{\ssf T}}\rpar /2$, which allows the detailed microscopic stress and strain-rate fields to be neglected and hence interactions between neighbouring grains to be disregarded. By moreover neglecting higher-order dependencies on the velocity gradient, and requiring that basal planes preserve their orientation when subject to simple shear (like a deck of cards), it has previously been proposed that (Svendsen and Hutter, Reference Svendsen and Hutter1996; Gödert and Hutter, Reference Gödert and Hutter1998)

for an arbitrary c-axis ${{\bf c} = {\bf c}{(\theta,\phi)} = [\sin(\theta)\cos(\phi)\comma\; \,\sin(\theta)\sin(\phi)\comma}$ ${\cos(\theta})]$.

Notice that (15)–(16) represents a separate kinematic deformation experiment which, given ${\bf \nabla }{\bf u}$, provides an idealized time-evolution of ${n(\theta\comma\; \phi)}$ (and hence 〈**c**^{k}〉) as input for $\langle {{\dot{\bf \epsilon }'}[ {\bf \tau }\rpar }\rangle$; that is, $\dot {{n}}{[ \theta \comma\; \phi \rpar }$ does not depend on the grain rheology. When dynamically modelling ice flow, the two should be coupled, but for the present purpose it suffices to consider the bulk enhancement factors *resulting* from a given orientation fabric distribution, ${n(\theta\comma\; \phi)}$.

In order to determine *E* _{vw} for any fabric, the distribution ${n(\theta\comma\; \phi)}$ must be able to represent structure that is fine enough to account for all the rheology-required structure tensors. For this reason, we expand ${n(\theta\comma\; \phi)}$ in terms of a spherical harmonic series

where $Y_{l}^{m}{[ \theta \comma\; \phi \rpar }$ are the spherical harmonic expansion functions, and ${n}^{m}_{l}$ are time-dependent expansion coefficients. Because $Y_{l}^{m}{[ \theta \comma\; \phi \rpar }$ form a complete orthonormal basis on *S* ^{2}, any distribution shape is in principal representable insofar as the expansion series is not truncated (*L* → ∞). In this sense, ${n(\theta\comma\; \phi)}$ is non-parametric since no distribution parameters exist constraining the possible range of distribution shapes. While the expansion (17) provides a more general approach for representing ${n(\theta\comma\; \phi)}$ compared to previous representations based on parametric distribution functions or low-order structure tensors (see Gagliardini and others (Reference Gagliardini, Gillel-Chaulet and Montagnat2009) and the discussion herein), it is first and foremost a mathematically convenient method that allows for higher-order structure tensors to easily be calculated: the entries of 〈**c**^{k}〉 are linear combinations of $n_l^{m}$ for *l* ≤ *k* (see Supplementary B).

Let us consider an initially isotropic fabric and truncate the expansion at *L* = 40 in order to reduce the influence of regularization on the higher-order structure tensors (see Supplementary B). Figure 4d shows the modelled fabric eigenvalues *a* _{1}, *a* _{2} and *a* _{3} of 〈**c**^{2}〉 as a function of the cumulative parcel strain ${\epsilon_{zz} = h/h_0{-1}}$, where *h* and *h* _{0} are the instantaneous and initial parcel heights, respectively. Here, we assume the usual ordering *a* _{1} ≥ *a* _{2} ≥ *a* _{3}, which implies *a* _{1} = *a* _{2} = *a* _{3} = 1/3 for isotropic fabrics, 1 > *a* _{1} > 1/3 > *a* _{2} = *a* _{3} ≥ 0 for axisymmetric preferred-direction fabrics, and 1 > *a* _{1} = *a* _{2} > 1/3 > *a* _{3} ≥ 0 for axisymmetric preferred-plane (girdle) fabrics (Woodcock, Reference Woodcock1977). For reference, Figures 4a to 4c show the specific parcel geometries and normalized ${n(\theta\comma\; \phi)}$ at ${\epsilon_{zz} = 1}$, ${\epsilon_{zz} = 0}$ and ${\epsilon_{zz} = -0.75}$, respectively.

In Figure 4e, the corresponding *E* _{mt} (dotted), *E* _{mm} (dashed) and *E* _{mt}/*E* _{pq} (full) are plotted for ${(n'\comma\;\ {{E}_{cc}'}\comma\;\ {{E}_{ca}'}) = (1\comma\;\ 1\comma\;\ 10^4) }$ in black lines and ${(n'\comma\;\ {{E}_{cc}'}\comma\;\ {{E}_{ca}'}) = (3\comma\;\ 1\comma\;\ 10^2) }$ in red lines (crosses in Figs 3a, 3b). Choosing ${{E}_{cc}'} = 1$ implies grains are equally hard to compress along **c** and **a**, which is a reasonable assumption (Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005), and ${{E}_{ca}'}$ was subsequently picked to reproduce *E* _{mt}/*E* _{pq} = 10^{4} for a unidirectional fabric (the ideal limit of ${\epsilon_{zz} = -1}$). For parcel strains of ${-0.4 \lt\epsilon_{zz} \lt 0.1}$, the linear (*n*′ = 1, or assuming (5)) and non-linear (*n*′ = 3) solutions are found to agree, while for stronger fabrics (larger absolute strains), the non-linear rheology produces bulk enhancements that scale faster with fabric strength compared to the linear rheology. Although *E* _{mt} (dotted lines) approaches the upper limits found for a unidirectional fabric as ${\epsilon_{zz}\rightarrow -1}$, the ratios *E* _{mt}/*E* _{pq} (full lines) remain far from 10^{4}. This discrepancy is the result of ${n(\theta\comma\; \phi)}$ never becoming $\delta [ \hat{\bi r}-{\bf m}\rpar$ due to non-negligible regularization for ${\epsilon_{zz}\,{\lt} -0.8}$ (see Supplementary B). In the limit *L* → ∞, no regularization is required, and we expect that ${n}{[ \theta \comma\; \phi \rpar }\rightarrow \delta [ \hat{\bi r}-{\bf m}\rpar$ as ${\epsilon_{zz}\rightarrow -1}$. Nonetheless, for the intermediate-to-strong fabrics that we *are* able to represent, we find that the bulk enhancement factors produced by a linear grain rheology (and therefore also (5)) scales more slowly with fabric strength compared to *n*′ = 3. Specifically, we find that the non-linear-to-linear ratios of *E* _{mt}/*E* _{pq} and *E* _{mm} can be *at least* 10 and 0.1, respectively.

## Discussion

A shortcoming of our study is our one-sided focus on the Sachs (homogeneous stress) assumption, which does not allow for *E* _{mt} ≃ 10 in accordance with Pimienta and others (Reference Pimienta, Duval, Lipenkov, Waddington and Walder1987). Although calculating bulk enhancement factors for both the Sachs and Taylor end-member cases would allow for a more complete assessment of the approximation (5), we note that the Sachs assumption is arguably the better of the two simple homogenization schemes (see above). Modelling efforts (Pettit and others, Reference Pettit, Thorsteinsson, Jacobson and Waddington2007; Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009; Martín and Gudmundsson, Reference Martín and Gudmundsson2012; Hruby and others, Reference Hruby2020) that take the Sachs grain-averaged rheology (9) with (5) as the bulk rheology might therefore under- or overestimate bulk strain-rate enhancements by at least an order of magnitude for intermediate-to-strong fabrics.

Other studies have adopted the orthotropic flow law proposed by Gillet-Chaulet and others (Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005), which relies on the viscoplastic self-consistent (VPSC) homogenization scheme (Meyssonnier and Philip, Reference Meyssonnier and Philip1996) to infer the bulk directional enhancement factors of the orthotropic rheology. The VPSC approach represents a compromise between the Sachs and Taylor assumptions and does not necessarily restrict *E* _{mt}, but requires an iterative numerical procedure to determine (11). Insofar as the VPSC approach is based on the grain rheology (1) with (5), such as proposed by Gillet-Chaulet and others (Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005) and adopted in several modelling studies (Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006; Durand and others, Reference Durand2007; Ma and others, Reference Ma2010), our above caveats apply, too, although to what extent is less clear. Castelnau and others (Reference Castelnau1998) showed that the VPSC approach based on a basal slip (Schmid) grain rheology with an orientation-dependent fluidity can reproduce observed strain-rate enhancements derived from mechanical tests on the GRIP ice core (single-maximum fabrics). The strongest GRIP fabric occurs at a depth of approximately 2500 m where *a* _{1} ≃ 0.95. For this fabric, the VPSC approach gave *E* _{mm} ≃ 0.03 and ${E}_{mt}^{{\dagger}ger } = 12$, compared to *E* _{mm} = 0.03 and ${E}_{mt}^{{\dagger}ger } = 65$ determined from mechanical tests, where ${E}_{mt}^{{\dagger}ger } = {E}_{mt}[ \tau _0] {\bf I}/3-{\bf p}{\bf p}\rsqb \rpar$ is the **m**–**t** shear enhancement under uniaxial compression along **p**. In comparison, for *n*′ = 3 we find *E* _{mm} = 0.05 and ${E}_{mt}^{{\dagger}ger } = 3$ (latter not shown) given ${n(\theta\comma\; \phi)}$ corresponding to *a* _{1} ≃ 0.95, although the enhancement-factor definition used by Castelnau and others (Reference Castelnau1998) is different from ours which might partly explain the discrepancy.

Given a unidirectional fabric and the Sachs assumption, Gödert and Hutter (Reference Gödert and Hutter1998) and Gagliardini and Meyssonnier (Reference Gagliardini and Meyssonnier1999) previously reported *E* _{mt} ≤ 2.5 using a linear grain rheology akin to the present. Thorsteinsson (Reference Thorsteinsson2001) considered the same configuration and reported that a basal slip (Schmid) model produces *E* _{mt} ≤ 4.375 for an orientation-dependent fluidity with a power-law exponent of three. For non-unidirectional fabrics, however, we emphasize that the rheology (1)–(2) and a Schmid model will generally not produce the same bulk directional enhancement factors (Thorsteinsson, Reference Thorsteinsson2001): in the limit of incompressibility along **c** (i.e. ${{E}_{cc}' = 0}$), longitudinal straining in the crystal basal plane is still possible with (1)–(2) (Supplementary A), unlike a Schmid model where only basal-plane shear is permitted.

We end with an outlook on the possibility of including our work in large-scale ice-flow models; i.e. by using (9)–(10) *as is*, or to determine bulk directional enhancement factors for another bulk rheology (e.g. the orthotropic flow law). The main obstacle is that ${n(\theta\comma\; \phi)}$ must be able to represent structure that is fine enough to account for all the rheology-required structure tensors. A popular fabric model is to represent ${n(\theta\comma\; \phi)}$ as a series expansion in terms of structure tensors (Gödert, Reference Gödert2003; Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005) by noting that any ${n(\theta\comma\; \phi)}$ can be expressed as a linear combination of even-ordered structure tensors (Advani and Tucker, Reference Advani and Tucker1987). In such models, calculating the time-evolution of ${n(\theta\comma\; \phi)}$ requires specifying ${\rm d}{\langle {{\bf c}^{k}}\rangle }\big /{{\rm d}t}$, but for closure reasons 〈**c**^{k}〉 for *k* ≥ 4 tend to be parameterized in terms of lower-order structure tensors. Unless 〈**c**^{k}〉 for *k* ≤ 8 are allowed to freely evolve (not parameterized), calculated bulk enhancement factors may not agree with ours. Implementing our spectral grain rotation model could therefore provide a useful path forward. To that end, we note that truncating ${n(\theta\comma\; \phi)}$ at e.g. *L* = 8, 16, or 32 implies solving a 45-, 153- or 561-dimensional linear advection–reaction–diffusion problem per grid point, respectively.

## Conclusion

We investigated how the orientation-dependent non-linear fluidity of a transversely isotropic grain rheology can influence the directional enhancement factors of a polycrystal by assuming a homogeneous stress field over the polycrystal scale (Sachs hypothesis). For polycrystals with intermediate-to-strong orientation fabrics, we find that the soft-to-hard shear-enhancement ratio, *E* _{mt}/*E* _{pq}, of axisymmetric polycrystals can be at least ten times larger compared to if the fluidity is orientation independent, while the compressional/extensional enhancement along the fabric symmetry direction, *E* _{mm}, can be at least ten times smaller. Care should therefore be exercised in large-scale ice-flow modelling efforts that either (i) identify the (Sachs) grain-averaged transversely isotropic rheology as the bulk rheology, or (ii) rely on it to derive bulk directional enhancement factors for another anisotropic bulk rheology. Our results thus extend previous work on grain-averaged rheologies with orientation-dependent fluidities (Thorsteinsson (Reference Thorsteinsson2001) and Pettit and others (Reference Pettit, Thorsteinsson, Jacobson and Waddington2007), among others) to the transversely isotropic grain rheology, which has implications for large-scale anisotropic ice-flow modelling. Finally, because the orientation-dependent fluidity causes the grain-averaged rheology to depend on higher-order fabric structure tensors, we proposed a new non-parametric spectral fabric model that allows for higher-order structure tensors to easily be calculated for arbitrary orientation fabrics.

## Acknowledgements

We thank two anonymous reviewers for comments that helped improve this manuscript greatly, and Sérgio H. Faria for his valuable comments on the spectral fabric model. The research leading to these results has received funding from the Villum Investigator grant No. 16572 as part of the IceFlow project and the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement 61 0055 as part of the ice2ice project.

## Supplementary material

The supplementary material for this article can be found at https://doi.org/10.1017/jog.2020.117.