A unified framework for large-scale fabric evolution models and anisotropic rheologies

Ice is anisotropic, with its viscosity varying by an order of magnitude in different directions when ice crystals align. However, how this variation affects ice flow is not well understood. This is because of a lack of a) models for fabric (the collective distribution of crystal orientations) evolution accurate enough to reproduce observations, and b) knowledge of which anisotropic rheology is most appropriate. Here we address both these problems. First, we review a range of previous models for fabric evolution and show they can be combined into a common differential equation. This incorporates a handful of parameters and an anisotropic rheology, which can be freely chosen. We apply this model, with a range of different anisotropic rheologies, to both an ice stream and an ice divide. For each rheology we choose the parameters to give the best possible fit to observations. We find these parameters are significantly different from those used previously. Best results come from assuming the grains rotate due to stress rather than deformation, with the stress calculated through an anisotropic rheology. By including grain rotation primarily due to stress, combined with a diffusion of the fabric, we can reproduce observations at both an ice divide and, for the first time, at an ice stream. We also compare and rank a range of anisotropic rheologies based on the accuracy of their fabric predictions. The rheologies which give the closest fit to observations have a tensor description of the anisotropy and assume that neighbouring ice grains experience approximately the same stress.

1 Introduction

The conventional approach to modelling the flow of glaciers and ice sheets is to consider ice as a continuum, viscous fluid obeying the Stokes equations (or some simplification thereof, e.g. Morland, 1987; Schoof, 2006). However, at a small enough scale ice is not a traditional fluid, but rather a polycrystalline solid constituted of grains with a regular, hexagonal, crystallographic structure (Paterson, 1999). At the grain scale, ice deforms in a very different way to liquids: while liquids deform by the free movement of molecules past one another, ice predominantly deforms through the movement of line defects along particular planes in the crystal lattice (e.g. Montagnat et al., 2014b). Because of the hexagonal lattice and molecular structure of natural ice, it is around 70 times easier for these defects to move within a particular plane – along the basal plane (Duval et al., 1983). Consequently, the collective orientation of these basal planes – the fabric – can cause the mechanical properties of ice to vary by an order of magnitude in different directions (e.g. Pimienta and Duval, 1987). This directional variation is called anisotropy.

Physically, the effect of anisotropy on ice flow can be conceptualised as follows: for a given ice velocity and stress field, the fabric will evolve as grains rotate due to the movement of dislocations within them. The new fabric then feeds back into the ice flow problem by altering locally the bulk mechanical properties of the ice. Mathematically, these physics are described by (i) an evolution equation for the fabric as a function of a given velocity field, and (ii) a constitutive relationship between stress and deformation (that is, a rheology or flow law) that depends on the fabric through a suitably defined anisotropic viscosity.

Despite the order of magnitude effect on the ease of flow, and the fact that fabrics affect flow everywhere, the number of studies including this coupled fabric-flow evolution has been limited to idealised models, either at ice divides (Martín et al., 2009; Martín and Gudmundsson, 2012; Lilien et al., 2023), simplified ice sheets (Gillet-Chaulet et al., 2006; Ma et al., 2010) or ice streams (Lilien et al., 2021; Rathmann and Lilien, 2021; Gerber et al., 2023). When considering real-world domains, the effect of fabric is usually only represented by adding to the viscosity an enhancement factor which softens the ice everywhere, for which values between 3 to 8 have been suggested (Graham et al., 2018; McCormack et al., 2022). This, rather than a coupled fabric-flow approach, is the approach taken by almost all large-scale ice sheet simulation codes whenever they account for anisotropic effects (e.g Winkelmann et al., 2011; Larour et al., 2012).

At the other end of the scale, there exist models which directly model grain dynamics or dislocation dynamics directly (see Montagnat et al., 2014b, for a review). These models can provide valuable insight into the underlying physics and some have been applied to model fabrics in ice sheets (Llorens et al., 2022). However they are too numerically expensive to be coupled to a large-scale ice flow model. In this contribution we focus solely on fabric evolution models which are computationally inexpensive enough to be included in large-scale models, which precludes simulating grain dynamics directly.

A large range of fabric measurements are available from deformation tests in the laboratory, extending from 1950s to the present day (e.g. Glen, 1952; Craw et al., 2018; Journaux et al., 2019, and references therein). However experiments need to take place at strain rates around 5 orders of magnitude greater than those seen in the natural world to allow for reasonable durations (days to weeks). Some recent work found that models for fabric evolution calibrated against laboratory experiments (Richards et al., 2021) could not reproduce fabrics observed in the natural world (Richards et al., 2023; Lilien et al., 2023), suggesting that different microstructural processes are occurring at the lower strain rates that take place in ice sheets. This partly motivates this contribution, where we seek to constrain models of anisotropy directly against observations from the natural world.

Historically there has been a lack of real-world observational data with which to compare models to. Despite recent advances in radar and seismics for measuring fabrics from the surface or air (e.g. Smith et al., 2017; Jordan et al., 2020; Kufner et al., 2023), ice cores provide the most complete dataset for fabrics. Unfortunately, ice cores have traditionally been drilled at ice divides (for a review see Faria et al., 2014a, and refs. therein), locations specifically chosen for their minimal and simple deformation in order to give accurate paleo-climate records. While we have lots of data on the fabrics produced at ice divides, there has historically been very limited data in other locations, with only a few ice cores available (Jackson and Kamb, 1997; Voigt, 2017).

Recently, the EGRIP (East Greenland Ice core project) (Stoll et al., 2021) has drilled an ice core in the much more dynamic environment of an ice stream. This represents the first complete fabric dataset for which the observed fabric is a product of a complex and changing deformation regime. Consequently, we now have full fabric data in a range of locations and deformation regimes in the natural world, which opens up new possibilities for developing and testing coupled anisotropy models.

Currently, there are two main problems in modelling large-scale anisotropy. Firstly, there is a lack of confidence in the ability of existing fabric evolution models to accurately predict fabrics in the slow (in comparison to the laboratory) and changing stress regimes seen in the natural world. This has been highlighted by the aforementioned difficulties of applying laboratory constrained models to the natural world (Richards et al., 2023; Lilien et al., 2023). Secondly, there has been no way to test whether an anisotropic rheology is accurately representing the anisotropic effect. Unlike an isotropic material, ice with a fabric will have a different mechanical response when deformed in different directions, and stresses in one direction can induce a range of deformations in other directions (what we refer to as the anisotropic effect). A laboratory experiment measures stress and deformation in the same direction giving only 1 out of the 36 components of the anisotropic viscosity tensor. This lack of information has led to a proliferation of competing rheologies (Gillet-Chaulet et al., 2005; Pettit et al., 2007; Martín et al., 2009; Placidi et al., 2010; Graham et al., 2018; Rathmann and Lilien, 2021, 2022), with no consensus on which may be most appropriate.

In this contribution, we aim to address these problems by providing both accurate modelling of ice fabrics in a range of natural conditions, and a test of the accuracy of competing rheologies. To do this, we first review and analyse the underlying physics in existing models of anisotropy, through the dynamics of an individual grain and the link between the grain scale and ice sheet scale. Physically, it is known that grains in ice are closer to the limit of experiencing uniform stress rather than uniform deformation (Castelnau et al., 1996a), and we trace this knowledge through to the large-scale models.

We develop a unified framework for fabric evolution, which can express a range of previous models of anisotropy as cases of a common equation (Sects. 2, 3, 4). This equation depends on a few parameters, and the anisotropic rheology used. This allows us to compare the anisotropic effect of different rheologies against each other, with the hypothesis that the most realistic rheology provides the most accurate fabric predictions. We apply the fabric predictions from these rheologies to two specific locations where ice cores have been drilled: an ice stream, from the aforementioned EGRIP site, and at an nearby ice divide from the earlier GRIP (Greenland Ice core Project) (Sects. 5 and 6), representing very different deformation conditions in an ice sheet.

2 Grain dynamics

The movement of glaciers and ice sheets, over centuries and millennia, is fundamentally governed by the movement of dislocations through grains (regions of similar lattice orientation) (Duval, 1981). Dislocations are line defects in the crystal lattice shape caused by impurities or a local excess or deficiency of hydrogen atoms (Glen, 1968). Under an applied stress, dislocations are generated and move in a certain direction along certain slip planes within the crystal lattice. Therefore, the deformation and thus movement of glaciers and ice sheets at the largest scales is governed by the dynamics of dislocations deforming individual grains. For ice, slip along the basal plane – with unit normal denoted by the c-axis – is approximately 70 times easier than other slip systems (Duval et al., 1983).

Table 1 List of recurring symbols used in this paper.

We start by examining a range of models that all treat the grain as the smallest physical element, i.e. all quantities are uniform within the grain. We will later make assumptions for these grain dynamics as we move to the ice sheet scale. A grain can be described by its mass m, its position x and its orientation. Due to the ease of dislocation glide along the basal plane, all the models we examine make the assumption that the basal plane is the only plane dislocations glide along. Consequently, the grains orientation can be described solely by the basal plane normal, the c-axis, c, and other crystallographic axes can be neglected. The direction of a c-axis can be described by 2 angles or, equivalently, by a unit vector in Cartesian space, and the space of possible c-axis orientations is described by the unit sphere.

2.1 Grain rheology

A constitutive law for the grain gives the deformation ${\dot{\mathit{\epsilon }}}^{\prime }$ of the grain under an applied deviatoric stress τ^′, or vice versa, where dashes denote grain quantities. The most complete approach to the grain rheology is to consider the contribution of all possible slip systems (e.g. Schmid and Casey, 1986). However, the models we discuss consider only dislocation movement along the basal plane. Therefore, with this assumption the most general model for grain rheology is to consider the grain as a transversely isotropic medium, which was first developed for metals by Johnson (1977) and first applied to ice by Meyssonnier and Philip (1996). Using the notation of Rathmann et al. (2021) the grain is described by two “enhancement factors” (Fig. 1): An enhancement factor for deformation in the basal plane, E_ca, and an enhancement factor for deformation normal to the basal plane, E_cc:

$$\begin{array}{}\text{(1)}& {\dot{\mathit{\epsilon }}}_{cc}^{\prime }={\mathit{\eta }}^{-\mathrm{1}}{E}_{cc}{\mathit{\tau }}_{cc}^{\prime },{\dot{\mathit{\epsilon }}}_{ca}^{\prime }={\mathit{\eta }}^{-\mathrm{1}}{E}_{ca}{\mathit{\tau }}_{ca}^{\prime }\end{array}$$

where the subscript c is the direction of the c-axis, and a is any direction perpendicular to this (Fig. 1), and η represents the isotropic response of the grain (i.e. some power-law). While for an isotropic material, E_cc and E_ca would both be unity, for ice E_cc≈1, and E_ca≫1, reflecting the ease of deformation along the basal plane. For a linear grain behaviour, the rheology is:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\dot{\mathit{\epsilon }}}_{ij}^{\prime }& ={\mathit{\eta }}^{-\mathrm{1}}\left({\mathit{\tau }}_{ij}^{\prime }+{\displaystyle \frac{\mathrm{3}(E_{cc}-\mathrm{1})-\mathrm{4}(E_{ca}-\mathrm{1})}{\mathrm{2}}}{\mathit{\tau }}_{kl}^{\prime }{c}_k{c}_l{c}_i{c}_j\right.\\ & \left.+(E_{ca}-\mathrm{1})({\mathit{\tau }}_{ik}^{\prime }{c}_k{c}_j+c_i{c}_k{\mathit{\tau }}_{kj}^{\prime })-{\displaystyle \frac{E_{cc}-\mathrm{1}}{\mathrm{2}}}{\mathit{\tau }}_{kl}^{\prime }{c}_k{c}_l{\mathit{\delta }}_{ij}\right)\end{array}\end{array}$$

where index notation is used and repeated indices imply summation. c is the c-axis and δ_ij is the Kronecker delta. This rheology has been used previously (Thorsteinsson, 2002; Gödert, 2003), and simplifies to the grain rheology used by Azuma and Goto-Azuma (1996) if the further assumption of keeping only terms of O(E_ca) is made. It has been extended to the non-linear case in Rathmann et al. (2021), however this introduces terms that make it computationally infeasible to include in a large-scale model. Consequently, all large-scale models which include grain dynamics have so far been limited to linear grains, and we do the same in this work. This does not preclude the use of power-law rheology at the large scale.

Figure 1 We model grains as transversely isotropic, shown on the right. This means we negelct over slip systems shown on the left, and consider only the basal plane, described by its normal c. The grain then has a “plane of isotropy” within the basal plane, where the response is the same, and the overall response can be described by an enhancement in compression E_cc≈1 and enhancement in shear E_ca≫1.

2.2 Grain rotation

Equation (2) gives a constitutive law for grain behaviour linking the deformation ${\dot{\mathit{\epsilon }}}^{\prime }$ and the stress τ^′. We must now consider how the grain's orientation c changes under an applied deformation. In reality, a grain will always be embedded inside a wider network of grains. We consider a grain experiencing a known deformation ${\dot{\mathit{\epsilon }}}^{\prime }$, as well as a known rotation, described by a rotation tensor ω^′, describing the grain's rotation relative to a fixed reference frame. The total change in rotation of a grain is a combination of the bulk rotation and the deformation d^′:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}{c}_i}{\mathrm{d}t}}=({\mathit{\omega }}_{ij}^{\prime }+d_{ij}^{\prime })c_j\end{array}$$

This rotation due to deformation is illustrated for an extensional deformation in Fig. 2. The diagonal lines represent basal planes, along which the ice deforms “like a deck of cards” (McConnel, 1891). This induces a rotation. By assuming basal planes remain parallel under simple shear (Gödert and Hutter, 1998), d^′ can be expressed as (Dafalias, 1984):

$$\begin{array}{}\text{(4)}& d_{ij}^{\prime }=c_i{c}_k{\dot{\mathit{\epsilon }}}_{kj}^{\prime }-{\dot{\mathit{\epsilon }}}_{ik}^{\prime }{c}_k{c}_j\end{array}$$

combining Eqs. (3) and (4) gives:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{c}_i}{\mathrm{d}t}}={\mathit{\omega }}_{ij}{c}_j-({\dot{\mathit{\epsilon }}}_{ij}^{\prime }{c}_j-{\dot{\mathit{\epsilon }}}_{jk}^{\prime }{c}_k{c}_j{c}_i),\end{array}$$

If a model of transversely isotropic grains is used then the expression for grain rheology (Eq. 2) can be substituted into Eq. (5) to express the lattice rotation in terms of grain stress and grain enhancements with no extra assumptions:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}{c}_i}{\mathrm{d}t}}={\mathit{\omega }}_{ij}^{\prime }{c}_j-E_{ca}{\mathit{\eta }}^{-\mathrm{1}}({\mathit{\tau }}_{ij}^{\prime }{c}_j-{\mathit{\tau }}_{jk}^{\prime }{c}_k{c}_j{c}_i).\end{array}$$

The derivation of this can be found in Appendix A.

Figure 2 Illustration of an ice grain deforming in tension, adapted from Placidi (2005). (a) the initial undeformed grain, with diagonal lines representing basal planes, and c-axis shown. (b) slide alone the basal planes deforms the grain under tension, however here the crystal is not in equilibrium as there is a couple acting on it. (c) the equilibrium position involves grain rotation as well, giving a change in the c-axis orientation.

2.3 Recrystallization

2.3.1 Migration recrystallization

We now consider how the mass of a grain m evolves. While grains in the firn or those undergoing no deformation may undergo normal grain growth (Arthern et al., 2010), in an ice sheet which is always deforming, the primary process for grain growth will be migration recrystallization (Faria et al., 2014b).

Placidi et al. (2010) considered the effect of migration recrystallization on a continuum description of the fabric. In this sub-section, we aim to more explicitly explore how the equation of Placidi et al. (2010) can be derived by considering the dynamics of an individual grain.

Migration recrystallization is the process where one grain grows at the expense of another. Hence, whereas previously we considered a single grain we must now consider a system of N grains. Migration recrystallization is driven by the minimisation of potential energy of the system. Grains with a higher dislocation density have higher potential energy. Consequently, the change in mass of a grain $\dot{m}$ is a function of this dislocation density ρ_D and the grains mass m:

$$\begin{array}{}\text{(7)}& \dot{m}=f({\mathit{\rho }}_{\mathrm{D}},m),{\displaystyle \frac{\partial \dot{m}}{\partial {\mathit{\rho }}_{\mathrm{D}}}}<\mathrm{0}.\end{array}$$

where f is some arbritary function. The total mass of the system must also be conserved:

$$\begin{array}{}\text{(8)}& \sum _{i=\mathrm{1}}^N{\dot{m}}_i=\sum _{i=\mathrm{1}}^{N}f({{\mathit{\rho }}_{\mathrm{D}}}_i,m_i)=\mathrm{0}\end{array}$$

Dislocations exist at the molecular scale so for any fabric evolution model aiming to be applied at the large scale they cannot be modelled directly. The dislocation density, and hence migration recrystallization, is also likely related to the history of the grain. However, the state-of-the-art for modelling migration recrystallization at the large scale assumes it is a function of the instantaneous properties (Placidi et al., 2010). To arrive at this result we introduce a proxy for dislocation density ρ_D such that is a function of the grain properties described previously (i.e. the stress τ^′, deformation ${\dot{\mathit{\epsilon }}}^{\prime }$, mass m and orientation c). Humphreys and Hatherly (2004) propose that the dislocation density is related to the direction of the stress acting on the grain, i.e. ${\mathit{\rho }}_{\mathrm{D}}=f({\mathit{\tau }}^{\prime },\mathit{c})$ (where f is again some arbitrary function, not the same as Eq. 7). For example, if the stress imposed on the grain is in the direction of the basal plane, the dislocations can move easily to deform the grain. However, if the stress imposed on the grain is perpendicular to the basal plane, any existing dislocations cannot move easily. Consequently, new dislocations will nucleate to accommodate the applied stress, increasing the dislocation density. We can express this mathematically as:

$$\begin{array}{}\text{(9)}& {\mathit{\rho }}_{\mathrm{D}}=f\left(s^{\prime }\right),{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{D}}}{\partial s^{\prime }}}<\mathrm{0}\end{array}$$

where s^′ is the resolved stress on the basal plane, and the dislocation density decreases as the resolved stress increases (i.e. the stress is in the direction of the basal planes). The resolved stress can be expressed as (Placidi et al., 2010):

$$\begin{array}{}\text{(10)}& s^{\prime }={\displaystyle \frac{\left({\mathit{\tau }}_{ij}^{\prime }{c}_j\right)\left({\mathit{\tau }}_{ik}^{\prime }{c}_k\right)-({\mathit{\tau }}_{ij}^{\prime }{c}_j{c}_i{)}^{\mathrm{2}}}{{\mathit{\tau }}_{mn}^{\prime }{\mathit{\tau }}_{nm}^{\prime }}},\end{array}$$

Combining Eqs. (7) and (9) gives:

$$\begin{array}{}\text{(11)}& \dot{m}=f(s^{\prime },m),{\displaystyle \frac{\partial \dot{m}}{\partial s^{\prime }}}>\mathrm{0}\end{array}$$

Based on Placidi et al. (2010), we then make the simplest possible assumption, that the functional relationship between s^′ and m can be approximated by multiplying the two variables together:

$$\begin{array}{}\text{(12)}& \dot{m}=\mathit{\beta }m{s}^{\prime }-C\end{array}$$

where β is a constant of proportionality and C can be found through conservation of mass, Eq. (11) becomes:

$$\begin{array}{}\text{(13)}& {\dot{m}}_i=\mathit{\beta }{m}_i\left(s_i^{\prime }-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N\left(m_i{s}_i^{\prime }\right)}{{\sum }_{i=\mathrm{1}}^N{m}_i}}\right)\end{array}$$

Where the subscript i denotes an individual grain. To simplify this equation, we can replace m with a mass fraction w:

$$\begin{array}{}\text{(14)}& w_i={\displaystyle \frac{m_i}{{\sum }_{i=\mathrm{1}}^N{m}_i}}\end{array}$$

to finally give an equation instead for the change in mass fraction $\mathrm{d}w/\mathrm{d}t$:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathrm{d}w}{\mathrm{d}t}}=\mathit{\beta }w(s^{\prime }-\langle s^{\prime }\rangle )\end{array}$$

where 〈s^′〉 denotes a weighted average over N grains ${\sum }_{i=\mathrm{1}}^N{w}_i{s}_i^{\prime }$. As we show in Sect. 3, this is equivalent to the term derived by Placidi et al. (2010) for migration recrystallization, with the exception that s^′ in Placidi et al. (2010) is a function of the large-scale deformation. The parameter β has been shown to be strongly dependent on temperature (Richards et al., 2021) and, to a lesser extent, strain rate.

2.3.2 Rotational recrystallization

We are reviewing models which assume a grain has a uniform stress τ^′ and deformation ${\dot{\mathit{\epsilon }}}^{\prime }$. However, this is a simplification. In reality, it has been observed that as a grain is deformed it tends to have a less stressed inner “core” and a more stressed outer “mantle” (Faria et al., 2009). Under stress, subgrains – small changes in orientation within the grain – form and grow in this outer region as a result of increased dislocation density and recovery (Drury and Urai, 1990). As strain increases, eventually a subgrain will have a sufficiently different (∼10°) orientation from the core grain that it can be called a new grain. This process is called rotational recrystallization. The subgrains orientation will depend on the orientation of its parent grain, the stress on the grain and the orientation of its neighbours (e.g. Piazolo et al., 2015). At the large scale, this acts to diffuse concentrations in the fabric.

However the models we are reviewing, which aim to be included in large-scale ice sheet models, all neglect heterogenities within a grain as well as the interactions of neighbouring grains for reasons of computational efficiency. Gödert (2003) first proposed incorporating the effect of rotational recrystallization by incorporating an additional Brownian motion term alongside the lattice rotation term $\mathrm{d}\mathit{c}/\mathrm{d}t$ on the change in grain orientation away from the parent grain:

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\mathrm{d}{c}_i}{\mathrm{d}t}}=v_i+\sqrt{\mathrm{2}\mathit{\lambda }}{W}_i,\end{array}$$

where v represents the contribution to $\mathrm{d}\mathit{c}/\mathrm{d}t$ from basal-slip deformation such as the right hand side in Eq. (6), λ is a parameter and W(t) is a multivariate Wiener process describing Brownian motion. This can be thought of as subgrains causing the average orientation of the grain as a whole to change in a “random” way. This Browninan motion corresponds, as in molecular dynamics, to a diffusion at a larger scale – in this case of fabric concentrations. While the change in orientation would not be random if grain interactions could be modelled, observations from ice cores Durand et al. (2008) suggest it appears random when compared to the large scale deformation.

3 Moving from a grain to an ice sheet

So far, we have reviewed models for an isolated grain, internally uniform and subjected to externally imposed stresses or deformations. We collected equations for the rheology (Meyssonnier and Philip, 1996; Rathmann et al., 2021), change in c-axis orientation due to deformation (Azuma and Goto-Azuma, 1996) or rotational recrystallization (Gödert, 2003), and change in grain mass (Placidi et al., 2010) into a common model for an ice grain:

$$\begin{array}{}\text{(17)}& {\displaystyle }\begin{array}{rl}{\dot{\mathit{\epsilon }}}_{ij}^{\prime }& ={\mathit{\eta }}^{-\mathrm{1}}\left({\mathit{\tau }}_{ij}^{\prime }+{\displaystyle \frac{\mathrm{3}(E_{cc}-\mathrm{1})-\mathrm{4}(E_{ca}-\mathrm{1})}{\mathrm{2}}}{c}_i{c}_j{c}_k{c}_l{\mathit{\tau }}_{kl}^{\prime }\right.\\ & \left.+(E_{ca}-\mathrm{1})({\mathit{\tau }}_{ik}^{\prime }{c}_k{c}_j+c_i{c}_k{\mathit{\tau }}_{kj}^{\prime })-{\displaystyle \frac{E_{cc}-\mathrm{1}}{\mathrm{2}}}{c}_k{c}_l{\mathit{\tau }}_{kl}^{\prime }{\mathit{\delta }}_{ij}\right),\end{array}\text{(18)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{c}_i}{\mathrm{d}t}}=v_i+\sqrt{\mathrm{2}\mathit{\lambda }}{W}_i,v_i& ={\mathit{\omega }}_{ij}^{\prime }-({\dot{\mathit{\epsilon }}}_{ij}^{\prime }{c}_j-{\dot{\mathit{\epsilon }}}_{jk}^{\prime }{c}_k{c}_j{c}_i)\\ & ={\mathit{\omega }}_{ij}^{\prime }-E_{ca}{\mathit{\eta }}^{-\mathrm{1}}({\mathit{\tau }}_{ij}^{\prime }{c}_j-{\mathit{\tau }}_{jk}^{\prime }{c}_k{c}_j{c}_i)\end{array}\text{(19)}& {\displaystyle \frac{\mathrm{d}w}{\mathrm{d}t}}=\mathit{\beta }w(s^{\prime }-\langle s^{\prime }\rangle ),s^{\prime }={\displaystyle \frac{\left({\mathit{\tau }}_{ij}^{\prime }{c}_j\right)\left({\mathit{\tau }}_{ik}^{\prime }{c}_k\right)-({\mathit{\tau }}_{ij}^{\prime }{c}_j{c}_i{)}^{\mathrm{2}}}{{\mathit{\tau }}_{mn}^{\prime }{\mathit{\tau }}_{nm}^{\prime }}}\end{array}$$

The above equations could be solved by placing the grains in a network in space, and deriving more equations for their interactions in a similar fashion to micro-mechanical models (e.g. Castelnau et al., 1996a; Piazolo et al., 2010).

However, as discussed previously, if we are concerned with an ice sheet or glacier, this is numerically infeasible as there are around 10¹² grains in each cubic kilometre. Instead, what is done is to make a continuum assumption, that at each point in the material there are a large number of grains. This is analogous to considering a glass of water as a continuum fluid rather than a collection of water molecules. Like the velocity of a fluid at a point is the mean velocity of all the molecules, so small that they can be considered all “at that point”, in our case most quantities can be found by taking a weighted average of all the grains at that point, e.g. for the deformation tensor:

$$\begin{array}{}\text{(20)}& \dot{\mathit{\epsilon }}=\langle {\dot{\mathit{\epsilon }}}^{\prime }\rangle =\sum _{i=\mathrm{1}}^N{w}_i{\dot{\mathit{\epsilon }}}_i\end{array}$$

However, continuing the analogy with a glass of water, when this continuum approximation is made emergent quantities occur: for example, the temperature of an individual water molecule is not a relevant concept but the temperature of the glass of water as a whole is a key control of its physics. The emergent quantity in this case is the fabric. The fabric is the distribution of grain orientations in a region. It can be represented by a probability distribution function, often called the orientation distribution function f. This is different from a usual distribution function in that it is defined over the space of possible orientations, i.e. over two angles prescribing the surface of a sphere. In the real world, fabrics are measured when ice cores are drilled and the orientation of each grain within the core is measured. These can be considered as samples from the distribution function. If the dynamics of one sample are described by Eqs. (18) and (19), the Feynman-Kac formula (Kac, 1949) gives the corresponding partial differential equation for the probability distribution function f:

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\mathrm{d}f\left(\mathit{n}\right)}{\mathrm{d}t}}=-{\mathrm{\nabla }}_i^{*}\left(f{v}_i\right)+\mathit{\lambda }{{\mathrm{\nabla }}^{*}}^{\mathrm{2}}\left(f\right)+\mathit{\beta }f(s^{\prime }-\langle s^{\prime }\rangle ),\end{array}$$

here n is the space of all possible orientations a single c-axis may have, and v and s^′ are functions over n and ∇^* is the gradient operator over n (not physical space x) restricted to the space of possible orientations:

$$\begin{array}{}\text{(22)}& {\mathrm{\nabla }}_i^{*}={\displaystyle \frac{\partial }{\partial n_i}}-{\displaystyle \frac{\partial }{\partial n_j}}{n}_j{n}_i.\end{array}$$

This equation has been used before a number of times, with differing forms for v, and sometimes without migration recrystallization (Gödert, 2003; Faria, 2006; Placidi et al., 2010; Richards et al., 2021; Rathmann et al., 2021).

The problem we have with this equation is that v and s^′ depend on grain quantities. However when viscous anisotropy – controlled by fabric – is included in ice sheet models, it is included by coupling the equations for anisotropy to a large scale flow. At each timestep, the input to the anisotropic model is the flow field, from which the deformation tensor $\dot{\mathit{\epsilon }}$ can be trivially calculated. The output from the model of anisotropy will be some viscosity, described by a fourth-rank tensor, such that the equation for the flow can be solved for the next timestep.

If we only have macroscopic quantities as inputs, we need to recast Eqs. (17)–(19) in terms of these macroscopic deformations and stresses rather than grain quantities. Equation (21) can then expressed in macroscopic quantities only.

3.1 The Taylor and Sachs approximations

To express Eqs. (18) and (19), and hence Eq. (21) in macroscopic quantities, there are two end members available: either all grains deform at the same rate ${\dot{\mathit{\epsilon }}}^{\prime }=\langle {\dot{\mathit{\epsilon }}}^{\prime }\rangle =\dot{\mathit{\epsilon }}$, or all grains experience the same stress ${\mathit{\tau }}^{\prime }=\langle {\mathit{\tau }}^{\prime }\rangle =\mathit{\tau }$ . These two approximations, which act as the two bounds that the true grain behaviour will lie between, are called the Taylor and Sachs bounds respectively. Polycrystalline materials which deform close to the Taylor bound have many easy slip systems, such as some metals (e.g Thamburaja and Jamshidian, 2014). This means that dislocations can move in different directions, and the applied deformation can be easily accommodated by the grains. However, ice has only 1 easy slip system: the basal plane. When applied to ice, the Taylor approximation predicts around 60 % of the strain is due to non basal slip, and weak fabrics (Castelnau et al., 1996b, a). Observations from ice sheets which predominantly show strong fabrics, and the fact that slip along the basal plane is around 70 times easier than other slip systems (Duval et al., 1983) provide evidence that ice is deforming closer to the Sachs bound than the Taylor.

3.1.1 Sachs model

Consequently, the Sachs approximation ${\mathit{\tau }}^{\prime }=\langle {\mathit{\tau }}^{\prime }\rangle =\mathit{\tau }$ more closely represents reality for ice, and is used by Azuma and Goto-Azuma (1996) and Gödert (2003) to derive a macroscopic rheology. The Sachs rheology can be derived by applying the averaging operater 〈⋅〉 to Eq. (17):

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{\dot{\mathit{\epsilon }}}_{ij}& ={\mathit{\eta }}^{-\mathrm{1}}\left({\mathit{\tau }}_{ij}+{\displaystyle \frac{\mathrm{3}(E_{cc}-\mathrm{1})-\mathrm{4}(E_{ca}-\mathrm{1})}{\mathrm{2}}}{A}_{ijkl}^{\left(\mathrm{4}\right)}{\mathit{\tau }}_{kl}\right.\\ & \left.+(E_{ca}-\mathrm{1})({\mathit{\tau }}_{ik}{A}_{kj}+A_{ik}{\mathit{\tau }}_{kj})-{\displaystyle \frac{E_{cc}-\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_{ij}{A}_{kl}{\mathit{\tau }}_{kl}\right)\end{array}\end{array}$$

where A and A⁽⁴⁾ are defined below in Eqs. (26) and (27). For linear grain behaviour, ${\mathit{\eta }}^{-\mathrm{1}}=B$, the commonly used rate factor (which we denote with B to avoid confusion with the orientation tensors). However Eq. (23) has been extended to include a power-law response by Pettit et al. (2007) by setting:

$$\begin{array}{}\text{(24)}& {\mathit{\eta }}^{-\mathrm{1}}=B({\mathit{\tau }}_{ij}{\mathit{\tau }}_{ji}{)}^{(n-\mathrm{1})/\mathrm{2}}\end{array}$$

and by Martín et al. (2009) by setting:

$$\begin{array}{}\text{(25)}& \mathit{\eta }=B^{-\mathrm{1}/n}({\dot{\mathit{\epsilon }}}_{ij}{\dot{\mathit{\epsilon }}}_{ji}{)}^{(\mathrm{1}-n)/\mathrm{2}n}\end{array}$$

allowing a macroscopic power law response through n even if we assume the grains behave with a linear response. In Eq. (23) A is the second-order orientation tensor:

$$\begin{array}{}\text{(26)}& A_{ij}=\langle c_i{c}_j\rangle =\sum _{g=\mathrm{1}}^N{w}_g{c_g}_i{c_g}_j=\underset{S^{\mathrm{2}}}{\int }f\left(\mathit{n}\right)n_i{n}_j\mathrm{d}\mathit{n}\end{array}$$

As can be seen from this equation, it is also the second moment of the orientation distribution function f, and is a second-rank tensor. The eigenvalues of A are an important measure of the fabric. Similarly, the fourth-order orientation tensor is:

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{A}_{ijkl}^{\left(\mathrm{4}\right)}& =\langle c_i{c}_j{c}_k{c}_l\rangle =\sum _{g=\mathrm{1}}^N{w}_g{c_g}_i{c_g}_j{c_g}_k{c_g}_l\\ & =\underset{S^{\mathrm{2}}}{\int }f\left(\mathit{n}\right)n_i{n}_j{n}_k{n}_l\mathrm{d}\mathit{n}\end{array}\end{array}$$

Note for the case E_cc=1 and E_ca=1, Eq. (23) simplifies to Glen's rheology.

While the Sachs rheology has been used before, in this contribution we use it in fabric predictions and compare it to other rheologies. To do this, its helpful to normalise each rheology by the magnitude of its isotropic response. For the Sachs rheology, the isotropic response can be found by substituting orientation tensors corresponding to an isotropic fabric ($A_{ij}={\mathit{\delta }}_{ij}/\mathrm{3}$ and $A_{ijkl}=({\mathit{\delta }}_{ij}{\mathit{\delta }}_{kl}+{\mathit{\delta }}_{ik}{\mathit{\delta }}_{jl}+{\mathit{\delta }}_{il}{\mathit{\delta }}_{jk})/\mathrm{15}$) into Eq. (23), giving:

$$\begin{array}{}\text{(28)}& {\dot{\mathit{\epsilon }}}_{ij}={\mathit{\eta }}^{-\mathrm{1}}(\mathrm{0.4}{E}_{ca}+\mathrm{0.2}{E}_{cc}+\mathrm{0.4}){\mathit{\tau }}_{ij}\end{array}$$

We introduce normalised deviatoric stress $\widehat{\mathit{\tau }}$, which is defined as being equal to the deformation tensor $\dot{\mathit{\epsilon }}$ for an isotropic fabric. This means that, if we start with the deformation tensor, the normalised stress captures only the changes due to fabric anisotropy, and does not include any isotropic softening. For this rheology,

$$\begin{array}{}\text{(29)}& {\widehat{\mathit{\tau }}}_{ij}={\mathit{\eta }}^{-\mathrm{1}}(\mathrm{0.4}{E}_{ca}+\mathrm{0.2}{E}_{cc}+\mathrm{0.4}){\mathit{\tau }}_{ij}\end{array}$$

and we can rewrite Eq. (23) in terms of the normalised stress:

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{\dot{\mathit{\epsilon }}}_{ij}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{0.4}{E}_{ca}+\mathrm{0.2}{E}_{cc}+\mathrm{0.4}}}\\ & \left({\widehat{\mathit{\tau }}}_{ij}+{\displaystyle \frac{\mathrm{3}(E_{cc}-\mathrm{1})-\mathrm{4}(E_{ca}-\mathrm{1})}{\mathrm{2}}}{A}_{ijkl}^{\left(\mathrm{4}\right)}{\widehat{\mathit{\tau }}}_{kl}\right.\\ & \left.+(E_{ca}-\mathrm{1})({\widehat{\mathit{\tau }}}_{ik}{A}_{kj}+A_{ik}{\widehat{\mathit{\tau }}}_{kj})-{\displaystyle \frac{E_{cc}-\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_{ij}{A}_{kl}{\widehat{\mathit{\tau }}}_{kl}\right)\end{array}\end{array}$$

This is a linear equation and if $\dot{\mathit{\epsilon }}$ is known, it can be inverted easily to find $\widehat{\mathit{\tau }}$. Note that Eq. (30) does not depend on η⁻¹, and hence it does not depend on n.

Equations (18) and (19) can also be expressed in macroscopic quantites using the Sachs approximation. It is also necessary to approximate the grain rotation ω^′. For this, the approximation ${\mathit{\omega }}^{\prime }=\mathit{\omega }=(\mathrm{\nabla }\mathit{u}-\mathrm{\nabla }{\mathit{u}}^T)/\mathrm{2}$, i.e. the grain rotation is equal to the bulk rotation, is made as no other information is available. Eq. (18) can be combined with Eq. (6) to express v in macroscopic quantities:

$$\begin{array}{}\text{(31)}& v_i={\mathit{\omega }}_{ij}{c}_j-\mathit{\iota }({\widehat{\mathit{\tau }}}_{ij}{c}_j-{\widehat{\mathit{\tau }}}_{jk}{c}_k{c}_j{c}_i)\end{array}$$

where

$$\begin{array}{}\text{(32)}& \mathit{\iota }={\displaystyle \frac{E_{ca}}{\mathrm{0.4}{E}_{ca}+\mathrm{0.2}{E}_{cc}+\mathrm{0.4}}}.\end{array}$$

Note that Eq. (31) now depends asymptotically on E_ca, as E_ca→∞ ι→2.5.

Similarly the expression for s^′ in Eq. (19) can be expressed in macroscopic quantities:

$$\begin{array}{}\text{(33)}& s^{\prime }={\displaystyle \frac{{\widehat{\mathit{\tau }}}_{ij}{c}_j{\widehat{\mathit{\tau }}}_{ik}{c}_k-({\widehat{\mathit{\tau }}}_{ij}{c}_j{c}_i{)}^{\mathrm{2}}}{{\widehat{\mathit{\tau }}}_{mn}{\widehat{\mathit{\tau }}}_{nm}}}\end{array}$$

Because both Eqs. (31) and (33) can be written solely depending on $\widehat{\mathit{\tau }}$ rather than τ, the fabric evolution predictions do not depend on the value of n.

3.1.2 Taylor model

If the Taylor approximation is used (${\dot{\mathit{\epsilon }}}^{\prime }=\langle {\dot{\mathit{\epsilon }}}^{\prime }\rangle =\dot{\mathit{\epsilon }}$) then v from Eq. (18) can be expressed in macroscopic quantities:

$$\begin{array}{}\text{(34)}& v_i={\mathit{\omega }}_{ij}{c}_j-({\dot{\mathit{\epsilon }}}_{ij}{c}_j-{\dot{\mathit{\epsilon }}}_{jk}{c}_k{c}_j{c}_i)\end{array}$$

A macroscopic equation for s^′ is more complicated with the Taylor approximation and requires a constitutive law. However, as we test the model in only low temperature environments, where migration recrystallization is insignificant such that we set β=0 (see Sect. 6.1) we have not derived this term.

4 Macroscopic models

In the previous section we have reviewed two models for anisotropy based on assumptions about grain dynamics: the Taylor and Sachs model. In particular, we have shown that for the Sachs model the fabric evolution equation Eq. (31) depends on the stress, which must be found through a rheology if we start with the velocity hence deformation $\dot{\mathit{\epsilon }}$.

We aim to test the set anisotropic rheologies aiming to be included in ice sheet models, through their effect on fabric predictions. Consequently, we now review these other models. Instead of considering the dynamics of grains directly, these models consider equations directly at the macroscopic scale.

4.1 Orthotropic rheologies

Gillet-Chaulet et al. (2005) and Rathmann and Lilien (2022) derive anisotropic rheologies not by assuming an ice grain has a certain stress or deformation, but by assuming that the rheology has certain symmetries, specifically that it is invariant under reflections along symmetry directions of the fabric. As these rheologies work from assumptions about the fabric rather than the grain, they can be called macroscopic.

4.1.1 General Orthotropic Linear flow law (GOLF)

The GOLF rheology (Gillet-Chaulet et al., 2005) can be expressed as:

$$\begin{array}{}\text{(35)}& \begin{array}{rl}{\widehat{\mathit{\tau }}}_{ij}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{r=\mathrm{1}}^{\mathrm{3}}\left[{\mathit{\lambda }}_r\left({\dot{\mathit{\epsilon }}}_{kl}{m}_l^r{m}_k^r\right)(m_i^r{m}_j^r{)}^{\mathrm{D}}\right.\\ & \left.+{\mathit{\lambda }}_{r+\mathrm{3}}(m_i^r{m}_k^r{\dot{\mathit{\epsilon }}}_{kj}+{\dot{\mathit{\epsilon }}}_{ik}{m}_k^r{m}_j^r{)}^{\mathrm{D}}\right],\end{array}\end{array}$$

where $\widehat{\mathit{\tau }}$ is the normalised stress and ${\mathbf{X}}^{\mathrm{D}}=\mathbf{X}-\mathrm{tr}\left(\mathbf{X}\right)\mathbf{I}/\mathrm{3}$ represents the deviatoric part of a tensor. The vectors ${\mathit{m}}^{\mathrm{1}},{\mathit{m}}^{\mathrm{2}},{\mathit{m}}^{\mathrm{3}}$, are the symmetry directions of the fabric, which are assumed to be along the eigenvectors of the second-order orientation tensor A. The values of the terms λ_r and λ_r+3 are drawn tabulated values, set to reproduce to behaviour of a micro-mechanical model including grain interactions. In this way the GOLF aims to incorporate grain physics between the Sachs and Taylor bounds and theoretically be more accurate. The micro-mechanical model the GOLF is constrained against (Castelnau et al., 1996a) depends on the grain parameters E_cc and E_ca¹. While, as the name suggests, the GOLF was initially a linear rheology, it has been coupled to the power law response in Eq. (24) in Pettit et al. (2007); Ma et al. (2010), although this does not affect the normalised stress.

4.1.2 Non-Linear Orthotropic (NLO) rheology

Similarly, Rathmann and Lilien (2022) derive an alternative orthotropic rheology:

$$\begin{array}{}\text{(36)}& {\dot{\mathit{\epsilon }}}_{ij}={\mathit{\eta }}^{-\mathrm{1}}\sum _{r=\mathrm{1}}^{\mathrm{3}}\left({\mathit{\lambda }}_r{\mathit{\tau }}_{kl}{M}_{kl}^r{M}_{ij}^r+{\mathit{\lambda }}_{r+\mathrm{3}}{\mathit{\tau }}_{kl}{M}_{kl}^{r+\mathrm{3}}{M}_{ij}^{r+\mathrm{3}}\right)\end{array}$$

where $M_{ij}^{\mathit{\{}r=\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}}}=(m_i^{\mathit{\{}\mathrm{2},\mathrm{3},\mathrm{1}\mathit{\}}}{m}_j^{\mathit{\{}\mathrm{2},\mathrm{3},\mathrm{1}\mathit{\}}}-m_i^{\mathit{\{}\mathrm{3},\mathrm{1},\mathrm{2}\mathit{\}}}{m}_j^{\mathit{\{}\mathrm{3},\mathrm{1},\mathrm{2}\mathit{\}}})/\mathrm{2}$ and $M_{r+\mathrm{3}}^{\mathit{\{}r=\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}}}=(m_i^{\mathit{\{}\mathrm{2},\mathrm{3},\mathrm{1}\mathit{\}}}{m}_j^{\mathit{\{}\mathrm{3},\mathrm{1},\mathrm{2}\mathit{\}}}+m_i^{\mathit{\{}\mathrm{3},\mathrm{1},\mathrm{2}\mathit{\}}}{m}_j^{\mathit{\{}\mathrm{2},\mathrm{3},\mathrm{1}\mathit{\}}})/\mathrm{2}$. For this rheology, the viscosities λ_r and λ_r+3 are constrained against a linear combination of the linear Sachs and Taylor rheologies. They use a ratio of 98.75:1.25 Sachs to Taylor weighting, so the rheology closely matches the linear Sachs rheology, but gives a larger maximum anisotropic effect for very strong fabrics.

Rathmann and Lilien (2022) draw a distinction between the grain power law response and the macroscopic power law response. The terms λ_i and λ_i+3 depend on the grain parameters E_cc and E_ca, the Sachs-Taylor weight, but also on the macroscopic power exponent n. They also provide an unapproximated form of η⁻¹ in the context of the orthotropic rheology:

$$\begin{array}{}\text{(37)}& \begin{array}{rl}{\mathit{\eta }}^{-\mathrm{1}}& =B\left(\sum _{r=\mathrm{1}}^{\mathrm{3}}\left[{\mathit{\lambda }}_r({\mathit{\tau }}_{ij}:M_{ji}^r{)}^{\mathrm{2}}\right.\right.\\ & {\left.\left.+{\mathit{\lambda }}_{r+\mathrm{3}}({\mathit{\tau }}_{ij}:M_{ji}^{r+\mathrm{3}}{)}^{\mathrm{2}}\right]\right)}^{(n-\mathrm{1})/\mathrm{2}}\end{array}\end{array}$$

The normalised stress is given by:

$$\begin{array}{}\text{(38)}& {\widehat{\mathit{\tau }}}_{ij}=B({\mathit{\tau }}_{kl}{\mathit{\tau }}_{lk}{)}^{(n-\mathrm{1})/\mathrm{2}}{\mathit{\tau }}_{ij}\end{array}$$

4.2 Other rheologies

We now review other rheologies which have been used to attempt to approximate the anisotropic effect.

4.2.1 CAFFE rheology

Another macroscopic rheology is the CAFFE rheology (Placidi et al., 2010). There has been some discussion in the literature over whether this rheology is truly anisotropic (Faria, 2008), as it has a co-linear relationship between stress and deformation:

$$\begin{array}{}\text{(39)}& {\dot{\mathit{\epsilon }}}_{ij}=BE\left(\mathcal{D}\right)({\mathit{\tau }}_{kl}{\mathit{\tau }}_{lk}{)}^{(n-\mathrm{1})/\mathrm{2}}{\mathit{\tau }}_{ij}\end{array}$$

where

$$\begin{array}{}\text{(40)}& \mathcal{D}=\mathrm{5}{\displaystyle \frac{A_{ij}{\dot{\mathit{\epsilon }}}_{ik}{\dot{\mathit{\epsilon }}}_{kj}-A_{ijkl}^{\left(\mathrm{4}\right)}{\dot{\mathit{\epsilon }}}_{kl}{\dot{\mathit{\epsilon }}}_{ij}}{{\dot{\mathit{\epsilon }}}_{mn}{\dot{\mathit{\epsilon }}}_{nm}}},\end{array}$$

and E is a monotonic function of 𝒟:

$$\begin{array}{}\text{(41)}& \begin{array}{rl}& E\left(\mathcal{D}\right)=\\ & \left\{\begin{array}{ll}(\mathrm{1}-E_{\text{min}}){\mathcal{D}}^k+E_{\text{min}},k=\frac{\mathrm{8}}{\mathrm{21}}\left(\frac{E_{\text{max}}-\mathrm{1}}{\mathrm{1}-E_{\text{min}}}\right),& \mathcal{D}\in [\mathrm{0},\mathrm{1}]\\ \frac{\mathrm{4}{\mathcal{D}}^{\mathrm{2}}(E_{\text{max}}-\mathrm{1})+\mathrm{25}-\mathrm{4}{E}_{\text{max}}}{\mathrm{21}}& \mathcal{D}\in [\mathrm{1},\frac{\mathrm{5}}{\mathrm{2}}].\end{array}\right.\end{array}\end{array}$$

For an isotropic fabric, 𝒟=1 and E(𝒟)=1, whereas for a perfect single-maximum fabric $\mathcal{D}=\mathrm{5}/\mathrm{2}$ and E(𝒟)=E_max. Therefore, expressing the rheology in terms of the normalised stress $\widehat{\mathit{\tau }}$, which is defined as $\widehat{\mathit{\tau }}=\dot{\mathit{\epsilon }}$ for an isotropic fabric only, gives:

$$\begin{array}{}\text{(42)}& {\dot{\mathit{\epsilon }}}_{ij}=E\left(\mathcal{D}\right){\widehat{\mathit{\tau }}}_{ij}\end{array}$$

i.e. E cannot be incorporated into the isotropic response. Faria (2008) shows because of this Eq. (39) is technically anisotropic, even though the relationship is co-linear.

4.2.2 Isotropic rheologies: Glen's and E^*

There have also been attempts to represent anisotropy through considering the velocity and deformation only, such as the Estar flow relation (Graham et al., 2018). This rheology can be expressed as

$$\begin{array}{}\text{(43)}& {\dot{\mathit{\epsilon }}}_{ij}={\mathit{\eta }}^{-\mathrm{1}}(\dot{\mathit{\epsilon }},\mathit{u}){\mathit{\tau }}_{ij}.\end{array}$$

For any rheology like this, η is not a function of the fabric. Therefore $\dot{\mathit{\epsilon }}=\widehat{\mathit{\tau }}$ for any fabric, and the Estar rheology will give the same fabric predictions as Glen's rheology.

4.3 A common equation for fabric evolution

As the above macroscopic models are not tied to either the Sachs or Taylor bound, and often aim to give a response between them, they can postulate different equations for fabric evolution. Azuma (1994) first postulated a linear blend between the Taylor and Sachs approximations:

$$\begin{array}{}\text{(44)}& \begin{array}{rl}{v}_i& ={\mathit{\omega }}_{ij}{c}_j-(\mathrm{1}-\mathit{\alpha })({\dot{\mathit{\epsilon }}}_{ij}{c}_j-{\dot{\mathit{\epsilon }}}_{jk}{c}_k{c}_j{c}_i)\\ & -\mathit{\alpha }{E}_{ca}{\mathit{\eta }}^{-\mathrm{1}}({\mathit{\tau }}_{ij}{c}_j-{\mathit{\tau }}_{jk}{c}_k{c}_j{c}_i),\end{array}\end{array}$$

where $\mathit{\alpha }\in [\mathrm{0},\mathrm{1}]$, with α=0 corresponding to the Taylor bound and α=1 corresponding to the Sachs bound. Although Eq. (44) seems to suggest that as E_ca increase the contribution to lattice rotation from the stress increases. However, if the Sachs approximation and rheology is used, Eq. (31) shows that rewriting Eq. (44) in terms of the normalised stress gives:

$$\begin{array}{}\text{(45)}& \begin{array}{rl}{v}_i& ={\mathit{\omega }}_{ij}{c}_j-(\mathrm{1}-\mathit{\alpha })({\dot{\mathit{\epsilon }}}_{ij}{c}_j-{\dot{\mathit{\epsilon }}}_{jk}{c}_k{c}_j{c}_i)\\ & -\mathit{\alpha }\mathit{\iota }({\widehat{\mathit{\tau }}}_{ij}{c}_j-{\widehat{\mathit{\tau }}}_{jk}{c}_k{c}_j{c}_i),\end{array}\end{array}$$

with ι given by Eq. (32). Gillet-Chaulet et al. (2006) and Ma et al. (2010) used Eq. (44) with the GOLF rheology (Eq. 35). However it should be noted that Eq. (31), i.e. Eq. (44) with α=1, is derived from the Sachs approximation, which results in the leading E_ca term, and which is cancelled out when considering the normalised stress $\widehat{\mathit{\tau }}$. It is more accurate to consider the expressions for the normalised stress $\widehat{\mathit{\tau }}$ if one is postulating this relationship with a rheology other than Eq. (23), and then use these expression in Eq. (45) rather than Eq. (44). Despite this, multiple studies have used values of ι=10 in Eq. (45) (e.g. Ma et al., 2010; Lilien et al., 2021). Work by Richards et al. (2021) has also used a free parameter in front of the contribution from $\dot{\mathit{\epsilon }}$. Consequently, the most general form for v which captures all previous work is:

$$\begin{array}{}\text{(46)}& \begin{array}{rl}{v}_i& ={\mathit{\omega }}_{ij}{c}_j-{\mathit{\alpha }}_{\mathrm{D}}({\dot{\mathit{\epsilon }}}_{ij}{c}_j-{\dot{\mathit{\epsilon }}}_{jk}{c}_k{c}_j{c}_i)\\ & -{\mathit{\alpha }}_{\mathrm{S}}\mathit{\iota }({\widehat{\mathit{\tau }}}_{ij}{c}_j-{\widehat{\mathit{\tau }}}_{jk}{c}_k{c}_j{c}_i),\end{array}\end{array}$$

where ι is given by Eq. (32) and ${\mathit{\alpha }}_{\mathrm{D}},{\mathit{\alpha }}_{\mathrm{S}}\ge \mathrm{0}$. This is the equation we use when modelling fabrics and testing the predictions of different rheologies. The lattice rotation v and hence fabric evolution implicitly depends on the rheology through the appearance of $\widehat{\mathit{\tau }}$ in Eq. (46).

5 Testing the model: Application to ice sheets

In Eqs. (21) and (46) we have an equation for fabric evolution which depends on a few parameters and on the anisotropic rheology used. Our goal in this work is twofold. Firstly to provide a model and set of parameters that can predict the fabrics observed in a range of locations in the natural world, including ice streams. Secondly, to provide insight into which anisotropic rheology is most accurate. To tackle both these problems we apply the fabric evolution model to two locations in the Greenland ice sheet: the GRIP site, at an ice divide, and the EGRIP site, at an ice stream (Fig. 3). These sites are chosen as they have extensive ice core data (Thorsteinsson et al., 1997; Stoll et al., 2021), giving all three eigenvalues of the 2nd-order orientation tensor, A, as a function of depth. They also have similar strain rates and temperatures, meaning the same parameter set will be valid for both locations.

At each location, we find the best fit parameters (by minimising the difference between modelled and measured eigenvalues) to observations for each rheology. Consequently, we can compare the predictions of different rheologies or models against each other. Our hypothesis is that a model and rheology combination which actually captures the underlying physics should be able to accurately predict the observed fabrics at both locations with the same parameter set, despite the very different deformation conditions. A model that does not sufficiently capture the physics may be able to give good predictions at one of the locations just by fitting a large number of parameters to the observations, but is less likely to be able to give accurate predictions at both locations.

The fabric is simulated in a Lagrangian sense, by tracking the evolution of the fabric attached to an “ice parcel” as it moves through the ice sheet. As the ice parcel moves through the ice sheet, we consider how the velocity gradient, and hence ω and $\dot{\mathit{\epsilon }}$ evolve. These represent the inputs to the fabric evolution model, with the output being the fabric and second-order orientation tensor to be compared to ice core observations.

We solve the fabric using a Monte-Carlo method, directly solving equations for $\mathrm{d}\mathit{c}/\mathrm{d}t$ (Eq. 46) and $\mathrm{d}w/\mathrm{d}t$ (Eq. 19) for a large number (2000) of samples, and then calculating the orientation tensors through Eq. (26).

5.1 Ice divide

Fabric models have commonly been compared to observations at ice divides (Castelnau et al., 1996b; Montagnat et al., 2012). To predict fabrics at a specific ice divide and hence compare to observational measurements from drilled ice cores, we need to approximate the deformation $\dot{\mathit{\epsilon }}$ as the input to the fabric evolution model. To do this a number of assumptions must be made. Firstly, it is assumed the core is drilled at a perfect ice divide, such that it only experiences unconfined compression:

$$\begin{array}{}\text{(47)}& \mathrm{\nabla }\mathit{u}=\left[\begin{array}{ccc}-{\dot{\mathit{\epsilon }}}_{zz}/\mathrm{2}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -{\dot{\mathit{\epsilon }}}_{zz}/\mathrm{2}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\dot{\mathit{\epsilon }}}_{zz}\end{array}\right]\end{array}$$

Secondly, The vertical shear rate ${\dot{\mathit{\epsilon }}}_{zz}$ is assumed to follow a Dansgard-Johnson profile (Dansgaard and Johnsen, 1969), being constant for the upper two-thirds of the ice divide, and then decreasing linearly to zero at the ice-bedrock interface. Once this is assumed, the value of ${\dot{\mathit{\epsilon }}}_{zz}$ in the upper two-thirds, ${\dot{\mathit{\epsilon }}}_{zz}^{\prime }$ can be calculated:

$$\begin{array}{}\text{(48)}& \underset{\mathrm{0}}{\overset{H}{\int }}{\dot{\mathit{\epsilon }}}_{zz}dz={\dot{\mathit{\epsilon }}}_{zz}^{\prime }{\displaystyle \frac{\mathrm{2}H}{\mathrm{3}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\dot{\mathit{\epsilon }}}_{zz}^{\prime }{\displaystyle \frac{H}{\mathrm{3}}}=\dot{a}\end{array}$$

where H is the ice thickness and $\dot{a}$ is the accumulation. We make a further assumption that the values for H and $\dot{a}$ are unchanged over time. For GRIP, H=3029 m and $\dot{a}=\mathrm{0.23}$ a⁻¹ (Thorsteinsson et al., 1997). As we use the same parameters throughout the entire depth of the ice divide, we are implicitly assuming that the temperature and deformation rate stay approximately constant. The temperature at GRIP is between −31 and −20 °C above a depth of 2500 m, but increases rapidly over the deepest 500 m to −8 °C at the bedrock (Johnsen, 2003). Therefore, we exclude this region from our analysis.

5.2 Ice stream

To predict fabrics at the EGRIP drill site, located in the centre of an ice stream and where ice fabric measurements are available (Stoll et al., 2021), we must estimate the three-dimensional path and deformation of the each “ice parcel” which ends up along the vertical depth of the drilled ice core (Fig. 3). The ice parcel undergoes a changing history of deformations as it moves through the ice sheet.

The same methodology as Richards et al. (2023) is used to estimate the deformation. To leading order, we can assume that the ice flow around EGRIP obeys the shallow stream approximation (MacAyeal, 1989), such that the ice flows with an approximate plug flow, with negligible vertical shear. This means that the horizontal velocities (and their corresponding horizontal gradients) can be applied into the ice sheet to leading order. These horizontal velocities are found using satellite velocity data (Joughin et al., 2016, 2018).

Because the horizontal velocities are applied into the depth of the ice sheet, all the three-dimensional particle paths sketched in Fig. 3 follow the same surface path. Different initial locations along this surface path create three-dimensional particle paths that intersect the EGRIP drill site at different depths. At the starting point of each path on the surface, the vertical velocity can be found through the kinematic condition:

$$\begin{array}{}\text{(49)}& w_s=-\dot{a}+u_{\mathrm{s}}{\displaystyle \frac{\partial h}{\partial s}},\end{array}$$

with $\dot{a}$ the accumulation rate, u_s the surface velocity in the streamline direction and $\partial h/\partial s$ the change in surface height in the streamline direction, calculated from surface height data (Morlighem et al., 2017). This can then be integrated forward to give a three-dimensional path from the surface to some depth at EGRIP. For each path we iteratively choose the accumulation rate such that the final age and depth of an ice parcel at EGRIP match the measured age-depth relationship (Gerber et al., 2021).

The measured temperature along the EGRIP core is approximately −30 °C to a depth of 1400 m, (Prior-Jones et al., 2021) – the depth to which the fabric data in this paper extends to. This temperature is similar to what is measured at GRIP at these depths (Johnsen, 2003). We assume this temperature-depth relationship is also valid upstream, along the particle paths. This means we can say the temperature is approximately constant throughout.

Figure 3 Schematic showing the location of GRIP (at an ice divide) and EGRIP (in an active ice stream). The cores drilled from GRIP all lie on the same streamline, whereas the cores drilled from EGRIP correspond to different streamlines all starting from different points upstream.

6 Results

We now move on to the results from applying described rheologies to these two locations in the Greenland ice sheet. We first present the results of the parameter fitting for each model/rheology, then show how these best parameters for each model agree with the ice core observations.

6.1 Parameter estimation

To re-summarise, there exists one equation for lattice rotation that contains all previously described models:

$$\begin{array}{}\text{(50)}& \begin{array}{rl}{v}_i& ={\mathit{\omega }}_{ij}{c}_j-{\mathit{\alpha }}_{\mathrm{D}}({\dot{\mathit{\epsilon }}}_{ij}{c}_j-{\dot{\mathit{\epsilon }}}_{jk}{c}_k{c}_j{c}_i)\\ & -{\mathit{\alpha }}_{\mathrm{S}}\mathit{\iota }({\widehat{\mathit{\tau }}}_{ij}{c}_j-{\widehat{\mathit{\tau }}}_{jk}{c}_k{c}_j{c}_i),\end{array}\end{array}$$

The deviatoric stress $\widehat{\mathit{\tau }}$ depends on the rheology used. Consequently, when the lattice rotation depends on the stress the fabric evolution becomes coupled to the rheology used. This means we can compare fabric predictions to ice cores to test the accuracy of different rheologies.

We review 6 different rheologies models: The Taylor and Sachs models derived from grain dynamics, as well as the 4 different macroscopic rheologies: GOLF, NLO, CAFFE and Glen's flow law.

The rheology of a grain, and consequently the large-scale, is described by two parameters E_cc and E_ca controlling the increased deformation experienced by a grain in compression or shear (Eq. 1). The parameters α_D and α_S control the contribution from deformation and stress to lattice rotation respectively in Eq. (50). There are two parameters controlling recrystallization processes: λ and β controlling rotational and migration recrystallization respectively. This describes fully most rheologies. The NLO also depends on the power-law exponent n, and the CAFFE rheology (Eq. 39) has unique maximum and minimum softening parameters E_max and E_min.

We non-dimensionalise the two recrystallization parameters by the macroscopic effective strain rate $\dot{\mathit{\gamma }}=\sqrt{{\dot{\mathit{\epsilon }}}_{ij}{\dot{\mathit{\epsilon }}}_{ji}/\mathrm{2}}$, such that $\mathit{\lambda }=\stackrel{\mathrm{̃}}{\mathit{\lambda }}\dot{\mathit{\gamma }}$, $\mathit{\beta }=\stackrel{\mathrm{̃}}{\mathit{\beta }}\dot{\mathit{\gamma }}$, as in Richards et al. (2022).

While there are a large number of parameters, there are a number of physical constraints available to reduce the parameter space. The Taylor model requires ${\mathit{\alpha }}_{\mathrm{D}}=\mathrm{1},{\mathit{\alpha }}_{\mathrm{S}}=\mathrm{0}$ in Eq. (50) (and hence does not involve a rheology), while the Sachs model requires ${\mathit{\alpha }}_{\mathrm{D}}=\mathrm{0},{\mathit{\alpha }}_{\mathrm{S}}=\mathrm{1}$. The other 4 rheologies place no restriction on the values of α_D and α_S. Furthermore, there is no evidence to suggest that compression is easier or harder in the direction of the c-axis or perpendicular to it (Gillet-Chaulet et al., 2005), which means E_cc=1. For E_ca, Rathmann and Lilien (2021) finds a value of 10⁴ best reproduces the expected polycrystal enhancement of 10. However we use E_ca=10² as this is the maximum value the GOLF rheology is constrained against, and we aim to make comparisons between rheologies. This has only a minimal effect on the results here.

The GRIP ice divide (above a depth of 2500 m) and the EGRIP ice stream are at $T<-\mathrm{20}$ °C and strain rates $<{\mathrm{10}}^{-\mathrm{11}}$ s⁻¹. Migration recrystallization is known to be predominantly a high temperature (Craw et al., 2018; Qi et al., 2019) or potentially also a low temperature high strain rate phenomena (Faria et al., 2014a). As both GRIP and EGRIP are at low temperature and low strain rates (relative to laboratory strain rates), we do not think migration recrystallization will be significant for either location and set it to zero.

As it is the normalised stress which appears in Eq. (46), most models do not depend on the value of n used. The exception is the non-linear orthotropic rheology, for which we use n=3. Furthermore, when rheologies have internal parameters we use those suggested by their authors. The non-linear orthotropic rheology depends on a Sachs-Taylor weight, which we use 98.75:1.25 following Rathmann and Lilien (2021). The CAFFE model depends on parameters E_max and E_min, for which we use the values of 10 and 0.1 respectively, as is suggested in Placidi et al. (2010). While the CAFFE model does not depend on E_cc and E_ca, we use the same value of ι in Eq. (46) for ease of comparison.

For the Sachs and Taylor models, this leaves only $\stackrel{\mathrm{̃}}{\mathit{\lambda }}$ to be constrained. For the GOLF, non-linear orthotropic and CAFFE rheologies 3 parameters must be determined: α_D, α_S and $\stackrel{\mathrm{̃}}{\mathit{\lambda }}$. Glen's rheology has only 2 free parameters: α_D and $\stackrel{\mathrm{̃}}{\mathit{\lambda }}$ because ${\dot{\mathit{\epsilon }}}_{ij}={\widehat{\mathit{\tau }}}_{ij}$.

For all models/rheologies the parameter space is explored and a set of best parameters are chosen (Table 2), aiming to give the best fit for both locations. There is not one unique set of parameters which gives the unambiguous best fit for each model, but the parameters do show the best possible fit for each model/rheology.

Table 2 Comparison best fit parameters of different model and rheologies. Parameters in italics were chosen by fitting, blank entries mean the parameter is not included in the model.

The fabric predictions from the parameters shown in Table 2, for each model/rheology, are compared to observations from the GRIP ice divide and the EGRIP ice stream.

6.2 GRIP ice divide

At the ice divide we assume the fabric undergoes unconfined compression as it moves vertically down. The measured largest eigenvalue from different depths at GRIP is shown in Fig. 4 as crosses (Thorsteinsson et al., 1997). The same data points are shown on both subplots. As the fabric at an ice divide has rotational symmetry around the vertical axis, the two smaller eigenvalues are equal. Consequently, we only plot the largest eigenvalue.

Figure 4 also shows the results from the suite of models/rheologies tested, plotted as solid lines. The figure has been split into two subplots to aid readability. The tensor rheologies are on the right-hand subplot. Because of the way we chose parameters, first finding a set that agreed well with this case and then improving the fit for the ice stream, all the models in the right subplot agree well with observations here. There is very little difference between the best predictions of the Sachs and GOLF. The NLO rheology has marginally worse agreement with ice cores above a depth of 1500 m.

In the left subplot of Fig. 4, the Glen rheology with α_D=2.6, α_S=0 and $\stackrel{\mathrm{̃}}{\mathit{\lambda }}=\mathrm{0.32}$ also gives good agreement with the observed eigenvalues. For this rheology $\widehat{\mathit{\tau }}=\dot{\mathit{\epsilon }}$, therefore this parameter set could equivalently be expressed as α_D=0, α_S=1.06, $\stackrel{\mathrm{̃}}{\mathit{\lambda }}=\mathrm{0.32}$, for the same value of ι=2.46 in Eq. (46), as in the other models. This is a similar parameter set to that found for the Sachs model, highlighting that for this simpler case of the ice divide, it is not necessary to have an anisotropic rheology contributing to the fabric evolution to get reasonably accurate predictions.

The CAFFE model exhibits a sharp increase of the largest eigenvalue to ≈1 at a depth of around 1500 m. This is not seen in observations for GRIP but has been observed to an extent at other ice divides (Montagnat et al., 2012). We find that this sharp increase is sensitive to the value of E_min, though detailed exploration of this is beyond the scope of this work.

Figure 4 Evolution of largest (vertical) eigenvalue of the second-order orientation tensor A at the GRIP drill site (an ice divide) with depth. Only the largest eigenvalue is shown as the fabric has rotational symmetry about the vertical, so the other two eigenvalues are equal. Measurements from ice cores are plotted as crosses, showing an increase with depth. Model results are plotted as lines, and have been split across two subplots to aid readability. The right hand plot shows the tensorial rheologies.

6.3 EGRIP ice stream

The flow at an ice stream is much more complex than an ice divide. Consequently, at each depth at the EGRIP ice core the ice there has experienced changing deformation history. In Fig. 5 we show the observed eigenvalues at EGRIP with depth as blue dots. Each eigenvalue is shown in a separate column. The data is also reproduced across each row to aid readability. Close to the surface at a depth of around 50 m, the observed eigenvalues at EGRIP have one larger eigenvalue (the vertical) and two smaller, indicating a single-maximum fabric produced as snow compresses into ice. We do not model this effect. Below this, the fabric is produced by ice stream deformation, and is characterised by one eigenvalue ≈0, and two larger eigenvalues, corresponding to a girdle fabric. This fabric is fairly invariant with depth below 550 m, with some noise.

We also show the results from the suite of models/rheologies tested in Fig. 5. The first row shows results from the Taylor, Glen and CAFFE models, and the second row shows the results from the tensorial rheologies. We do not include the effect of firn in our analysis, which significantly affects the observed fabric above a depth of around 550 m. Hence, we only compare the models to observations below a depth of 550 m. Although different initial conditions for the fabric could perhaps fit the observations better, we find that the results are insensitive to initial condition below a depth of 400 m (Fig. S13).

The Taylor model does not give good agreement with the EGRIP observations, showing one much larger eigenvalue of ≈0.8 over the depth we are interested in (below 550 m). The best parameters with Glen's rheology and the GOLF rheology also do not agree, not predicting that the smallest eigenvalue is much larger than what is observed. From this figure, the CAFFE rheology, Sachs rheology, and the NLO rheology give reasonably good agreement with observations, predicting a smaller eigenvalue close to zero and two larger eigenvalues.

To provide a more detailed comparison, Fig. 6 shows a comparison of modelled to observed eigenvalues in a depth-averaged sense. Observed eigenvalues from EGRIP are averaged over a depth of 550 to 1425 m, and plotted as dashed lines with the standard deviation shown by the shaded region. In this figure, the colour coding simply shows which eigenvalue (largest, middle or smallest) is being referred to. Then for each model/rheology (x-axis), we average the modelled eigenvalues over the same depth and plot them as as crosses, with the same colour coding as observations. Distance of the crosses from the respective dashed lines then allows us to assess agreement of models with observations in this deeper area of the ice core.

Figure 6 shows that the CAFFE, Sachs, and NLO rheologies all have their predictions for the largest two eigenvalues within 1 standard deviation of the mean observed value at EGRIP. They also all have the smallest eigenvalue close to zero, with the Sachs model being closest to observations. This figure also shows that while Glen's and the GOLF rheologies with their respective best parameters are less accurate, they are both significantly more accurate than the Taylor model.

Previously, we have limited our analysis to solely the eigenvalues. We now extend this to examine the pole figures – showing the complete orientation distribution function (ODF) – directly in Fig. 7. The ODF represents the probability of finding a grain at each orientation. Brighter regions indicate a high probability of a grain being orientated towards that orientation. The space of possible orientations exist on the surface of a sphere. There is also no possibility of distinguishing a grain's top from bottom, so the ODF is symmetrical between its north and south hemispheres. Consequently, we plot the ODF over a single hemisphere, looking vertically downwards on the fabric. Figure 7 shows two pole figures observed from EGRIP on the left. While we have eigenvalue data over a large depth range, we only have pole figures for these two depths. There is variance between these two depths. The corresponding eigenvalues to these EGRIP pole figures are shown in Fig. 5 as black crosses, showing that these two pole figures represent end members of the variance seen at EGRIP. Despite this, they both show a girdle fabric pattern, indicating most grains have c-axes pointed along the bright band.

We also show pole figures for each model in Fig. 7. The Taylor model had one large eigenvalue and the other two ≈0. This corresponds to what is seen in the pole figure: a cluster in only 1 direction, which is repeated due to the symmetry of the ODF. Glen's rheology and GOLF rheologies both predicted the smallest eigenvalue was larger compared to observations, and this corresponds to a generally too diffuse pole figure, which was also seen in Richards et al. (2023), which was the previous state-of-the-art for modelling ice stream fabrics.

The remaining models all provide improvements on this previous state-of-the-art. The pole figures for the Sachs model and NLO rheology lie within the variance of what is seen in the EGRIP pole figures, predicting girdle fabrics with a slight bias towards grains being orientated towards the edges of the pole figure. The CAFFE model also predicts a girdle fabric, though the detail features of the pattern in Figure 7 show it is slightly different from what is seen at EGRIP, with peaks closer to the centre of the figure. These detail features can only be examined through looking directly at the pole figures.

These results are summarised in Fig. 8. Here we show the summed difference between the modelled and observed eigenvalues for each model/rheology. This error is then normalised against the results for Glen's rheology (so the results for Glen are by definition 100 %). This again highlights that the Sachs model has the closest agreement for both the ice divide and ice stream. Figure 8 also shows how the anisotropic rheologies provide the significant improvements in fabric predictions for the ice stream case (EGRIP) only.

Figure 5 This figure shows the measured eigenvalues at EGRIP, plotted as dots, alongside the results for each model/rheology (lines). Each column represents a single eigenvalue (smallest, middle, largest) and the figure has also been split over two rows to aid readability, with the tensorial rheologies in the second panel. Highlighted with black crosses are the eigenvalues of the pole figures shown in Fig. 7.

Figure 6 This figure shows, for each model/rheology, the mean eigenvalues over a depth of 550 to 1425 m. Colours indicate the different eigenvalues. The horizontal line shows the mean of the observed eigenvalues over this depth, the shaded region corresponding to ±1 standard deviation. Along the x-axis, the mean eigenvalues for each model are plotted as crosses.

Figure 7 Comparison of pole figures, for all models and measured from an ice core at two depths, for which the corresponding eigenvalues are highlighted in Fig. 5. The pole figure shows a hemisphere of the orientation distribution function, which is sufficient to show all the information. Pole figures are shown looking vertically downwards, with North at the top. For example, the observed fabric at EGRIP at 1048 m shows that almost all grains lie close to a plane normal to the North-East direction.

Figure 8 Normalised (relative to best parameters with Glen's rheology) error for both GRIP and EGRIP. This is calculated by considering the sum of the difference between observed and measured eigenvalues at each depth.

6.4 Sensitivity to changes in the ice stream over time

In our previous analysis, as explained in Sect. 5.2, we use the present-day satellite velocity data to trace ice parcels upstream. However tracing a particle path upstream also means tracing a particle path back in time. Consequently, this method implicitly assumes the ice velocity field has not changed over the past 16 000 years (the time it takes for ice to travel from the furthest point upstream to EGRIP).

Recent work has raised the possibility that the North-East Greenland Ice Stream (NEGIS) may be relatively young, existing only for 2000 years (Jansen et al., 2024), before which no ice stream existed here. To test the robustness of our results to a “young NEGIS”, we assume that before a certain time, there was no ice stream, i.e. the ice was frozen to the bed. Consequently, the deformation would be dominated by vertical shear and a single-maximum would be produced, as is observed at locations like NEEM and NGRIP (Montagnat et al., 2014a). Ice cores from these locations reveal a vertical single maximum with a largest eigenvalue of around 0.75 over the depth we examine. Therefore we assume this as the pre-existing fabric over the whole depth. The ice stream deformation is then switched on at a certain time ago, and the fabric evolution is calculated from then onwards, with the vertical single-maximum as the initial fabric.

We explore this in Fig. 9. This again shows how the average eigenvalues evolve as in Fig. 6, but this time depending on the age of the NEGIS. The results are shown only for the Sachs model with the parameters shown in Table 2, though the trend does not change between models. Figure 9 shows that if the ice stream was extremely young (e.g. 100 years) the fabric would still be a vertical single maximum. If the NEGIS is was only 2000 years old, our results change slightly: a girdle fabric is still predicted, but the eigenvalues are different, consequently the best parameters would be slightly different. However, even for a young NEGIS of only 3000 years or older, the results are unchanged. Consequently, we can say our results are robust to an ice stream age of 3000 years or older. If the ice stream is found to be only 2000 years old, our best fit parameters change, but the general conclusions and ranking of the different models do not. Reproductions of Table 2, and Figs. 4, 6 and 7 for a 2000 year old ice stream can be found in the Supplement.

Figure 9 This figure shows how the results for the Sachs model change at EGRIP, for an ice stream that only started x000 years ago (x-axis). Plotted is the modelled mean eigenvalue as the solid line. For comparison, the mean and standard deviation of the observed eigenvalues at EGRIP is also shown as in Fig. 6, as the dotted and shaded region respectively. Colours indicate the eigenvalue (i.e. smallest, middle, largest).

7 Discussion

In the introduction, we set out two main problems facing accurate modelling of viscous anisotropy in ice sheet models. Firstly, a model for fabric evolution which can give accurate predictions of fabrics observed in the natural world, specifically in ice streams, is missing. Secondly, there has been no way to test whether an anisotropic rheology is accurately representing the anisotropic effect, leading to a large number of competing rheologies.

In the first part of this paper (Sects. 2, 3 and 4) we examined a range of existing models of anisotropy in ice, starting from the dynamics of an ice grain. We showed that previous models for fabric evolution (e.g. Azuma and Goto-Azuma, 1996; Gillet-Chaulet et al., 2006; Richards et al., 2021; Lilien et al., 2021) can be expressed as a common equation (Eq. 46) depending on a few parameters and also the rheology used. Therefore, the choice of anisotropic rheology becomes part of the fabric evolution model.

In Sects. 5 and 6 we applied these models to two locations in the Greenland ice sheet: an ice divide and ice stream, resulting in significantly more accurate results than previous work (Richards et al., 2023). We found the best fit parameters for each rheology. These best fit parameters (Table 2) suggest that lattice rotation due to stress, consistent with the assumption that all grains are experiencing the same stress (Eq. 31), is the most physically accurate. Of the rheologies examined, we found that the rheologies with best agreement to observations in both the ice stream (Fig. 5) and ice divide (Fig. 4) were the Sachs rheology (Eq. 23) and the non-linear orthotropic rheology (Eq. 36). Both these rheologies describe anisotropy through a fourth-rank tensor, rather than attempting to represent it through some scalar enhancement factor. They both also assume the ice grains are experiencing approximately the same stress.

7.1 How is the ice grain deforming?

The results we have shown support the case for ice grains deforming closer to the Sachs hypothesis – that all grains experience the same stress – than the Taylor hypothesis – that all grains experience the same deformation. We conclude this because of the two models derived directly from grain dynamics, the Sachs model predicts eigenvalues much closer to observations, for both the divide (Fig. 4) and the ice stream (Fig. 5), than the Taylor model. The Taylor model is also not able to predict the girdle fabric observed in the ice stream (Fig. 7). Furthermore, we observed that the Sachs model had the closest agreement to observations of all the examined models in this work. This is despite the fact that the other rheologies have more free parameters, which would be expected to result in better predictions. Indeed, for the anisotropic rheologies where α_D and α_S (controlling lattice rotation due to deformation and stress respectively) can take any value, the best predictions come from choosing values close those required by the Sachs model.

One caveat to this is the variable amount of diffusion/Browninan motion (controlled by λ) required to give accurate predictions. If there is no diffusion acting on the fabric (λ=0), the Sachs model vastly over predicts the fabric strength at both locations. This diffusional process is meant to mimic rotational recrystallization. Observations of grain boundaries suggest some rotational recrystallization is indeed occurring in the EGRIP (Stoll et al., 2021). However it is possible that the some amount of the modelled diffusion of the fabric is acting to ameliorate an over-prediction of lattice rotation from the Sachs hypothesis alone. Castelnau et al. (1996b), Thorsteinsson (2002) showed that including nearest neighbour interactions on grains can have a diffusional-like effect on the fabric. Consequently, we are cautious about using this model to infer rates of rotational recrystallization.

7.2 What is needed for accurate fabric modelling?

Previous work including only lattice rotation due to deformation has indeed succeeded at ice divides (Lilien et al., 2021; Richards et al., 2023). Yet it is worth noting that at divides the deformation is steady over time, and the observed fabric is simple – a single maximum. This makes it fairly easy for a model with a number of free parameters to fit to observations. In contrast, ice streams have a changing deformation history and the observed fabric at EGRIP is more complex. Models which take only lattice rotation due to deformation into account have not been able to accurately model this fabric (Gerber et al., 2023; Richards et al., 2023). Conversely, we observe that lattice rotation due to stress rather than deformation is key, combined with an anisotropic rheology like the Sachs or the NLO rheologies (Fig. 7). The parameter fitting in Sect. 6.1 finds that for all models with free parameters controlling both lattice rotation due to deformation (α_D) and due to stress (α_S), the best agreement comes from lattice rotation due to stress only (Table 2).

Furthermore, the fabric evolution model is only able to accurately predict the fabric observed at EGRIP when it is linked to either the Sachs or NLO rheologies. What these rheologies have in common is that they (a) have a fourth-rank tensor description of the anisotropy and (b) assume the grains deform with approximately the same stress. The other models are not able to predict the observed fabrics at EGRIP for any set of parameters chosen. The Sachs and NLO rheologies more accurately predict the form of the stress tensor, and are consequently able to accurately predict the eigenvalue evolution and correct pattern for both an ice stream and an ice divide, two very different deformation conditions. Consequently, a fabric model with either the Sachs or NLO rheology, along with the parameters shown in Table 2, can be applied more generally and can be used in ice sheet models, at least for low temperature conditions ($T<-\mathrm{20}$ °C) which apply to GRIP and EGRIP.

As explored in Sect. 7.4, we have tested the models at only two locations in the ice sheet. Further in-situ observations across a range of deformation and temperature conditions are crucial for accurate modelling of fabrics and anisotropy. Recent advances with radar (e.g. Young et al., 2021) and seismics (e.g. Lutz et al., 2020) are welcome, but ice cores, which provide a full fabric dataset, are still needed. These would be particularly beneficial in areas with deformation similar to EGRIP, but at higher temperatures.

7.3 Which rheologies are most appropriate?

By incorporating the rheology into the fabric evolution model, we are able to compare the accuracy of different rheologies and provide a test of them for the first time. All the anisotropic rheologies provide an improvement over Glen's rheology. For each rheology, we have chosen the set of parameters which minimise the difference between the modelled and observed eigenvalues. Despite this, only the Sachs and non-linear orthotropic rheology are able to predict the observed girdle fabric at EGRIP, and have the minimum difference between the modelled and the observed eigenvalues.

Despite providing accurate fabric predictions, the Sachs rheology under-predicts the maximum possible effect of fabric on viscosity. Experimental tests (Pimienta and Duval, 1987), have shown a vertical single maximum fabric has approximately 10 times higher viscosity in the vertical direction compared to an isotropic fabric (i.e. a fabric with the grains randomly orientated). However, the Sachs rheology (with linear grain behaviour as used here) predicts a maximum ratio in viscosity between these two fabrics of only 2.5 (Gagliardini et al., 2009). This effect does not impact fabric predictions due to our use of a normalised stress throughout (Eq. 29).

Fortunately, this disadvantage does not affect the NLO rheology. To give the correct maximum enhancement with the NLO rheology, an E_ca value of 10³ should be used, along with the Sachs-Taylor weight used here. As is seen in Fig. S9 of the Supplement, there is only a small effect of increasing E_ca above 100 on fabric predictions, as its effect on the eigenvalues asymptotically approaches some value as E_ca→∞. Therefore, we recommend using the NLO rheology in ice-sheet models.

It is surprising that the GOLF rheology performs poorly. The GOLF reproduces the response (through tabulated values) of the Visco-Plastic Self-Consistent (VPSC) model, a microscopic model which includes grain interactions, and should theoretically be more accurate than the Sachs or Taylor bounds, or a linear combination of them. A possible explanation for this is that the underlying VPSC model is outdated: Montagnat et al. (2012) notes that this VPSC model cannot accurately predict fabrics in simple shear, which occurs in the ice stream. This has been improved upon (Signorelli and Tommasi, 2015), so it would be interesting to see how the GOLF performs with an updated underlying model.

Finally, the CAFFE model (Placidi et al., 2010) gives some interesting predictions. It gives good predictions for the eigenvalues at EGRIP. However, for the ice divide case it predicts a sharp increase of the largest eigenvalue at 1500 m which is not seen in other models or the ice core observations. This sharp increase is sensitive to the parameter E_min (see Fig. S12 of the Supplement) whereas Placidi et al. (2010) claim this parameter can be freely chosen between 0 and 0.1. Nevertheless, it provides a significant increase in accuracy for fabric predictions over Glen's rheology, and is easier to incorporate into ice-sheet models than a tensor flow relation.

7.4 Limitations of this analysis

One limitation of this work is that we have only explored one ice divide and ice stream each. There exist other observations from different divides (see Faria et al. (2014a) for a review). Fortunately, the observed fabrics at different ice divides are generally similar, so our results applied to GRIP should also apply elsewhere.

Unfortunately, there is only one ice core drilled from an active ice stream: EGRIP. Recent advances with radar (e.g. Young et al., 2021) and seismics (e.g. Lutz et al., 2020) provide some information on fabrics at other ice streams, however they do not provide the full eigenvalue and pattern information with depth that is available from an ice core. Now that we have a fabric model we believe is accurate, it can be applied to locations where radar and seismic measurements are available to see if it agrees with these observations. Ideally, there would be another full ice core available from an active ice stream.

We have also not considered the effect of temperature. In fact, while the parameters we have found give good predictions for both EGRIP and GRIP, and can be applied to low temperature regions ($T<-\mathrm{20}$ °C), further work is required to constrain these parameters at higher temperatures. At the lattice scale, far too small to resolve in this model, all the processes are dependent on temperature. At the larger scale, this means the parameters controlling recrystallization processes and lattice rotation will be temperature dependent. Migration recrystallization is known to be highly temperature dependent (Richards et al., 2021). However, the comparison to GRIP and EGRIP only covers a temperature range of −40 to −20 °C. Consequently, the parameters shown do not cover the highest temperature seen in ice sheets. Although we have neglected migration recrystallization in this study, it is expected to be significant at higher temperatures (Faria et al., 2014b; Richards et al., 2021). These high temperature regions include ice shelves and the ice-bedrock interface, both areas of critical importance to the large-scale flow.

As discussed in Sect. 6.4, as part of this analysis we have assumed the deformation around EGRIP is unchanged going back around 16 000 years. However, there is open discussion in the literature over whether the ice stream at EGRIP is a recent formation. Some recent work (Jansen et al., 2024) has suggested it may be as young as 2000 years, although a region of convergent ice flow could have been present before this which would produce a similar fabric. We have shown in Fig. 9 that the results used here are unchanged for an ice stream age as young as 3000 years. For an age of 2000 years, the best fit parameters are changed slightly (Table S1), but the overall conclusions and ranking of best rheologies are unchanged.

8 Conclusions

In the introduction, we mentioned that there are currently two main problems in modelling anisotropy: (i) what is the correct physics for modelling fabric evolution and (ii) how can we know if competing anisotropic rheologies in the literature are accurately representing the anisotropic effect. In this work we have made a step towards addressing both these questions. This is done by first showing that many previous fabric evolution models for ice can be written as one common equation, subject to different rheologies and parameters. These different models and rheologies are then tested against observations and their results compared.

We find that only with certain physical assumptions can accurate predictions of ice fabrics be obtained across a range of natural conditions, including predicting the girdle fabric observed at EGRIP. The physics required is lattice rotation due to stress, where the stress is calculated using an anisotropic rheology. This is either the Sachs rheology (which assumes all grains experience the same stress) or the non-linear orthotropic rheology (Rathmann and Lilien, 2022) (which is close to this approximation). Alongside this, it is necessary to include a diffusive process acting on the fabric modelling rotational recrystallization.

We also provide a comparison of competing anisotropic rheologies, by testing which rheologies give the best fabric predictions. As mentioned above, the Sachs rheology and the non-linear orthotropic rheology give the best predictions, suggesting they accurately represent the anisotropic effect. However, as the Sachs rheology under predicts the maximum softening observed for strong fabrics, we suggest the rheology of Rathmann and Lilien (2022) is most physically accurate overall.

While further work is needed to find the correct parameters for fabric evolution at higher temperatures (where migration recrystallization will become important) this works provides significant insight into the physical assumptions to use when modelling anisotropy in ice-sheet models.

Appendix A: Expressing grain rotation in terms of stress

Here we derive Eq. (6) from combining Eqs. (2) and (5). Restating Eqs. (2) and (5):

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}{\dot{\mathit{\epsilon }}}_{ij}^{\prime }& ={\mathit{\eta }}^{-\mathrm{1}}\left({\mathit{\tau }}_{ij}^{\prime }+{\displaystyle \frac{\mathrm{3}(E_{cc}-\mathrm{1})-\mathrm{4}(E_{ca}-\mathrm{1})}{\mathrm{2}}}{\mathit{\tau }}_{kl}^{\prime }{c}_k{c}_l{c}_i{c}_j\right.\\ & \left.+(E_{ca}-\mathrm{1})({\mathit{\tau }}_{ik}^{\prime }{c}_k{c}_j+c_i{c}_k{\mathit{\tau }}_{kj}^{\prime })-{\displaystyle \frac{E_{cc}-\mathrm{1}}{\mathrm{2}}}{\mathit{\tau }}_{kl}^{\prime }{c}_k{c}_l{\mathit{\delta }}_{ij}\right)\end{array}\text{(A2)}& {\displaystyle \frac{\mathrm{d}{c}_i}{\mathrm{d}t}}={\mathit{\omega }}_{ij}^{\prime }{c}_j-({\dot{\mathit{\epsilon }}}_{ij}^{\prime }{c}_j-{\dot{\mathit{\epsilon }}}_{kj}^{\prime }{c}_k{c}_j{c}_i),\end{array}$$

Taking advantage of the fact that c_ic_i=1, and defining ${\mathit{\tau }}_c^{\prime }={\mathit{\tau }}_{ij}^{\prime }{c}_j{c}_i={\mathit{\tau }}_{ji}^{\prime }{c}_j{c}_i={\mathit{\tau }}_{kl}^{\prime }{c}_k{c}_l$ etc.

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{\dot{\mathit{\epsilon }}}_{ij}^{\prime }{c}_j& ={\mathit{\eta }}^{-\mathrm{1}}\left({\mathit{\tau }}_{ij}^{\prime }{c}_j+{\displaystyle \frac{\mathrm{3}(E_{cc}-\mathrm{1})-\mathrm{4}(E_{ca}-\mathrm{1})}{\mathrm{2}}}{\mathit{\tau }}_c^{\prime }{c}_i\right.\\ & \left.+(E_{ca}-\mathrm{1})({\mathit{\tau }}_{ik}^{\prime }{c}_k+c_i{\mathit{\tau }}_c^{\prime })-{\displaystyle \frac{E_{cc}-\mathrm{1}}{\mathrm{2}}}{\mathit{\tau }}_c^{\prime }{c}_i\right)\end{array}\end{array}$$

and

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{\dot{\mathit{\epsilon }}}_{ij}^{\prime }{c}_i{c}_j& ={\mathit{\eta }}^{-\mathrm{1}}\left({\mathit{\tau }}_c^{\prime }+{\displaystyle \frac{\mathrm{3}(E_{cc}-\mathrm{1})-\mathrm{4}(E_{ca}-\mathrm{1})}{\mathrm{2}}}{\mathit{\tau }}_c^{\prime }\right.\\ & \left.+(E_{ca}-\mathrm{1})({\mathit{\tau }}_c^{\prime }+{\mathit{\tau }}_c^{\prime })-{\displaystyle \frac{E_{cc}-\mathrm{1}}{\mathrm{2}}}{\mathit{\tau }}_c^{\prime }\right)\end{array}\end{array}$$

hence

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}{\dot{\mathit{\epsilon }}}_{kj}^{\prime }{c}_k{c}_j{c}_i& ={\mathit{\eta }}^{-\mathrm{1}}\left({\mathit{\tau }}_c^{\prime }{c}_i+{\displaystyle \frac{\mathrm{3}(E_{cc}-\mathrm{1})-\mathrm{4}(E_{ca}-\mathrm{1})}{\mathrm{2}}}{\mathit{\tau }}_c^{\prime }{c}_i\right.\\ & \left.+(E_{ca}-\mathrm{1})({\mathit{\tau }}_c^{\prime }{c}_i+{\mathit{\tau }}_c^{\prime }{c}_i)-{\displaystyle \frac{E_{cc}-\mathrm{1}}{\mathrm{2}}}{\mathit{\tau }}_c^{\prime }{c}_i\right)\end{array}\end{array}$$

Then, subtracting Eq. (A5) from (A3), and renaming indices in Eq. (A3) so that ${\mathit{\tau }}_{ik}^{\prime }{c}_k\to {\mathit{\tau }}_{ij}^{\prime }{c}_j$:

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& {\dot{\mathit{\epsilon }}}_{ij}^{\prime }{c}_j-{\dot{\mathit{\epsilon }}}_{kj}^{\prime }{c}_k{c}_j{c}_i=\\ & {\mathit{\eta }}^{-\mathrm{1}}\left(({\mathit{\tau }}_{ij}^{\prime }{c}_j-{\mathit{\tau }}_c^{\prime }{c}_i)+(E_{ca}-\mathrm{1})({\mathit{\tau }}_{ij}^{\prime }{c}_j-{\mathit{\tau }}_c^{\prime }{c}_i)\right)\end{array}\end{array}$$

finally giving:

$$\begin{array}{}\text{(A7)}& {\dot{\mathit{\epsilon }}}_{ij}^{\prime }{c}_j-{\dot{\mathit{\epsilon }}}_{kj}^{\prime }{c}_k{c}_j{c}_i={\mathit{\eta }}^{-\mathrm{1}}{E}_{ca}({\mathit{\tau }}_{ij}^{\prime }{c}_j-{\mathit{\tau }}_{kj}^{\prime }{c}_k{c}_j{c}_i)\end{array}$$

which is Eq. (6).