# Rheological control of crystal fabrics on Antarctic ice shelves

Nicholas M Rathmann, David A Lilien, Daniel H Richards, Felicity S

Mccormack, Maurine Montagnat

► To cite this version:

Nicholas M Rathmann, David A Lilien, Daniel H Richards, Felicity S McCormack, Maurine Montagnat. Rheological control of crystal fabrics on Antarctic ice shelves. Journal of Glaciology, 2025, 71, pp.e110. 10.1017/jog.2025.10070 . hal-05377209

HAL Id: hal-05377209

https://hal.science/hal-05377209v1

Submitted on 21 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Article

**Cite this article:** Rathmann NM, Lilien DA, Richards DH, McCormack FS, Montagnat M (2025) Rheological control of crystal fabrics on Antarctic ice shelves. *Journal of Glaciology* **71**, e110, 1–21. https://doi.org/10.1017/jog.2025.10070

Received: 27 March 2025

Revised: 28 May 2025

Accepted: 4 July 2025

**Keywords:**

anisotropic ice; ice dynamics; ice physics; ice
shelves; recrystallization

**Corresponding author:**

Nicholas M. Rathmann;

Email: rathmann@nbi.ku.dk

Rheological control of crystal fabrics on
Antarctic ice shelves

Nicholas M. Rathmann<fnref n="1" /> ID, David A. Lilien² ID, Daniel H. Richards³,<fnref n="4" />,

Felicity S. McCormack<fnref n="5" /> ID and Maurine Montagnat⁶ ID

Abstract

Ice crystal fabrics can exert significant rheological control on ice sheets and ice shelves, poten-
tially softening or hardening anisotropic ice by several orders of magnitude compared to isotropic
ice. We introduce an anisotropic extension of the Shallow Shelf Approximation (SSA), allowing
for fabric-induced viscous anisotropy to affect the flow of ice shelves in coupled, transient simu-
lations. We show that the viscous anisotropy of synthetic ice shelves can be parameterized using
an isotropic flow enhancement factor, suggesting that existing SSA flow models could, with little
effort, approximate the effect of fabric on flow. Next, we propose a new way to directly solve for SSA
fabric fields using satellite-derived velocities, assuming velocities are approximately steady and that
fabric evolution is dominated by lattice rotation with or without discontinuous dynamic recrys-
tallization. We apply our method to the Ross and Pine Island ice shelves, Antarctica, suggesting
that these regions might experience significant fabric-induced hardening and softening depending
on the relative strength of lattice rotation and recrystallization. Our results emphasize the ice-
dynamical relevance of needing to better constrain the strength of fabric processes. This calls for
more widespread fabric and temperature measurements from the field, since measurements are
currently too sparse for model validation.

1. Introduction

Nearly three-quarters of Antarctica's coastline is surrounded by floating ice shelves. The buttressing provided by these ice shelves can influence the stability and mass flux of marine-terminating glaciers across their grounding lines, but quantifying buttressing requires an accurate representation of ice shelves in large-scale flow models (Rignot and others, 2004; Rathmann and others, 2017; Reese and others, 2018; Joughin and others, 2021; Otosaka and others, 2023; Sun and Gudmundsson, 2023) with some caveats (Hulbe and Fahnestock, 2004; Hill and others, 2021). Unlike modeling the flow of grounded ice which depends on information about ice geometry, boundary conditions, thermal state, and material (rheological) properties, floating ice shelves represent a somewhat simpler problem: their boundary conditions are better known (no drag on the bottom ocean boundary), they deform by plug flow (negligible vertical shear), and their geometry can be approximated from surface elevation measurements if assumed to fulfill the flotation criterion.

Despite their simpler setting, many factors complicate the modeling of ice shelves, such as (hydro)fracturing causing locally softer ice (Borstad and others, 2017; Sun and others, 2017), heterogeneous sub-shelf melting (Adusumilli and others, 2018; Wang and others, 2025), and the accretion of marine ice with different thermal and rheological properties (Treverrow and others, 2010; Craw and others, 2023; Miles and others, 2025). Another rheological phenomenon that is often neglected (or at least not explicitly modeled) is the heterogeneous viscous anisotropy caused by an evolving crystal-orientation fabric (henceforth fabric) (Diez and others, 2016; LeDoux and others, 2017; Fan and others, 2021), with potential analogies to continental ice streams insofar as they are facilitated by negligible basal drag due to subglacial melting (Ma and others, 2010; Gerber and others, 2023).

Of the rheological properties that evolve spatiotemporally (Figure 1), ice temperature exerts
in general the most important control. For the range of temperatures observed in ice sheets,
-55°C to 0°C, ice fluidity (ease of deformation) can vary by up to a factor of 1000 (Cuffey
and Paterson, 2010) and is sensitive to liquid water content (De La Chapelle and others, 1999;
Fowler and Iverson, 2023). If deformation is dominated by dislocation creep (typically the case
at high stress), the crystal fabric is arguably the second most important rheological property

© The Author(s), 2025. Published by
Cambridge University Press on behalf of
International Glaciological Society. This is an
Open Access article, distributed under the
terms of the Creative Commons Attribution
licence (http://creativecommons.org/licenses/
by/4.0), which permits unrestricted re-use,
distribution and reproduction, provided the
original article is properly cited.

<footnote n="1">¹Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark; ²Department of Earth and Atmospheric Sciences, Indiana University, Bloomington, Indiana, USA; ³British Antarctic Survey, Cambridge, England, UK;</footnote>

<footnote n="4">⁴The Australian Centre for Excellence in Antarctic Science, University of Tasmania, Hobart, Lutruwita, Australia;</footnote>

<footnote n="5">⁵Securing Antarctica's Environmental Future, School of Earth, Atmosphere and Environment, Monash University, Clayton, Victoria, Australia and ⁶University Grenoble Alpes, CNRS, Grenoble INP, IGE, Grenoble, France</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Qualitative overview of strain-rate enhancing mechanisms in ice, where fadings indicate that some end-member uncertainty exists. In the case of temperature and mean grain size, the offset (reference value at $E = 1$) does not have any physical significance and is therefore centrally aligned. The grain size range shown broadly covers that typically found in ice sheets, and the corresponding enhancements are calculated using the rheology by Goldsby and Kohlstedt (2001) with a grain size exponent of 1.4.</figcaption></figure>

(Duval and others, 1983; Castelnau and others, 1998; Fan and Prior, 2023), unless the ice mass is very fractured. Strongly developed anisotropic ice can be up to 10 times easier to shear and 100 times harder to compress compared to isotropic ice, depending on how favorable grain orientations are to basal glide (Shoji and Langway, 1985; Pimienta and others, 1987; Jacka and Budd, 1989; Azuma, 1994, 1995; Jacka and Jun, 2000) (Figure 1). At sufficiently low stresses, grain boundary sliding is often considered the dominant deformation mechanism, where ice fluidity increases with decreasing grain size (Goldsby and Kohlstedt, 1997, 2001). Although this mechanism can also cause fluidity to vary by orders of magnitude, strain weakening at high homologous temperatures is predominantly controlled by fabric development (Fan and others, 2021). Finally, the impurity content of ice can also exert control; for example, dusty ice from glacial periods is typically twice as soft as interglacial ice (Dahl-Jensen and Gundestrup, 1987).

Modeling the coupled evolution between fabric and flow has received a lot of attention in the literature, although limited to flow-line models primarily due to computational expense (Mangeney and others, 1996; Castelnau and others, 1998; Gillet-Chaulet and others, 2006; Durand and others, 2007; Pettit and others, 2007; Martín and others, 2009; Ma and others, 2010; Martín and Gudmundsson, 2012; Lilien and others, 2021, 2023; Rathmann and Lilien, 2021). The full-Stokes problem (unapproximated momentum balance) for isotropic ice is likewise computationally expensive, making it unfeasible for transient simulations over long time scales. Vertically integrated momentum balance approximations are therefore frequently used (Seroussi and others, 2023), relying on scaling arguments and the shallow aspect ratio of polar ice masses.

The Shallow Shelf Approximation (SSA), loosely referred to as plug flow, was originally proposed to model the flow of ice shelves, assuming negligible vertical shear and basal drag (Morland, 1987; MacAyeal, 1989). But ice shelves are predominantly in a state of longitudinal extension, which can lead to the development of girdle fabrics (Treverrow and others, 2010) that harden the ice for

further along-flow extension compared to isotropic ice (Pimienta and others, 1987; Castelnau and others, 1996). Moreover, horizontal single-maximum fabrics can develop in ice-shelf shear margins similar to those in ice streams (Monz and others, 2021; Thomas and others, 2021; Gerber and others, 2023) that are softer for shear than isotropic ice (Pimienta and others, 1987; Echelmeyer and others, 1994; Jackson and Kamb, 1997), potentially leading to faster interior (trunk) flow than predicted by isotropic ice flow models.

In this paper, we introduce an anisotropic version of the SSA equations and explore the extent to which the viscous anisotropy of ice shelves can be parameterized by an equivalent scalar flow enhancement factor, thus making it easier to account for the effect of crystal fabrics in existing state-of-art isotropic ice flow models. This paper therefore shares the goals of, and builds on, previous work by Ma and others (2010), Placidi and others (2010), Treverrow and others (2015), and Graham and others (2018).

In the following, we begin by introducing the anisotropic SSA problem and how fabric and flow are coupled based on recent developments of ours (Rathmann and Lilien, 2021; Richards and others, 2021; Lilien and others, 2023). The reader less interested in technical details may choose to focus on the subsequent section “Equivalent isotropic enhancement” where we review how the equivalent enhancement factor has previously been defined and introduce a new definition based on the effective strain rate, including testing the validity thereof in an idealized shelf model. Finally, in section “Application to Antarctic ice shelves”, we apply our methodology to the Ross and Pine Island Glacier (PIG) ice shelves, Antarctica, suggesting that significant fabric-induced hardening and softening may be present in Antarctic ice shelves.

## 2. Anisotropic SSA

The coupled problem between ice flow and an evolving crystal fabric requires solving a sequence of consecutive problems (Figure 2). Solving the momentum balance gives the large-scale flow of ice, which, together with the thermal state, allows to drive fabric evolution forward in time. The new fabric state permits an updated estimate of the fabric-induced viscous anisotropy, calculated using a micromechanical model. This, in turn, informs the bulk anisotropic rheology, needed by the momentum balance.

In what follows, we address each of these subproblems before considering how the coupled problem can be simplified by reducing the effect of fabric to a scalar flow enhancement. That is, we seek to simplify the blue and orange pieces in Figure 2 by supposing that the Glen isotropic flow law is sufficient and to test the accuracy of doing so in idealized SSA modeling.

The difference between our approach and previous work that presents the coupled problem in this way (e.g. Gillet-Chaulet and others, 2006; Lilien and others, 2023) is mainly the micromechanical model used and the presumption of plug flow (SSA).

### 2.1. Momentum balance

The SSA momentum balance is (Morland, 1987; MacAyeal, 1989)

$$ \nabla \cdot (H\mathbf{R}) = \rho g H \nabla S, \quad (1) $$

where S is the surface height, H is the ice thickness, $\rho = 917 \text{ kg m}^{-3}$ is the density of ice, $g = 9.8 \text{ m s}^{-2}$ is the gravitational acceleration, and the resistive viscous stress tensor is defined as

$$ \mathbf{R} = \begin{bmatrix} 2\tau_{xx} + \tau_{yy} & \tau_{xy} \\ \tau_{xy} & \tau_{xx} + 2\tau_{yy} \end{bmatrix}. \quad (2) $$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Schematic of the two-way coupled problem between flow and fabric evolution. Anisotropic ice flow modelling requires solving the momentum balance for an anisotropic bulk rheology, the solution of which provides the means to calculate the evolution of the crystal fabric, which in turn allows for an updated estimate of the fabric-induced viscous anisotropy that informs the bulk rheology.</figcaption></figure>

The stress deviator components $\tau_{xx}$, $\tau_{yy}$ and $\tau_{xy}$ follow from the power-law rheology of ice assuming dislocation creep is the dominant deformation mechanism,

$$ \tau = A^{-1/n} \dot{\epsilon}_e^{(1-n)/n} \mathbf{C}, \quad (3) $$

$$ \dot{\epsilon}_e = \sqrt{\frac{\mathbf{C} : \dot{\mathbf{\epsilon}}}{2}}, \quad (4) $$

where A is the depth-averaged flow rate factor, $\dot{\epsilon}$ is the strain rate tensor, $\dot{\epsilon}_e$ is the effective strain rate, n is the flow exponent, and C is the tensorial part of the rheology. Assuming incompressibility, the strain-rate tensor is

$$ \dot{\mathbf{\epsilon}} = \begin{bmatrix} \dot{\epsilon}_{xx} & \dot{\epsilon}_{xy} & 0 \\ \dot{\epsilon}_{xy} & \dot{\epsilon}_{yy} & 0 \\ 0 & 0 & -\dot{\epsilon}_{xx} - \dot{\epsilon}_{yy} \end{bmatrix}, \quad (5) $$

following usual scaling arguments for ice masses with a small aspect ratio, allowing vertical shear to be neglected (MacAyeal, 1989).

The difference between the isotropic and anisotropic SSA equations is presented in Supplementary section A. In effect, vertical velocity gradients vanish only to first order in the ice-mass aspect ratio, whereas for isotropic ice they also vanish to second order. Therefore, care should be taken when applying the anisotropic SSA to ice masses with an aspect ratio less than ~100.

## 2.2. Bulk rheology

In the case of the Glen isotropic rheology, the tensorial part is $\mathbf{C} = \dot{\mathbf{\epsilon}}$. For anisotropic rheologies, however, the strain-rate components should be scaled according to the local fabric-induced viscous anisotropy.

The orthotropic rheology (Gillet-Chaulet and others, 2005) has three mutually perpendicular planes of reflection symmetry with normals $\mathbf{m}_1$, $\mathbf{m}_2$, and $\mathbf{m}_3$, allowing each strain rate component

to be uniquely scaled to account for the effect of fabric development (Figure 2). In the notation of Lilien and others 2023, the tensorial part of the orthotropic rheology is

$$ \mathbf{C} = \sum_{i=1}^{3} (\eta_i \mathbf{P}_i \mathbf{P}_i + \eta_{i+3} \mathbf{P}_{i+3} \mathbf{P}_{i+3}) : \dot{\mathbf{\epsilon}}, \quad (6) $$

where $\mathbf{P}_i = \mathbf{I}/3 - \mathbf{m}_1\mathbf{m}_i$ and $\mathbf{P}_{i+3} = (\mathbf{m}_{j_i}\mathbf{m}_{k_i} + \mathbf{m}_{k_i}\mathbf{m}_{j_i})/2$ are projectors such that $\mathbf{P}_i : \dot{\mathbf{\epsilon}}$ and $\mathbf{P}_{i+3} : \dot{\mathbf{\epsilon}}$ give the longitudinal and shear strain-rate components in the rheological symmetry frame $\mathbf{m}_i$, respectively, and $\mathbf{P}_i\mathbf{P}_i$, $\mathbf{m}_i\mathbf{m}_i$ and $\mathbf{m}_{j_i}\mathbf{m}_{k_i}$ denote outer products. The index tuples are defined as $(j_1, j_2, j_3) = (2, 3, 1)$ and $(k_1, k_2, k_3) = (3, 1, 2)$.

The coefficients $\eta_i$ and $\eta_{i+3}$ are dimensionless directional viscosities that scale the strain rate components in the symmetry frame $\mathbf{m}_i$. These can be conveniently written in terms of the directional enhancement factors $E_{ij}$ (Rathmann and Lilien, 2022)

$$ \eta_i = 3 \frac{E_{j_i j_i}^{2/(n+1)} + E_{k_i k_i}^{2/(n+1)} - E_{ii}^{2/(n+1)}}{\sum_{i=1}^{3} (2E_{j_i j_i}^{2/(n+1)} E_{k_i k_i}^{2/(n+1)} - E_{ii}^{4/(n+1)})}, \quad \eta_{i+3} = 2E_{j_i k_i}^{-2/(n+1)}. \quad (7) $$

Here, $E_{ij}$ are to be understood as the size of the strain rate components in the symmetry frame relative to that if the fabric was isotropic (Glen rheology)

$$ E_{ij} = \frac{\mathbf{m}_i \cdot \dot{\mathbf{\epsilon}}(\hat{\mathbf{\tau}}) \cdot \mathbf{m}_j}{\mathbf{m}_i \cdot \dot{\mathbf{\epsilon}}_{\text{iso}}(\hat{\mathbf{\tau}}) \cdot \mathbf{m}_j}, \quad (8) $$

for a longitudinal ($i=j$) or shear ($i \neq j$) stress state aligned with the symmetry frame

$$ \hat{\mathbf{\tau}} = \tau_0 \begin{cases} \mathbf{I}/3 - \mathbf{m}_i\mathbf{m}_i & \text{if } i = j \\ (\mathbf{m}_i\mathbf{m}_j + \mathbf{m}_j\mathbf{m}_i)/2 & \text{if } i \neq j \end{cases}. \quad (9) $$

In this way, $E_{11}$ is the longitudinal strain-rate enhancement along $\mathbf{m}_1$ when subject to compression/tension along $\mathbf{m}_1$, $E_{12}$ is the $\mathbf{m}_1$-$\mathbf{m}_2$ shear strain-rate enhancement when subject to shear in the $\mathbf{m}_1$-$\mathbf{m}_2$ plane, etc. To be clear: $E_{ij} > 1$ implies a softened material response due to fabric development compared to isotropic ice, whereas $E_{ij} < 1$ implies hardening. Notice that (6)–(7) reduce to the Glen isotropic rheology for $E_{ij} \to 1$.

We mention in passing that modeling ice flow using an orthotropic rheology is not new (Staroszczyk and Gagliardini, 1999; Gillet-Chaulet and others, 2006; Ma and others, 2010; Lilien and others, 2021; 2023) and is capable of (i) reproducing the viscous anisotropy of Dye 3 ice core samples (Rathmann and Lilien 2022; revisited in Supplementary section B), (ii) modeling the fabric profile of the EPICA Dome C ice core (Durand and others, 2007), and (iii) explaining the formation of Raymond bumps with double peaks (Martín and others, 2009).

## 2.3. Fabric evolution

The orientation fabric of ice evolves as a function of several crystallographic processes, such as lattice rotation, continuous dynamic recrystallization (CDRX; or rotation recrystallization), and discontinuous dynamic recrystallization (DDRX), depending on the thermomechanical setting (Figure 3). For brevity, we focus on introducing the models of crystal processes used here, referring the reader to e.g. Montagnat and others (2014) or Faria and others (2014a, 2014b) for an introduction to the physics and observations of crystal processes.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Ice crystal processes affecting orientation fabric development: strain-induced rotation of crystal lattices (lattice rotation) and mass transfer between grains with different orientations (recrystallization; DDRX and CDRX).</figcaption></figure>

In our vertically integrated orthotropic rheology, the fabric field is to be understood as depth-averaged (elaborated on below). We therefore argue that it is sufficient to model the dominant crystallographic process in a depth-averaged sense, which is approximately lattice rotation for cold ice but must include DDRX if warmer than ~-15°C (De La Chapelle and others, 1998; Samyn and others, 2008; Qi and others, 2019).

CDRX accounts for the division of grains along subgrain boundaries resulting from local strain incompatibilities and is generally understood to not change grain orientations much (sub-dominant to lattice rotation), although it may impact the reduction of stored strain energy (Montagnat and Duval, 2000). For this reason, and because some work (Lilien and others, 2023; Richards and others, 2023) suggests little or no need to include this process when trying to reproduce ice core fabrics, we neglect CDRX for simplification.

### 2.3.1. Lattice rotation

Following a popular approach first proposed by Castelnau and others (1996) and Svendsen and Hutter (1996), lattice rotation is modeled as an advection process on the surface of the unit sphere ($S^2$), which is supposed to capture the strain-induced rotation of crystal lattices (c axes). In this model, the c-axis velocity field on $S^2$ depends solely on the bulk stretching $\dot{\epsilon}$ and spin $\omega$ tensors and is given by $\dot{c} = (\omega + \omega_p) \cdot c$, where c is an arbitrary c-axis, $\omega_p = \iota(c^2 \cdot \dot{\epsilon} - \dot{\epsilon} \cdot c^2)$ is the plastic spin, and $c^2$ is the outer product of c with itself. We take $\iota = 1$ to ensure that basal slip systems behave like a deck of cards; that is, slip planes tend to align with the bulk shear plane. Although $\iota$ could take another value to mimic the activity of non-basal slip systems (Gillet-Chaulet and others, 2006; Richards and others, 2021), this phenomenological model is concerned with the activation of the dominant slip system.

Figure 4a-c shows the predicted c-axis velocity field (arrows) and speed (colored contours) for ice subject to three different deformation kinematics (strain geometries) relevant to SSA flows.

### 2.3.2. DDRX

Following Placidi and others (2010), DDRX is modeled as a spontaneous mass decay-production process in orientation space, intended to represent the combined effect of nucleation and grain boundary migration. That is, mass is spontaneously exchanged between grains with different orientations depending on the local stress state, strain rate, and temperature, in a statistical sense.

The decay-production rate is defined as $\Gamma = \Gamma_0(D - \langle D \rangle)$, where the prefactor $\Gamma_0 = \sqrt{\dot{\epsilon} : \dot{\epsilon}}/2A_\Gamma \exp(-Q_\Gamma/RT)$ accounts for the preferential (Arrhenius) activation at warm temperatures and the effect of strain-rate magnitude (Richards and others, 2021;

Lilien and others, 2023). Here, $R$ is the gas constant, $T$ is the temperature, and $A_\Gamma = 1.91 \times 10^7$ and $Q_\Gamma = 3.36 \times 10^4$ J mol$^{-1}$ according to a recent calibration by Lilien and others 2023. The deformability $D$ is the square of the basal-plane resolved shear stress, $(\tau \cdot \tau) : c^2 - \tau : c^4 : \tau$, relative to the average value if the fabric were isotropic

$$D = \frac{(\tau \cdot \tau) : c^2 - \tau : c^4 : \tau}{\tau : \tau/5}, \quad (10)$$

where c is an arbitrary c-axis. Because $D$ is largest for grains with an orientation favorable to basal glide, mass is spontaneously created/added to grains with such preferred orientations (in a statistical sense). Conversely, mass spontaneously decays if $D < \langle D \rangle$, corresponding to grains with an unfavorable orientation being consumed by grains with a more favorable orientation to basal glide. Here, $\langle D \rangle$ is the average grain deformability, where the average $\langle \cdot \rangle$ is defined in Equation (17), and the total mass of the polycrystal is conserved since $\Gamma$ is zero on average, $\langle \Gamma \rangle = \Gamma_0 \langle D - \langle D \rangle \rangle = 0$ (equal amounts of grain mass are created and destroyed per time).

The primary driver for grain boundary migration is generally regarded to be the contrast in dislocation density (stored strain energy) between grains (De La Chapelle and others, 1998; Hamann and others, 2007). In this sense, $1/D$ is a parameterization of the dislocation density: the larger $1/D$ is, the more dislocations a grain has, and the more likely it is to be consumed (decay).

Figure 4d-f shows the normalized decay-production rate $\Gamma/\Gamma_0$ of an isotropic fabric when subject to three different stress states relevant to SSA flows, i.e., mass is transferred from grains with c-axis orientations in red-shaded regions to gray-shaded regions. Notice that dividing by $\Gamma_0$ removes the effect of temperature which otherwise uniformly scales the patterns in Figure 4d-f to increase/decrease the rate of DDRX.

### 2.3.3. Representation

The crystal fabric is represented by the distribution function of grain mass-density in the space of possible c-axis orientations, $\psi$, following the theory of mixtures of continuous diversity (Faria, 2001) (also referred to as $\rho^*$ in the literature, but here $^*$ is reserved for complex conjugation). We expand $\psi$ in terms of a spherical harmonic series following recent work (Rathmann and others, 2021; Richards and others, 2021)

$$\psi(\mathbf{x}, t, \theta, \phi) = \sum_{l=0}^{L} \sum_{m=-l}^{l} \psi_l^m(\mathbf{x}, t) Y_l^m(\theta, \phi), \quad (11)$$

where $\psi_l^m$ are complex-valued expansion coefficients, $Y_l^m$ are harmonic basis functions, and $L$ is the wave mode truncation above which finer-scale structure in $\psi$ is unresolved. Since $\psi$ is the distribution of mass density in orientation space, integrating over
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Fabric dynamics for different deformation kinematics and stress states relevant to SSA flows: c-axis velocity field for lattice rotation (a-c), and DDRX decay-production rate for an isotropic fabric (d-f).</figcaption></figure>

orientation space gives the usual mass density of glacier ice $\rho = \int_{S2} \psi d\Omega$, where $d\Omega = \sin(\theta)d\theta d\phi$ is the infinitesimal solid angle. Note that $\psi_i^{-m} = (-1)^m(\psi_i^m)^*$ by identity since $\psi$ is real, and that the mass orientation distribution function (MODF) is defined as the normalized distribution function $\psi/\rho$.

The combined effect of advection, lattice rotation, and DDRX on $\psi$ is (see e.g. Svendsen and Hutter, 1996; Richards and others, 2021)

$$ \frac{D\psi}{Dt} = -\nabla_{S2} \cdot (\psi\dot{\mathbf{c}}) + \Gamma\psi, \qquad (12) $$

where $\nabla_{S2}$ is the gradient on $S^2$ and $D/Dt$ is the material derivative. The harmonic expansion (11) transforms the fabric model (12) into a high-dimensional matrix problem (Rathmann and Lilien, 2021; Rathmann and others, 2021),

$$ \frac{D\mathbf{s}}{Dt} = (\mathbf{M}_{\text{LROT}} + \mathbf{M}_{\text{DDRX}}) \cdot \mathbf{s}, \qquad (13) $$

where $\mathbf{M}_{\text{LROT}}$ and $\mathbf{M}_{\text{DDRX}}$ are matrices of harmonic interaction coefficients, and

$$ \mathbf{s} = [\psi_0^0, \psi_2^0, \psi_2^1, \psi_2^2, \psi_4^0, \dots, \psi_4^4, \dots, \psi_L^0, \dots, \psi_L^L] \in \mathbb{C}^{(L+2)^2/4} \qquad (14) $$

is the unknown fabric state vector field to be solved for. That is, the harmonic expansion transforms the problem of fabric evolution in $\mathbb{R}^3 \times S^2$ to a high-dimensional advection-reaction vector field problem in $\mathbb{R}^3 \times \mathbb{C}^{(L+2)^2/4}$.

### 2.3.4. SSA fabric model

The vertically averaged model of fabric evolution, consistent with the SSA approximation, is derived in Supplementary section C. Denoting the depth average of $\mathbf{s}(x, y, z, t)$ by $\bar{\mathbf{s}}(x, y, t)$, and implicitly assuming that the velocity field and gradients refer to their horizontal $xy$ parts, $\bar{\mathbf{s}}$ is governed by

$$ \frac{D\bar{\mathbf{s}}}{Dt} = (\mathbf{M}_{\text{LROT}} + \mathbf{M}_{\text{DDRX}}) \cdot \bar{\mathbf{s}} + \frac{a_{\text{sfc}}}{H}(\mathbf{s}_{\text{sfc}} - \bar{\mathbf{s}}) + \frac{a_{\text{sub}}}{H}(\mathbf{s}_{\text{sub}} - \bar{\mathbf{s}}), \qquad (15) $$

where $H$ is the ice thickness (Figure 5). Compared to the above three-dimensional problem, the vertically averaged fabric problem (15) contains additional terms on the right-hand side. The first term is simply the depth-averaged effect of lattice rotation and DDRX. The second and third terms are state-space attractors, causing $\bar{\mathbf{s}}$ to tend towards the characteristic fabric states of ice that accumulates on the surface $\mathbf{s}_{\text{sfc}}$ or subglacially $\mathbf{s}_{\text{sub}}$ (likely isotropic), depending on the positively-defined accumulation rates $a_{\text{sfc}}$ and $a_{\text{sub}}$.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** SSA fabric model and harmonic expansion series of the crystal orientation fabric. SSA fabrics evolve due to the combined effect of fabric advection, englacial crystal processes, and surface/subglacial accumulation of ice. Complex conjugation is denoted by “c.c.”</figcaption></figure>

### 2.3.5. Truncation error

The spectral formulation suffers from a closure problem in the sense that to represent increasing fabric development (increasing alignment of grains), harmonic modes of increasing degree $l$ must be included. For practical reasons, the expansion series is truncated at degree $L=10$ so that the finer scale structure is ignored. To avoid backscatter effects in spectral space at $l=L$ (causing quasi periodic noise to dominate $\psi$), a spectrally sharpened Laplacian term in orientation space is applied (regularization) that disproportionately causes high wave-number modes to decay (for details see Rathmann and Lilien, 2021).

## 2.4. Viscous anisotropy

Given the fabric state $\mathbf{s}$ or its depth average $\bar{\mathbf{s}}$, we are left to provide the local rheological symmetry axes $\mathbf{m}_i$ and associated enhancement factors $E_{ij}$ for the bulk rheology.

We follow Rathmann and Lilien 2021 and Lilien and others 2023 who calculated $E_{ij}$ by substituting $\dot{\epsilon}$ in (8) for the grain-averaged strain-rate tensor, subject to a linear combination of stress (Sachs) and strain-rate (Taylor) homogenization assumptions over the polycrystal scale:

$$E_{ij} = (1 - \alpha) \frac{\mathbf{m}_i \cdot \dot{\epsilon}^{\text{Sachs}}(\bar{\tau}) \cdot \mathbf{m}_j}{\mathbf{m}_i \cdot \dot{\epsilon}_{\text{iso}}^{\text{Sachs}}(\bar{\tau}) \cdot \mathbf{m}_j} + \alpha \frac{\mathbf{m}_i \cdot \dot{\epsilon}^{\text{Taylor}}(\bar{\tau}) \cdot \mathbf{m}_j}{\mathbf{m}_i \cdot \dot{\epsilon}_{\text{iso}}^{\text{Taylor}}(\bar{\tau}) \cdot \mathbf{m}_j}, \quad (16)$$

where $\alpha \in [0; 1]$. The Sachs and Taylor contributions are $\dot{\epsilon}^{\text{Sachs}} = \langle \dot{\epsilon}'(\tau) \rangle$ and $\dot{\epsilon}^{\text{Taylor}} = \langle \tau'(\dot{\epsilon}) \rangle^{-1}$, where $\tau'$ and $\dot{\epsilon}'$ are the microscopic (grain scale) stress and strain-rate tensors, and $\langle \cdot \rangle^{-1}$ inverts the tensorial relationship. The operation for calculating fabric-averaged quantities is defined as

$$\langle \cdot \rangle = \frac{1}{\rho} \int_{S^2} \langle \cdot \rangle \psi \, d\Omega. \quad (17)$$

In effect, grains are modeled as interactionless, which reduce the Sachs and Taylor homogenizations to that of calculating the ensemble-averaged monocrystal rheology.

Here, monocrystals are approximated as linear-viscous transversely isotropic following Rathmann and Lilien 2021:

$$\dot{\epsilon}' = A' \left( \tau' - \frac{E'_{cc} - 1}{2} (\tau' : \mathbf{c}^2) \mathbf{I} + \frac{3(E'_{cc} - 1) - 4(E'_{ca} - 1)}{2} \tau' : (\mathbf{c}^4) + (E'_{ca} - 1)(\tau' : \mathbf{c}^2 + \mathbf{c}^2 : \tau') \right), \quad (18)$$

where $A'$ is the isotropic part of the fluidity, $E'_{ca}A'$ is the fluidity specific for basal plane shear, and $E'_{cc}A'$ is the fluidity specific for compression/tension along $\mathbf{c}$ (Figure 2). The inverse $\tau'(\dot{\epsilon}')$ is simply (18) where $\dot{\epsilon}'$ and $\tau'$ are swapped, $E'_{ij} \to 1/E'_{ij}$, and $A' \to 1/A'$.

The free parameters of our micromechanical model ($A'$, $E'_{cc}$, $E'_{ca}$, $\alpha$) should be chosen such that modeled $E_{ij}$ reproduce bulk deformation tests. In the interest of brevity, we refer the reader to Supplementary section B for details on calibrated values.

From (18) it is clear that $\dot{\epsilon}^{\text{Sachs}}$ and $\dot{\epsilon}^{\text{Taylor}}$ end up depending on the second- and fourth-order structure tensors, $\langle \mathbf{c}^2 \rangle$ and $\langle \mathbf{c}^4 \rangle$. These depend, in turn, exclusively on the low-order harmonics $\psi_2^m$ and $\psi_4^m$ (Rathmann and others, 2022), implying that choosing a small truncation $L$ is sufficient as long as regularization does not affect harmonics of degree $l \le 4$ too much. This is the primary motivation for considering a linear viscous monocrystal rheology in spite of observations that indicate nonlinear viscous behavior (Duval and others, 1983). In the nonlinear-viscous case, the Sachs and Taylor averages depend on higher-order structure tensors (i.e. $l > 4$ harmonics) which is not computationally feasible in our coupled flow-fabric simulations when also needing to dedicate spectral width for regularization, in effect requiring $L \ge 10$ and therefore $\dim(\mathbf{s}) \ge 36$ (Rathmann and others, 2021).

### 2.4.1. Eigenenhancements

The rheological symmetry axes $\mathbf{m}_i$ are taken to align with the fabric principal directions, defined as the eigenvectors of $\langle \mathbf{c}^2 \rangle$ (Gillet-Chaulet and others, 2005). In the following, we therefore refer to the bulk enhancement factors in the frame $\mathbf{m}_i$ as eigenenhancements, denoted by subscripts $i, j = 1, 2, 3$, e.g. $E_{11}, E_{22}, E_{12}$, and so on. These are to be distinguished from the Cartesian enhancements where $\mathbf{m}_i$ are replaced by the Cartesian basis vectors in (8), denoted by subscripts $i, j = x, y, z$, e.g. $E_{xx}, E_{yy}, E_{xy}$, and so on, primarily useful for diagnostic purposes.

## 3. Equivalent isotropic enhancement

Following previous work (e.g. Ma and others, 2010; Placidi and others, 2010; Treverrow and others, 2015; Graham and others, 2018), we posit that the effect of fabric on ice viscosity can be approximated by an equivalent isotropic enhancement factor, understood in the usual sense as a scalar prefactor $E$ to the flow rate factor $A$ in the Glen rheology (3). If additional softening/hardening effects are present, the combined effect is given by the product of the respective enhancement factors. For example, damaged (fractured) ice is softer than undamaged, and dusty ice from glacial periods is typically twice as soft as interglacial ice (Dahl-Jensen and Gundestrup, 1987) (although this effect is to some extent indirectly
<figure><img src="image_7.png" /><figcaption>Figure 6. Biases resulting from approximating the viscous anisotropy of ice using a scalar enhancement factor. (a): Normalized strain-rate components of the orthotropic (black lines) and isotropic (grey and colored lines) rheology when subject to a fixed xy shear stress that is increasingly unfavorably aligned with a horizontal single-maximum fabric (decreasing compatibility). (b): Same as (a) but for a fixed horizontal single-maximum fabric aligned with the y axis, subject to a stress state that varies linearly between xy shear and uniaxial tension along y (varying stress superposition). Colored lines show predictions for the Glen rheology when using either CAFFE (purple) or EIE (green) to calculate E.</figcaption></figure>

caused by the dust impeding DDRX; Durand and others (2007)), implying

$$E = E_{\text{fabric}} E_{\text{damage}} E_{\text{dust}} \dots \quad (19)$$

In this work, the role of fabric is considered in isolation, so $E = E_{\text{fabric}}$ is henceforth assumed implicit.

Before introducing specific models of E, it is helpful to quantify two important biases that can occur when approximating the effect of fabric as a scalar flow enhancement. In short, we first set up a horizontal single-maximum fabric and apply a fixed horizontal shear stress, then run a series of controlled experiments, gradually rotating the orientation of the single maximum from perfectly aligned with the shear stress direction to misaligned and determine how closely a scalar enhancement of the Glen rheology can approximate the orthotropic strain rate components. Second, using the same horizontal single-maximum fabric, we gradually shift the applied stress from being horizontal simple shear, perfectly aligned with the single maximum, to being extensional in the horizontal plane and again determine how well a scalar enhancement factor performs.

### 3.1. Biases

The claim that viscous anisotropy can be represented by (reduced to) a scalar enhancement E implies that some viscosity tensor components matter more than others. This is not unreasonable: if ice is (say) subject to xy shear stress, setting $E = E_{xy}$ will probably give the most accurate strain-rate prediction, whereas the remaining directional enhancements matter less insofar as the fabric is favorably aligned to the shear stress, $\mathbf{m}_1, \mathbf{m}_2, \mathbf{m}_3 = \hat{\mathbf{x}}, \hat{\mathbf{y}}, \hat{\mathbf{z}}$ (or some permutation thereof). Figure 6a makes this point more clear by showing $\dot{\epsilon}_{xy}$ and $\dot{\epsilon}_{yy}$ (black lines) of the orthotropic rheology

(using our micromechanical model and $n = 3$), normalized by the effective isotropic strain rate

$$\dot{\epsilon}_{e,\text{iso}} = \sqrt{\frac{\dot{\epsilon} : \dot{\epsilon}}{2}}, \quad (20)$$

for a strong horizontal single-maximum fabric that is increasingly unfavorably aligned to (rotated away from) a fixed xy shear stress, denoted by the angle $\vartheta$ between the fabric principal direction $\mathbf{m}_1$ and direction of shear. Note that both the stress magnitude and the rate factor A cancel out in the division. As anticipated, the orthotropic xy strain rate (solid black line) is found to decrease from $E_{12}$ when aligned favorably at $\vartheta = 0^{\circ}$ (softer than predicted by the Glen rheology; solid gray line) to $\sim E_{11}$ when unfavorably aligned at $\vartheta = 45^{\circ}$ (harder than predicted by the Glen rheology). However, using a scalar enhancement E to account for the effect of fabric is complicated by the fact that non-coaxial behavior is permitted by the orthotropic rheology, unlike the Glen isotropic rheology where $\tau$ and $\dot{\epsilon}$ are coaxial. In the present example of simple shear, coaxiality is fulfilled only when the fabric and stress states are favorably aligned ($\vartheta = 0^{\circ}$), in which case setting $E = E_{12}$ allows the Glen rheology to reproduce the orthotropic rheology and thus give the sought strain rate tensor. But as the misalignment increases, other strain-rate components may become nonzero (e.g. $\epsilon_{yy}$; dashed black line) which a scalar enhancement of the Glen rheology cannot capture (dashed grey line). The degree to which a scalar enhancement E can accurately represent fabric-induced softening/hardening depends, therefore, partly on the alignment between the local fabric and the stress (or strain rate) geometry, henceforth referred to as fabric compatibility.

In addition to fabric compatibility, another bias may occur whenever $\tau$ is a superposition of two or more stress states. Consider
a horizontal single-maximum fabric, perfectly aligned with the y-axis. Such an ice fabric is soft for xy shear but hard for compression/tension along x and y. Hence, it is not clear which enhancement-factor component $E_{ij}$ to apply if shear and longitudinal stresses are superimposed—is the ice soft or hard? Figure 6b makes this point clear by showing $\dot{\epsilon}_{xy}$ and $\dot{\epsilon}_{yy}$ of the orthotropic rheology (black lines), normalized by $\dot{\epsilon}_{e,iso}$, for a stress state that varies linearly between xy shear and uniaxial tension along y. Scaling all strain rate components in the Glen rheology (gray lines) by the same factor, E, to match the orthotropic response is not possible when both shear and tension are simultaneously important. The degree to which E can accurately represent fabric-induced softening/hardening depends, therefore, also on the number and strength of stress (or strain rate) modes that contribute to $\tau$ (or $\dot{\epsilon}$), henceforth referred to as stress (or strain-rate) superposition.

In what follows, we model E using the CAFFE method (Continuum-mechanical, Anisotropic Flow model, based on an anisotropic Flow Enhancement factor; Placidi and others (2010)) together with a new method that we find useful based on the effective strain rate. Other methods such as ESTAR (Empirical Scalar Tertiary Anisotropy Regime; Graham and others (2018); McCormack and others (2022)) also exist, but are not applied in our work due to strong assumptions about fabric compatibility that we wish to relax. The interested reader is referred to Supplementary section D for further details on ESTAR.

## 3.2. CAFFE

Placidi and others (2010) proposed parameterizing an equivalent isotropic enhancement in terms of the average grain deformability $\langle D \rangle$. Isotropic fabrics are then characterized by $\langle D \rangle = 1$, while fabrics that are minimally and maximally favorable to basal glide are characterized by $\langle D \rangle = 0$ and $\langle D \rangle = 5/2$, respectively. CAFFE assigns the bulk enhancements $E_c$ and $E_s$ to these minimally and maximally favorable states, respectively, whereas for intermediate states a semi-empirical scaling is proposed based on Azuma (1995):

$$E = \begin{cases} \frac{(1 - E_c)\langle D \rangle^{8/21(E_s-1)/(1-E_c)} + E_c}{4\langle D \rangle^2(E_s - 1) + 25 - 4E_s} & \text{for } 0 \le \langle D \rangle \le 1, \\ \frac{21}{4\langle D \rangle^2(E_s - 1) + 25 - 4E_s} & \text{for } 1 < \langle D \rangle \le 5/2. \end{cases} \quad (21)$$

Here, we follow Placidi and others (2010) by setting $E_c = 0.1$ and $E_s = 10$, in agreement with bulk deformation tests made on strong single-maximum fabrics subject to compression or shear that is aligned with the preferred c-axis direction, respectively. Because CAFFE assumes that the Glen isotropic rheology is a suitable description at large scale, $\tau$ can be exchanged for $\dot{\epsilon}$ when calculating the deformability D.

## 3.3. Our definition (EIE)

Formally, E is the strain rate response of a parcel of anisotropic ice, divided by that of isotropic ice, for a given stress state. This ratio is only meaningful for scalar quantities, so we conjecture that the effective strain rate $\dot{\epsilon}_e$ in (4) is a useful quantity to consider, which fulfills the effective Norton-Bailey creep law (Naumenko and Altenbach, 2007, p. 21) for both isotropic and orthotropic rheologies:

$$\dot{\epsilon}_e = A\tau_e^n, \quad (22)$$

where $\tau_e$ is the effective deviatoric stress. The motivation for considering $\dot{\epsilon}_e$ is straight forward: it is a natural measure (contraction) of the strain-rate tensor that weighs each component according to $E_{ij}$. Further, $\dot{\epsilon}_e$ satisfies the requirement that if deformation is (say) simple shear aligned with $m_2$ and $m_3$ so that $\dot{\epsilon} \propto P_4$, then $\dot{\epsilon}_e$ is dominated by the term involving the corresponding shear enhancement $E_{23}$ since $P_4 : \dot{\epsilon}$ is the only nonzero strain-rate projection in (6).

The problem is, however, not as easy as defining $E = \dot{\epsilon}_e/\dot{\epsilon}_{e,iso}$, where $\dot{\epsilon}_e$ and $\dot{\epsilon}_{e,iso}$ are given by the contraction (4) with C replaced by (6) and $\dot{\epsilon}$, respectively. In the aforementioned case of simple shear, it can be shown that

$$\frac{\dot{\epsilon}_e}{\dot{\epsilon}_{e,iso}} = \sqrt{\frac{\eta_4(\mathbf{P}_4 : \dot{\mathbf{e}})^2}{\dot{\mathbf{e}} : \dot{\mathbf{e}}}} = E_{23}^{-1/(n+1)}. \quad (23)$$

Hence $\dot{\epsilon}_e/\dot{\epsilon}_{e,iso}$ should be raised to the power of $-n-1$ to give the desired result $E = E_{23}$. We therefore propose defining E according to the scaling relation (henceforth Equivalent Isotropic Enhancement, or EIE)

$$E = \left( \frac{\dot{\epsilon}_e}{\dot{\epsilon}_{e,iso}} \right)^q, \quad (24)$$

where $q = -n-1$. The fact that q is negative can also be understood intuitively: if an observed flow field is (say) soft for shear due to fabric, the stress required to facilitate that deformation is less than that required for isotropic ice and therefore $\tau_e < \tau_{e,iso} \Rightarrow \dot{\epsilon}_e < \dot{\epsilon}_{e,iso}$ by virtue of (22).

In the general case where deformation is a superposition of different kinematic modes, or if fabric compatibility is poor, $\dot{\epsilon}_e$ will consist of multiple contributions and the appropriate exponent q might differ from $-n-1$. However, in the following comparison between CAFFE and EIE, we find that $q = -n-1$ is a good approximation and is therefore taken to suffice.

## 3.4. Method comparison

Figure 7 shows E predicted by CAFFE and EIE (colored lines) for a parcel of initially-isotropic ice subject to different deformation kinematics relevant to SSA flows, assuming that DDRX is either negligible (cold ice limit; Figure 7a-c), strong (warm ice limit; Figure 7d-f), or very strong (very warm ice limit; Figure 7g-i), imposed by setting $\Gamma_0 = 0$, $\Gamma_0 = \Gamma_0(-15^{\circ}\text{C})$, and $\Gamma_0 = \Gamma_0(-5^{\circ}\text{C})$, respectively. In general, EIE is found to approximately recover the most important component of $E_{ij}$ (black line). CAFFE predicts similar behavior, but with more significant magnitudes. However, setting $E_c = 0.2$ and $E_s = 3.2$ enables CAFFE to approximately reproduce the behavior of EIE, henceforth referred to as CAFFE† (dashed purple line). Note that modeled rates of DDRX depend exponentially on temperature, so steady-state fabrics are quickly reached for very warm ice, which fulfills the compatibility assumption of ESTAR (Supplementary section D).

The discrepancy between EIE and CAFFE is partly due to spectral resolution (increasing the truncation L permits larger E for fabrics of intermediate strength; not shown) and partly due to the use of a linear grain rheology rather than nonlinear (Rathmann and others, 2021). Insofar as CAFFE is able to correctly estimate the degree of hardening and softening, adopting CAFFE over EIE or ESTAR is preferable since it allows for both an evolving fabric and is able to predict strong enhancements not limited by grain-rheological assumptions.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Fabric-induced enhancement factors for different deformation kinematics relevant to SSA flows, depending on whether DDRX is negligible (cold ice limit; panels a-c), strong (warm ice limit; panels d-f) or very strong (very warm ice limit; panels g-i). In each panel, the equivalent enhancement *E* is shown for each method (colored lines) compared to the most relevant component of *E*<sub>ij</sub> (black line). MODF insets show the modeled fabrics at selected strains.</figcaption></figure>

Regarding the biases mentioned above, caused by fabric incom-
patibility and strain-rate superposition, none of these are relevant
in Figure 7 except in the case of simple shear where the fabric
and strain-rate principal frames are slightly misaligned (Figure
7a). To quantify how these biases affect CAFFE and EIE, Figures
6a,b therefore also show the strain rate components predicted
by the Glen rheology using either method (purple and green

lines, respectively). Overall, both CAFFE and EIE demonstrate
a reasonable ability to capture the orthotropic shear enhance-
ment (solid lines approximately agree in Figure 6), but signif-
icantly overestimate longitudinal strain rates when the stress
state is a superposition of shear and tension (dashed pur-
ple and green lines disagree with the dashed black line in
Figure 6b).
<figure><img src="image_9.png" /><figcaption>Figure 8. Schematic of the idealized half-width ice shelf model.</figcaption></figure>

## 3.5. Accuracy of assuming isotropic viscosity

Although the simple comparison above provides some confidence in the ability of CAFFE and EIE to parameterize the most important component of viscous anisotropy as a scalar enhancement, the error made by neglecting the tensorial viscosity structure in more realistic flow problems remains to be understood. Before proceeding to calculate E over Antarctic ice shelves, we therefore estimate the velocity error resulting from replacing $E_{ij}$ with E in a simple ice shelf model, thus aiming to quantify (albeit not comprehensively) potential flow biases introduced by neglecting the tensorial viscosity structure. More precisely, we tested how well the isotropic SSA rheology can reproduce the steady state velocities of an idealized ice shelf model that builds on the orthotropic SSA rheology (i.e., full fabric-flow coupling) when calculating E using CAFFE and EIE.

### 3.5.1. Model setup

The orthotropic model considers the transient evolution of a laterally-confined shelf of half width $W = 15$ km, length $L = 100$ km, and uniform initial thickness $H_0 = 500$ m, subject to no-slip conditions at its lateral boundaries (Figure 8). The ice thickness $H$ is allowed to evolve according to the local mass divergence plus a Laplacian term for regularization but neglects accumulation and subshelf melting,

$$ \frac{\partial H}{\partial t} = -\nabla \cdot (H\mathbf{u}) + \beta_H \nabla^2 H, \quad (25) $$

where $\mathbf{u}$ is the velocity field and $\beta_H = 1 \times 10^{-4}$ is the strength of Laplacian regularization (chosen to be as small as possible, yet allowing for numerical convergence). The thickness is taken to be constant on the inflow boundary, equal to the initial thickness $H_0$. On the outflow boundary, calving is assumed instantaneous so that the shelf length is preserved, and is subject to the usual stress boundary condition according to the seawater pressure gradient (Huth and others, 2021) where the freeboard height is determined from the floatation criterion.

The ice mass is taken to be isothermal with a temperature of $T = -15^\circ C$ to make the effect of fabric on modeled viscosity unambiguous. The canonical rate factor $A(T)$ by Cuffey and Paterson (2010, p. 72–74) is used, but we recognize that recent work (Fan and others, 2025) provides a welcome set of updated flow parameters for both grain size sensitive and insensitive deformation. The bulk flow law exponent $n$ is set equal to 3 following Cuffey and Paterson (2010, p. 72–74), a popular choice in large-scale flow modeling.

This exponent is supposed to capture the combined effect of grain boundary sliding and dislocation creep in a single power law (Behn and others, 2021) despite recent work suggesting that $n \simeq 4$ (Bons and others, 2018; Ranganathan and Minchew, 2024; Fan and others, 2025) or $n < 3$ (Wang and others, 2025c) might be better suited. Changing $n$ in our box model should affect only the degree of strain-softening and hence the shape of the flow profile; that is, the fabric type that develops in the shear margin and the trunk should remain unchanged, although the strength might be different.

The crystal fabric on the inflow boundary is set to be isotropic but is unconstrained elsewhere. Fabric evolution is assumed to be dominated by the depth-averaged effect of crystal processes; that is, dominated by englacial fabric evolution owing to lattice rotation and DDRX. We therefore neglect contributions to (15) from surface and subglacial accumulation that add new ice with a different fabric signature. An additional real-space Laplacian term $\beta_s \nabla^2 \mathbf{s}$ is added to (15) for numerical stability following Rathmann and Lilien 2021 and Lilien and others 2023 at the expense of slightly limiting how large spatial fabric gradients are permitted (more sophisticated upwinding techniques were not attempted). The strength of regularization $\beta_s = 1 \times 10^{-5}$ was chosen to be as small as possible, yet to allow numerical convergence at each time step. A spectral resolution of $L = 10$ was chosen to be consistent with the above Lagrangian parcel modeling.

### 3.5.2. Numerics

Both the CAFFE and the EIE model can be used directly as a prefactor to $A$ in the Glen rheology, dynamically calculated at each time step given the fabric state of the previous time step. If solved fully implicit w.r.t. $\dot{\epsilon}$, this changes, however, the nonlinearity of the momentum balance that may require modifying to the Stokes solver in existing ice flow models. In practice, $E$ can instead be calculated from the strain rate of the previous time step, similar to how we solve the coupled anisotropic problem (Figure 2), and is the approach taken here.

The whole problem was solved using FEniCS (Logg and others, 2012), an open source software to solve PDEs using the finite element method. Our implementation relies on Newton’s method to solve the nonlinear momentum balance. The Jacobian of the residual form (required for Newton iterations) was calculated using the unified form language (UFL) (Alnæs and others, 2014), which is used by FEniCS to specify weak forms and supports automatic symbolic differentiation. For brevity, the reader is referred to Rathmann and Lilien 2021 for details on the weak form of the equations solved.

### 3.5.3. Results

The steady state of the orthotropic model is shown in Figures 9a–d and 10a–d for DDRX being negligible $\Gamma_0 = 0$ or strong $\Gamma_0 = \Gamma_0(-15^\circ C)$, respectively. As expected, the shelf is found to thin and accelerate towards the calving front in both cases, and shear margins develop around the trunk that otherwise deforms mostly by extension.

When DDRX is negligible (Figure 9), significant anisotropy develops in the shear margins (horizontal single-maximum fabrics) as measured by the pole figure index $J = \int_{S_2} |\psi|^2 d\Omega$ (Ismail and Mainprice, 1998) (Figure 9b; a value of 1 represents isotropic ice, larger values represent increasing grain alignment), while girdle fabrics develop in the trunk that become particularly strong near the calving front. Modeled MODFs are shown as insets for selected locations, denoted by markers 1–6. When DDRX is strong
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Orthotropic model results in steady state when lattice rotation is the dominant crystal process. (a): Ice speed (colored contours) and thickness (white contours). (b): Strength of fabric anisotropy as measured by the pole figure *J* index. (c) and (d): Fabric-induced shear and longitudinal enhancement factors. (e) and (f): *E* predicted by EIE and corresponding velocity misfit, respectively. (g) and (h): Same as EIE panels but for the CAFFE† method. Isotropic and free fabric boundaries are shown as cyan and magenta lines, respectively. Examples of MODFs are shown at selected locations denoted by markers 1–6. Dashed contours in panel (f) and (h) show the velocity misfits resulting from entirely disregarding the effect of fabric (naively applying the Glen rheology with *E* = 1 for the steady-state ice geometry).</figcaption></figure>
Steady state orthotropic model — Lattice rotation + DDRX

Warm ice, T = -15°C

<figure><img src="image_11.png" /></figure>

EIE

<figure><img src="image_12.png" /></figure>

CAFFE†

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Orthotropic model results in steady state when DDRX is the dominant crystal process. Caption same as in Figure 9.</figcaption></figure>

(Figure 10), the horizontal anisotropy and induced softening that develops in the shear margins is more pronounced than if DDRX is negligible. In the trunk where flow is primarily extensional, a double maximum structure is found, making along-flow extension softer than if DDRX is negligible.

In both orthotropic models, simulated MODFs roughly match the patterns anticipated from the idealized deformation experiments shown in Figure 7.

Figures 9c-d and 10c-d show the Cartesian enhancements expected to be most important for the simulated flows (Figure 8).
In the model without DDRX, both EIE and CAFFE† are found to approximate $E_{xy}$ in the shear margin but $E_{xx}$ in the trunk (Figure 9e,g), resulting in velocity errors that are generally below 10% (Figure 9f,h; see Supplementary section F for actual velocity maps). The same applies to the model with strong DDRX, except in the trunk where hardening is somewhat underestimated and that CAFFE† predicts slightly too small $E$ in the shear margin, which explains why the upstream flow speeds are too slow by more than 15% (Figure 10f versus 10h). This bias is likely caused by CAFFE† underestimating $E$ for DDRX fabrics in simple shear (Figure 7d,g). Whether there are more optimal values of $E_c$ and $E_s$ that improve the CAFFE† misfit for DDRX fabrics was not explored.

Note that we considered CAFFE†, rather than CAFFE (no dager), for these experiments, since comparing misfits is meaningful only if $E_c$ and $E_s$ are chosen such that the enhancement magnitudes of the orthotropic SSA rheology are approximately reproduced, and thus overall flow speed. Note also that if the effect of fabric is neglected altogether by setting $E = 1$, velocity errors as large as 30% to 60% are found (dashed contours in Figures 9f,h and 10f,h).

### 3.5.4. Summary

In summary, combining the Glen isotropic rheology with either CAFFE† or EIE is able to reproduce the behavior of the orthotropic rheology with acceptable velocity errors (at least in our idealized model), suggesting that either model of $E$ is potentially useful for understanding the effect of fabric development on real ice shelves. However, when applied to real ice shelves (such as in the following), the original CAFFE model is preferred; it allows the full range of observed fabric hardening and softening to be realized, unaffected by asymptotic behavior like that in our grain average model on which EIE depends (see Figure 7 and Supplementary section B).

Although $E$ is found to map to the components of $E_{ij}$ largely following intuition built from the above Lagrangian parcel experiments, we warn that such intuition might fail if advection is important so that fabric incompatibility is significant or if recent ice-dynamical changes cause a significant superposition of stresses (see the above discussion on biases).

## 4. Application to Antarctic ice shelves

We investigate the importance of accounting for fabric-induced viscous heterogeneities over ice shelves by using CAFFE to calculate the equivalent enhancement $E$ over the Ross and Pine Island Glacier (PIG) ice shelves, Antarctica. Both ice shelves and their upstream glaciers are presumed to exhibit plug flow such that the local strain rate tensor can be determined from satellite-derived surface velocities (MEaSURES; Rignot and others 2017). To close the SSA fabric problem, we further assume that all flows are in approximately steady state and therefore fabric is, too, $\partial\bar{s}/\partial t = 0$. Since it is unclear where DDRX is important (observations of fabric and temperature are few and sparse), we calculate the steady-state fabric field assuming that DDRX is either negligible or strong to explore the range of possible predictions for $E$, which are supposed to represent the limits of cold and warm ice, respectively. Uncertainties and caveats introduced by the above assumptions are addressed in the discussion.

If lattice rotation is the dominant crystal process (no DDRX), fabric evolution (15) reduces to a boundary value problem w.r.t. $\bar{s}$.

that can be solved given surface velocities alone

$$ (\mathbf{u} \cdot \nabla)\bar{s} = \mathbf{M}_{\text{LROT}} \cdot \bar{s} + \beta_s \nabla^2 \bar{s}. \quad (26) $$

That is, no velocity field modeling is required, hence nothing is assumed about the ice rheology, ice geometry, or boundary conditions, other than that SSA is valid. Disregarding regularization, the fabric field $\bar{s}$ that solves (26) is simply a balance between advective effects and lattice rotation. If velocity gradients are small, then fabric is advected along flow lines (long memory timescale), whereas if velocity gradients are large, then upstream fabrics matter less compared to the local strain-induced rotation of c-axes (short memory timescale).

If DDRX is non-negligible, then $\mathbf{M}_{\text{LROT}} \rightarrow \mathbf{M}_{\text{LROT}} + \mathbf{M}_{\text{DDRX}}$ in (26) and fabric evolution depends on $\tau$ and $T$ in addition to $\dot{\epsilon}$. To close the DDRX problem, we approximate $\tau$ as coaxial to $\dot{\epsilon}$ and assume isothermal ice at $T = -15^{\circ}\text{C}$ over both Ross and PIG. Coaxiality is not guaranteed unless fabric is favorably aligned with the local strain rate geometry (i.e., advective effects are small and no recent ice-dynamical changes have occurred; see above discussion). Temperatures as high as $-15^{\circ}\text{C}$ are probably unlikely over large parts of Ross and PIG, especially in upstream regions where the ice is not in contact with the warm ocean. But we are primarily interested in the potential effect of DDRX, which at $-15^{\circ}\text{C}$ is strong enough that developed fabrics can differ substantially from those generated by lattice rotation alone (Figure 7). If temperatures were better known over Ross or PIG, the importance of DDRX in controlling fabric fields (and thus $E$) could be determined more precisely. Without this information, the best we can do is to provide end-member solutions where the ice masses are assumed to be either very cold or very warm so that DDRX is negligible or strong, respectively.

The steady-state Ross and PIG fabric problems were solved numerically using FEniCS, and the method was validated by reproducing the steady-state fabrics of the idealized shelf model above (see Supplementary section E). For both Ross and PIG, we considered a rectangular bounding model geometry with a mesh resolution of approximately 8 km and 2.5 km, respectively, and assumed isotropic ice on upstream inflow boundaries but unconstrained on outflow boundaries. Because lattice rotation and DDRX are modeled to be more important than advection over grounded ice, assuming isotropic ice on inflow boundaries affects the solution only in the neighborhood near the inflow boundaries. That is, our method and results are robust against uncertainties in fabric boundary conditions if the model boundaries are set to be sufficiently far from the region of interest (see Supplementary section G).

### 4.1. Results

Figures 11 and 12 show our model results for Ross and PIG, respectively, including surface speeds in panel a and the corresponding strain rate magnitude $\dot{\epsilon}_{e,iso}$ in panel b.

The strength of fabric horizontal anisotropy is shown in panels d and f for the model with and without DDRX, defined as the horizontal eigenvalue difference $\Delta\lambda$ of $\langle c^2 \rangle$, a measure frequently reported in radioglaciological field studies (e.g. Gerber and others, 2023; Zeising and others, 2023). In the model of Ross and PIG that neglects DDRX (lattice rotation only), horizontal anisotropy is found to be well developed throughout, and horizontal single-maximum fabrics that are (sub)perpendicular to the shear plane are prevalent in shear zones. In contrast, girdle-like
<figure><img src="image_14.png" /><figcaption>Lattice rotation</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Cold ice</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Lattice rotation + DDRX</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Warm ice, T = -15 °C</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 11.** SSA fabric model results for the Ross ice shelf. (a): Satellite-derived surface velocities. (b): Effective strain rate. (c): *E* estimated using CAFFE assuming lattice rotation is the dominant crystal process (cold ice limit). (d): Fabric horizontal eigenvalue difference. (e,f): Same as (c,d) but assuming DDRX is strong (warm ice limit). Isotropic and free model boundaries are shown as cyan and magenta lines, respectively, and floating ice is delineated by green contours. Modeled MODFs are shown at selected locations, denoted by markers 1–4.</figcaption></figure>

fabrics are modeled in regions of longitudinal or trunk flow (see MODF insets). In the DDRX-activated model of Ross and PIG, horizontal anisotropy is also well developed, especially in shear margins where $\Delta\lambda$ is much larger compared to the model that neglects DDRX.

The CAFFE enhancement *E* is shown in panels *c* and *e*. In the case of the Ross ice shelf, we find that lattice rotation alone can

give rise to widespread hardening of the shelf (*E* < 1), but if DDRX is strong then zones of softening are also present (*E* > 1), especially near the grounding line and in shear margins. For PIG, both models with and without DDRX suggest that the shelf and its near-upstream region are softened due to fabric, particularly the shear margins, while hardening is prevalent over some grounded parts of PIG if DDRX is negligible.
<figure><img src="image_19.png" /><figcaption>**Figure 12.** PIG model results. Caption same as in Figure 11.</figcaption></figure>

# 5. Discussion

## 5.1. Ross and PIG results

The fact that fabric development results in hard ice over most of Ross when lattice rotation is dominant, but that significant softening can occur when DDRX is strong, leaves us with little ability to

conclude how fabric might affect the shelf viscosity without knowing more about the relative importance of the two processes. In contrast, fabric development over PIG seems to generally soften the shelf and upstream shear margins regardless of the strength of DDRX. PIG is a dynamically active part of West Antarctica, which has contributed ~ 1.5 mm to global sea level between 1995 and
2017 (Shepherd and others, 2019) and may undergo irreversible retreat should its grounding line continue to migrate further backward (Favier and others, 2014). Accounting for the effect of fabric development over PIG in large-scale flow models (particularly the effect of shear-margin weakening) might therefore be important for accurately predicting the near-term mass loss from Antarctica, as well as other key outlet systems in Antarctica and Greenland.

To make further progress on understanding the effect of crys-
tal fabrics on the viscosity structure of ice shelves, we emphasize
the need for more field observations of fabric and temperature to
help inform and validate fabric modeling. Useful locations would
include transects on ice shelves where fabric gradients are expected
to be large, such as across shear margins, along flow lines with
extensional flow, and near pinning points. Only with widespread
ice temperature measurements (allowing for the rate of DDRX
to be accurately modeled) and widespread fabric measurements
(allowing modeled fabric fields to be validated) do we have con-
fidence in relaxing the “cold” and “warm” ice bounds considered
here. Ultimately, this is needed for large-scale ice flow models to
accurately account for the effect of fabric development in prognos-
tic simulations of future ice mass loss.

5.1.1. Steady state assumption

The Ross and PIG ice shelves have experienced large area changes of -1.3% and -16.6%, respectively, between 2009 and 2019 (Andreasen and others, 2023), casting doubt on whether crystal fabrics can indeed be predicted assuming steady flow as is done here (assuming area changes are not caused by melting alone). Moreover, in the grounded parts of the Western Ross, there is evidence of deceleration on Whillans ice stream (Joughin and others, 2002), that Kamb ice stream might have shut down ~ 150 to 170 years ago (Retzlaff and Bentley, 1993; Catania and others, 2006), and that Bindschadler ice stream shut down ~ 450 years ago despite flowing rapidly today (Conway and others, 2002). Determining how these dynamical changes affect fabric fields, and in turn flow, is a major challenge that only forward ice flow modeling can answer when combined with field measurements (e.g. Leung and others, 2025). We emphasize, however, that our aim is *not* to provide accurate present-day maps of *E* over Ross and PIG, directly useful for large-scale ice flow modeling. Rather, we wish to highlight the potential effect that developed fabrics can have on controlling ice viscosity in regions of floating ice, which may constitute a large rheological uncertainty when attempting to estimate future ice mass loss.

5.1.2. Fabric compatibility and strain rate superposition

We evaluated the extent to which fabric incompatibility and strain-rate superposition could make our estimates of *E* less applicable to the Ross and PIG systems by calculating the shear fraction (Graham and others, 2018) and a new measure of fabric compatibility. We refer the reader to Supplementary section H for details, but note that fabric compatibility and strain-rate superposition over Ross and PIG is no worse than in our idealized shelf model where velocity errors are relatively small when using CAFFE (Figures 9 and 10). We therefore argue that our maps of *E* over Ross and PIG are sufficiently unaffected by these biases to be useful for understanding the potential effect of fabric on ice viscosity.

5.1.3. Ice damage

As indicated in Figure 1, strongly damaged (fractured) ice is much
softer than undamaged. In regions where damage is prevalent, the

main rheological uncertainty is probably the ability to model ice
damage and its effect on ice viscosity. Between 2017 to 2020, three
large calving events took place at PIG, and the resulting veloc-
ity speedup (Joughin and others, 2021) has been suggested to be
partly facilitated by the damage that evolved in the ice-shelf shear
margins (Sun and Gudmundsson, 2023). On the other hand, more
recent work (Gerli and others, 2024) finds that remotely sensed
damage maps over the PIG and Filchner–Ronne shelves are only
modestly related to patterns in the inferred flow rate factor A,
suggesting that remotely detected damage fields are not of direct
relevance to present-day ice shelf flow. Either way, if the effect of
fabric on ice viscosity is multiplicative as assumed here (19), any
damaged-induced softening should be further modified (multi-
plied) by the effect of fabric; that is, the fact that ice is damaged does
not cause fabric development to become irrelevant. Indeed, recent
modeling has suggested that recrystallization (causing larger grain
sizes in addition to fabric development) can increase the vulnera-
bility to fracture of ice shelves by decreasing the tensile strength
of shear margins by up to ~ 75% (Ranganathan and others,
2021).

5.1.4. Suture zones

Suture zones, defined as the boundaries between provinces with different upstream sources of ice, are believed to impede fracture propagation, possibly due to elevated seawater content (Kulessa and others, 2019), thus enhancing ice-shelf stability. The Ross ice shelf is fed by several different ice streams, and the long narrow stipes with large values of Δλ in Figure 11d seem to collocate with the suture zones identified by LeDoux and others (2017). Recent work suggests that these suture zones are regions of low horizontal viscosity due to fabric or aligned crevasses (Wang and others, 2025c). Our results indicate that the effect of fabric is predominantly to harden the ice in these regions (although DDRX might to some extent soften them; Figure 11e), suggesting that low-viscosity suture zones on Ross are possibly the result of aligned crevasses.

## 5.2. Implications for large-scale modeling

The agreement between velocities modeled using the orthotropic, EIE, and CAFFE† rheologies suggests that a scalar, isotropic enhancement can sufficiently approximate the full orthotropic rheology if the fabric field is known. Typical ice flow models, including those that solve vertically-integrated stress balance problems like the SSA, already permit spatially varying fields of A, so incorporating the EIE or CAFFE methods into large-scale models is trivial if the fabric is known. However, the sparsity of fabric measurements limits the use of fabric as a constraint on ice flow models and validation of fabric development models that could be used as a stopgap. For example, some models that rely on laboratory-calibrated rates of DDRX struggle to reproduce the fabric development observed in ice cores (Lilien and others, 2023; Richards and others, 2023), suggesting care must be taken to ensure that the rate of DDRX is modeled correctly. Considering the variation in enhancement that can result from DDRX being weak or strong (e.g. over Ross; Figure 11), this uncertainty in fabric directly impacts the putative effect on ice flow. In certain situations, the flow enhancement appears independent of the process dominating fabric development (e.g. over PIG; Figure 12), implying that, in isolated areas, there may be sufficient confidence in the effect of fabric on flow to include it in large-scale models.
Our results show that with better validation of fabric models, or better spatial constraints on fabric, the effect of fabric could be included in existing large-scale, transient models of ice shelf flow without much effort or computational expense. In practice, ice flow models that rely on initialization-by-inversion may already implicitly incorporate some effect of fabric by inverting for *E* (or equivalently, *A*) using remote sensing data (e.g. MacAyeal, 1993; Ranganathan and others, 2020). However, inferring *E* this way combines multiple rheological effects into a single factor in the sense that $E = E_{\text{fabric}} E_{\text{damage}} E_{\text{dust}} \dots$, making it difficult to separate individual contributions (see also Minchew and others, 2018). This work suggests new ways in which individual contributions might be disentangled. For example, estimates of *E*, inferred from inversions using an ice flow model, could be divided by $E_{\text{fabric}}$ as calculated above for Ross and PIG, allowing one to isolate how other mechanisms than fabric cause spatial variations in viscosity. Alternatively, $E_{\text{fabric}}$ might be used as a prior to inform Bayesian inversions of *E*.

### 5.3. Relevance for ice streams

This initial application of our anisotropic SSA model focuses on ice shelves due to their simplified stress state and boundary conditions, but SSA models are commonly applied to ice streams as well. Applying our method to grounded ice requires more caution than for ice shelves. Asserting a depth-constant fabric is problematic near ice-stream sticky spots or topographical bumps where fabrics can develop that facilitate flows in potential in violation of the SSA assumptions (Rathmann and Lilien, 2021; Zhang and others, 2024). For example, if the fabric develops in a way that is favorable for vertical shear, the assumption of plug flow may be less valid, similar to the effect of localized shear heating above topographical bumps (Liu and others, 2024). On the other hand, the undisturbed internal-layer structure of the North-East Greenland ice stream (NEGIS) (Jansen and others, 2024), together with surface velocity observations, indicates that the ice stream deforms by plug flow (Gerber and others, 2023). Although NEGIS is unique in being relatively flat-bedded (at least upstream), this suggests that anisotropic SSA modeling might be relevant in some parts of some ice streams. On that note, the near-bed fabrics of NEGIS (Stoll and others, 2024) are affected by DDRX which would support faster vertical shearing compared to shallower fabrics if vertical shear stresses are important, similar to the above-mentioned situation over bed bumps or sticky spots. But the fabric-induced softness for vertical shear is irrelevant if the stress regime is predominantly longitudinal (assuming fabric compatibility is good), which the observed plug flow would seem to indicate.

Because it is difficult to judge where applications of our method to ice streams will fail without comparing to a higher-order or full-Stokes model, it may be prudent simply to use a more complex flow model for grounded ice. That is not to say our results are irrelevant for such regions; considerable computational expense could be avoided by calculating *E* using EIE or CAFFE, rather than adopting an anisotropic rheology, regardless of stress balance approximation.

### 5.4. Marine ice accretion

The ice shelf models considered here rely on the assumption of a bulk rheology relevant to meteoric ice. In real ice shelves, there

are places where this assumption is less valid, such as when significant marine ice layers are present. Marine ice forms because of the accretion of platelet crystals on the ice shelf base. The microstructural, chemical, and thermal characteristics of marine ice differ significantly from meteoric ice. In particular, it is rheologically weaker than standard ice when in the tertiary creep phase (Craw, 2023), which may be related to the presence of liquid on the grain boundaries of marine ice crystals and hence the activation of grain boundary sliding (Barnes and others, 1971; Dash and others, 2006) that could facilitate lower viscosity and lead to faster ice flow compared to meteoric ice.

Marine ice is unlikely to be widespread in Antarctica, but has been inferred—either directly from borehole measurements or indirectly from satellite observations or ocean modeling—on the Amery Ice Shelf (Fricker and others, 2001; Craven and others, 2005; Treverrow and others, 2010; Galton-Fenzi and others, 2012), the Filchner-Ronne Ice Shelf (Moore and others, 1994), and the Nansen Ice Shelf (Khazendar and Jenkins, 2003; Dierckx and Tison, 2013; Dierckx and others, 2014). On these ice shelves, a bulk rheology that matches that of meteoric ice may not suffice, and changes in the distribution of marine ice over time may also impact the relevance of the steady-state fabric assumption. It is also worth noting that modeling of an idealized ice shelf using ESTAR (Craw and others, 2023) indicates that while fabric differences between marine and meteoric ice may be significant, variations in the vertical temperature distribution through an ice shelf had an order of magnitude greater impact on flow dynamics.

### 5.5. Fabric model shortcomings

Simulating ice fabric evolution, whether in laboratory conditions or in glaciers and ice sheets, requires accounting for both lattice rotation due to dislocation glide (intracrystalline slip) and accommodation processes such as dynamic recrystallization (DRX). Both processes are related to the strong viscoplastic anisotropy of ice. First, because the predominance of dislocation glide in the basal plane results in strongly anisotropic fabrics. Second, because DRX mechanisms are driven by the stored strain energy, which is related to deformation incompatibilities between polycrystal grains. Any modeling approach that does not provide an accurate representation of these processes will result in biased predictions of fabric evolution that somehow need to be compensated for.

Castelnau and others (1996) demonstrated this in detail, showing that the Sachs (homogeneous stress) and Taylor (homogeneous strain) approximations are bound estimates, unable to sufficiently simulate the mechanical response of ice polycrystals and hence fabric evolution. By neglecting grain interactions, the two bounds lead to an over- and underestimation of the mechanical anisotropy, respectively, and therefore also of lattice rotation due to intracrystalline slip. In the coupled model approach presented here, fabric evolution relies on the Taylor approximation, but this is (to some extent) compensated for by the $\alpha$ parameter, which controls how much the Sachs approximation contributes to the fabric-induced mechanical anisotropy (and hence large-scale flow, which drives fabric evolution).

The formulation of Placidi and others (2010) used here relates the occurrence of DDRX to the deformability of an individual grain, defined as the square of the basal-plane resolved shear stress that depends on the basal Schmid factor. This approach implicitly assumes that each grain is subjected to the same macro (bulk) stress, regardless of the stress redistribution resulting from intergranular interactions. The impact of these interactions on the
redistribution of stress and strain has been investigated experimentally and using full-field modeling (Grennerat and others, 2012; Piazolo and others, 2015). Results suggest that there is no correlation between the Schmid factor (as used here) and the amount of strain experienced by a grain, but also that local stresses can significantly deviate from the applied (bulk) stress. This stress heterogeneity probably explains why DDRX fabrics observed in the laboratory (e.g. Bouchez and Duval, 1982; Jacka, 1984; Montagnat and others, 2015; Journaux and others, 2019; Qi and others, 2019) are close to, but not exactly, the most favorable fabric for an applied stress. To some extent, this effect can be viewed as a dispersion of grain orientations that could be taken into account by including CDRX as a diffusive process in orientation space (Richards and others, 2021).

Most of the simulated fabrics presented here are in good agreement with recent full-field finite-element model predictions made by Chauve and others 2024, where the local (finite element scale) resolved shear stress is used to calculate an attractor toward which c-axis rotate due to DDRX. Using this approach, fabrics observed under various laboratory conditions can be reproduced by appropriately adjusting one single DDRX kinetic parameter. In particular, the double-maximum transitory fabrics that develop under simple shear (Hudleston, 1977; Bouchez and Duval, 1982; Qi and others, 2019) are well reproduced. Despite the increased physical realism of such full-field models compared to our fabric model that treats grains as interactionless and assumes crude stress/strain homogenizations, our approach has also been shown able to reproduce complicated transitory fabrics (Richards and others, 2021; Wang and others, 2025b). But how exactly the rate of DDRX is to be parameterized as a function of temperature and stress/strain-rate magnitude is not settled; recent work suggests that lab-calibrated rate functions are not necessarily applicable to large-scale ice flow modeling in East Antarctica (Lilien and others, 2023).

## 6. Conclusions

We investigated the rheological control (ease of deformation) exerted by well-developed crystal orientation fabrics on the flow of the Ross and Pine Island Glacier (PIG) ice shelves, Antarctica, by calculating maps of fabric-induced flow enhancement factors using the CAFFE model (Placidi and others, 2010). To do so, we estimated steady-state fabric fields over Ross and PIG by solving a high-dimensional boundary value problem, assuming depth-independent horizontal velocities (plug flow) and that the ice accumulation rate is small compared to the local ice thickness. In effect, the problem is closed by prescribing satellite-derived surface velocities and ice temperatures. Since the latter is not known in sufficient spatial detail, we considered the two end-member cases in which the Ross and PIG regions are either very cold or warm, so that dynamic recrystallization is either negligible or strong, respectively. The two solutions show that significant ice-shelf hardening or softening can occur depending on whether dynamic recrystallization is important or not. This emphasizes the ice-dynamical relevance of needing to better constrain the strength of fabric processes, which in turn calls on widespread fabric and temperature measurements from the field that are currently missing or too sparsely available for model validation.

To increase our confidence in the relevance of our calculated enhancement factor maps, we tested how well CAFFE, when combined with the Glen isotropic rheology, can reproduce the velocity field of an idealized ice-shelf box model based on an orthotropic

plug flow rheology (i.e., full fabric-flow coupling). We find that replacing the tensorial viscosity structure by a scalar flow enhancement factor leads to velocity errors on the order of ten percent, but if DDRX is strong, the error might be a bit larger. Given that neglecting the effect of ice fabric altogether leads to modeled velocity errors between thirty to sixty percent, this is much of an improvement. We therefore argue that in some instances it might be reasonable to simply use CAFFE in large-scale flow models to account for the rheological control exerted by fabric development, which can easily be implemented when combined with the spectral fabric model presented here.

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

**Data availability.** The FEniCS model code and scripts for generating all plots are available at https://github.com/nicholasmr/specfab.

**Acknowledgements.** We thank Editor Alan Rempel, Sheng Fan, and an anonymous reviewer for providing thoughtful comments that helped improve our work.

**Funding statement.** NMR was supported by the Independent Research Fund Denmark (DFF) grant no. 2032-00364B and the Novo Nordisk Foundation Challenge grant no. NNF23OC0081251.

DHR was supported by the Australian Research Council Special Research Initiative, Australian Centre for Excellence in Antarctic Science (SR200100008), and the UK Natural Environment Research Council (grant no. NE/X014991/1)

FSM was supported by an Australian Research Council Discovery Early Career Research Award (DE210101433), the Special Research Initiative Securing Antarctica's Environmental Future (SR200100005), and a French-Australian Cooperative Research Grant.

MM was supported by the European Research Council (ERC) grant no. 882450 under the European Union's Horizon 2020 Research and Innovation program (ERC RhEoVOLUTION).

**Competing interests.** The authors have no competing interests.

## References

Adusumilli S, Fricker HA, Siegfried MR, Padman L, Paolo FS and Ligtenberg SRM (2018) Variable basal melt rates of Antarctic Peninsula ice shelves, 1994–2016. *Geophysical Research Letters* 45(9), 4086–4095. doi:10.1002/2017GL076652

Alnæs MS, Logg A, Ølgaard KB, Rognes ME and Wells GN (2014) Unified Form Language: A domain-specific language for weak formulations of partial differential equations. *ACM Transactions on Mathematical Software* 40(2), 1–37. doi:10.1145/2566630

Andreasen JR, Hogg AE and Selley HL (2023) Change in Antarctic ice shelf area from 2009 to 2019. *The Cryosphere* 17(5), 2059–2072. doi:10.5194/tc-17-2059-2023

Azuma N (1994) A flow law for anisotropic ice and its application to ice sheets. *Earth and Planetary Science Letters* 128(3), 601–614. doi:10.1016/0012-821X(94)90173-2

Azuma N (1995) A flow law for anisotropic polycrystalline ice under uniaxial compressive deformation. *Cold Regions Science and Technology* 23(2), 137–147. doi:10.1016/0165-232X(94)00011-L

Barnes P, Tabor D and Walker JCF (1971) The friction and creep of polycrystalline ice. *Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences* 324(1557), 127–155. doi:10.1098/rspa.1971.0132

Behn MD, Goldsby DL and Hirth G (2021) The role of grain size evolution in the rheology of ice: implications for reconciling laboratory creep data and the Glen flow law. *The Cryosphere* 15(9), 4589–4605. doi:10.5194/tc-15-4589-2021

Bons PD and 6 others (2018) Greenland ice sheet: higher nonlinearity of ice flow significantly reduces estimated basal motion. *Geophysical Research Letters* 45(13), 6542–6548. doi:10.1029/2018GL078356
Borstad C, McGrath D and Pope A (2017) Fracture propagation and stability of ice shelves governed by ice shelf heterogeneity. *Geophysical Research Letters* **44**(9), 4186–4194. doi:10.1002/2017GL072648

Bouchez J and Duval P (1982) The fabric of polycrystalline ice deformed in simple shear: experiments in torsion, natural deformation and geometrical interpretation. *Texture, Stress, and Microstructure* **5**(3), 171–190. doi:10.1155/TSM.5.171

Castelnau O and 7 others (1998) Anisotropic behavior of GRIP ices and flow in Central Greenland. *Earth and Planetary Science Letters* **154**(1), 307–322. doi:10.1016/S0012-821X(97)00193-3

Castelnau O, Duval P, Lebensohn RA and Canova GR (1996) Viscoplastic modeling of texture development in polycrystalline ice with a self-consistent approach: Comparison with bound estimates. *Journal of Geophysical Research: Solid Earth* **101**(B6), 13851–13868. doi:10.1029/96JB00412

Catania GA, Scambos TA, Conway H and Raymond CF (2006) Sequential stagnation of Kamb Ice Stream, West Antarctica. *Geophysical Research Letters* **33**(14), L14502. doi:10.1029/2006GL026430

Chauve T, Montagnat M, Dansereau V, Saramito P, Fourteau K and Tommasi A (2024) A physically-based formulation for texture evolution during dynamic recrystallization. A case study for ice. *Article soumis* **352**, 99–134. doi:10.5194/egusphere-egu24-7333

Conway H, Catania G, Raymond C, Gades A, Scambos T and Engelhardt H (2002) Switch of flow direction in an Antarctic ice stream. *Nature* **419**(6906), 465–467. doi:10.1038/nature01081

Craven M and 7 others (2005) Borehole imagery of meteoric and marine ice layers in the Amery Ice Shelf, East Antarctica. *Journal of Glaciology* **51**(172), 75–84. doi:10.3189/172756505781829511

Craw L (2023) *University of Tasmania*, The influence of marine ice on ice shelf dynamics and stability. *PhD Thesis*, doi:10.25959/25209671.v1

Craw L, McCormack FS, Cook S, Roberts J and Treverrow A (2023) Modelling the influence of marine ice on the dynamics of an idealised ice shelf. *Journal of Glaciology* **69**(274), 342–352. doi:10.1017/jog.2022.66

Cuffey KM and Paterson WSB (2010) *The Physics Of Glaciers*. Amsterdam: Academic Press

Dahl-Jensen D and Gundestrup N (1987) Constitutive properties of ice at Dye 3, Greenland. *International Association of Hydrological Sciences Publication* **170**, 31–43

Dash JG, Rempel AW and Wettlaufer JS (2006) The physics of premelted ice and its geophysical consequences. *Rev. Mod. Phys.* **78**, 695–741. doi:10.1103/RevModPhys.78.695

De La Chapelle S, Castelnau O, Lipenkov V and Duval P (1998) Dynamic recrystallization and texture development in ice as revealed by the study of deep ice cores in Antarctica and Greenland. *Journal of Geophysical Research: Solid Earth* **103**(B3), 5091–5105. doi:10.1029/97JB02621

De La Chapelle S, Milsch H, Castelnau O and Duval P (1999) Compressive creep of ice containing a liquid intergranular phase: Rate-controlling processes in the dislocation creep regime. *Geophysical Research Letters* **26**(2), 251–254. doi:10.1029/1998GL900289

Dierckx M, Peternell M, Schroeder C and Tison JL (2014) Influence of pre-existing microstructure on mechanical properties of marine ice during compression experiments. *Journal of Glaciology* **60**(221), 576–586. doi:10.3189/2014JoG13j154

Dierckx M and Tison JL (2013) Marine ice deformation experiments: an empirical validation of creep parameters. *Geophysical Research Letters* **40**(1), 134–138. doi:10.1029/2012GL054197

Diez A and 8 others (2016) Ice shelf structure derived from dispersion curve analysis of ambient seismic noise, Ross Ice Shelf, Antarctica. *Geophysical Journal International* **205**(2), 785–795. doi:10.1093/gji/ggw036

Durand G and 8 others (2007) Change in ice rheology during climate variations – implications for ice flow modelling and dating of the EPICA Dome C core. *Climate of the Past* **3**(1), 155–167. doi:10.5194/cp-3-155-2007

Duval P, Ashby MF and Anderman I (1983) Rate-controlling processes in the creep of polycrystalline ice. *The Journal of Physical Chemistry* **87**(21), 4066–4074. doi:10.1021/j100244a014

Echelmeyer KA, Harrison WD, Larsen C and Mitchell JE (1994) The role of the margins in the dynamics of an active ice stream. *Journal of Glaciology* **40**(136), 527–538. doi:10.3189/S0022143000012417

Fan S and 6 others (2021) Crystallographic Preferred Orientation (CPO) development governs strain weakening in ice: insights from high-temperature deformation experiments. *Journal of Geophysical Research: Solid Earth* **126**(12), e2021JB023173. https://doi.org/10.1029/2021JB023173

Fan S and Prior DJ (2023) Cool ice with hot properties. *Nature Geoscience* **16**(12), 1073–1073. doi:10.1038/s41561–01330-z

Fan S, Wang T, Prior DJ, Breithaupt T, Hager TF and Wallis D (2025) Flow laws for ice constrained by 70 years of laboratory experiments. *Nature Geoscience* **18**, 296–304. doi:10.1038/s41561–01661-z

Faria S (2001) Mixtures with continuous diversity: general theory and application to polymer solutions. *Continuum Mechanics and Thermodynamics* **13**(2), 91–120. doi:10.1007/s001610100043

Faria SH, Weikusat I and Azuma N (2014a) The microstructure of polar ice. Part I: Highlights from ice core research. *Journal of Structural Geology* **61**, 2–20. doi:10.1016/j.jsg.2013.09.010

Faria SH, Weikusat I and Azuma N (2014b) The microstructure of polar ice. Part II: State of the art. *Journal of Structural Geology* **61**, 21–49. doi:10.1016/j.jsg.2013.11.003

Favier L and 8 others (2014) Retreat of Pine Island Glacier controlled by marine ice-sheet instability. *Nature Climate Change* **4**(2), 117–121. doi:10.1038/nclimate2094

Fowler JR and Iverson NR (2023) The relationship between the permeability and liquid water content of polycrystalline temperate ice. *Journal of Glaciology* **69**(278), 2154–2162. doi:10.1017/jog.2023.91

Fricker HA, Popov S, Allison I and Young N (2001) Distribution of marine ice beneath the Amery Ice Shelf. *Geophysical Research Letters* **28**(11), 2241–2244. doi:10.1029/2000GL012461

Galton-Fenzi BK, Hunter JR, Coleman R, Marsland SJ and Warner RC (2012) Modeling the basal melting and marine ice accretion of the Amery Ice Shelf. *Journal of Geophysical Research: Oceans* **117**(C9), C09031. doi:10.1029/2012JC008214

Gerber TA and 10 others (2023) Crystal orientation fabric anisotropy causes directional hardening of the Northeast Greenland Ice Stream. *Nature Communications* **14**(1), 2653. doi:10.1038/s41467-023-38139-8

Gerli C, Rosier S, Gudmundsson GH and Sun S (2024) Weak relationship between remotely detected crevasses and inferred ice rheological parameters on Antarctic ice shelves. *The Cryosphere* **18**(6), 2677–2689. doi:10.5194/tc-18-2677-2024

Gillet-Chaulet F, Gagliardini O, Meyssonnier J, Montagnat M and Castelnau O (2005) A user-friendly anisotropic flow law for ice-sheet modeling. *Journal of Glaciology* **51**(172), 3–14. doi:10.3189/172756505781829584

Gillet-Chaulet F, Gagliardini O, Meyssonnier J, Zwinger T and Ruokolainen J (2006) Flow-induced anisotropy in polar ice and related ice-sheet flow modelling. *Journal of Non-Newtonian Fluid Mechanics* **134**(1), 33–43. doi:10.1016/j.jnnfm.2005.11.005

Goldsby D and Kohlstedt D (1997) Grain boundary sliding in fine-grained Ice I. *Scripta Materialia* **37**(9), 1399–1406. doi:10.1016/S1359-6462(97)00246-7

Goldsby DL and Kohlstedt DL (2001) Superplastic deformation of ice: Experimental observations. *Journal of Geophysical Research: Solid Earth* **106**(B6), 11017–11030. doi:10.1029/2000JB900336

Graham FS, Morlighem M, Warner RC and Treverrow A (2018) Implementing an empirical scalar constitutive relation for ice with flow-induced polycrystalline anisotropy in large-scale ice sheet models. *The Cryosphere* **12**(3), 1047–1067. doi:10.5194/tc-12-1047-2018

Grennerat F and 6 others (2012) Experimental characterization of the intra-granular strain field in columnar ice during transient creep. *Acta Materialia* **60**(8), 3655–3666. doi:10.1016/j.actamat.2012.03.025

Hamann I, Weikusat C, Azuma N and Kipfstuhl S (2007) Evolution of ice crystal microstructure during creep experiments. *Journal of Glaciology* **53**(182), 479–489. doi:10.3189/002214307783258341

Hill EA, Gudmundsson GH, Carr JR, Stokes CR and King HM (2021) Twenty-first century response of Petermann Glacier, northwest Greenland to ice shelf loss. *Journal of Glaciology* **67**(261), 147–157. doi:10.1017/jog.2020.97

Hudleston PJ (1977) *Progressive Deformation and Development of Fabric Across Zones of Shear in Glacial ice*, Berlin Heidelberg, Berlin, Heidelberg: Springer. pp. 121–150. doi:10.1007/978-3-642-86574-9_7
Hulbe CL and Fahnestock MA (2004) West Antarctic ice-stream discharge variability: mechanism, controls and pattern of grounding-line retreat. *Journal of Glaciology* **50**(171), 471–484. doi:10.3189/172756504781829738

Huth A, Duddu R and Smith B (2021) A generalized interpolation material point method for shallow ice shelves. 2: Anisotropic nonlocal damage mechanics and rotor propagation. *Journal of Advances in Modeling Earth Systems* **13**(8), e2020MS002292. doi:10.1029/2020MS002292

Ismail WB and Mainprice D (1998) An olivine fabric database: an overview of upper mantle fabrics and seismic anisotropy. *Tectonophysics* **296**(1), 145–157. doi:10.1016/S0040–(98)00141–

Jacka T (1984) Laboratory studies on relationships between ice crystal size and flow rate. *Cold Regions Science and Technology* **10**(1), 31–42. doi:10.1016/0165-232X(84)90031-4

Jacka T and Budd W (1989) Isotropic and anisotropic flow relations for ice dynamics. *Annals of Glaciology* **12**, 81–84. doi:10.3189/S026030550006996

Jacka TH and L Jun (2000). Flow rates and crystal orientation fabrics in compression of polycrystalline ice at low temperatures and stresses. In *Physics of Ice Core Records*, Hokkaido University Press, pp.83–102.

Jackson M and Kamb B (1997) The marginal shear stress of Ice Stream B, West Antarctica. *Journal of Glaciology* **43**(145), 415–426. doi:10.3189/S0022143000035000

Jansen D and 9 others (2024) Shear margins in upper half of Northeast Greenland Ice Stream were established two millennia ago. *Nature Communications* **15**(1), 1193. doi:10.1038/s41467-024-45021-8

Joughin I, Shapero D, Smith B, Dutrieux P and Barham M (2021) Ice-shelf retreat drives recent Pine Island Glacier speedup. *Science Advances* **7**(24), eabg3080. doi:10.1126/sciadv.abg3080

Joughin I, Tulaczyk S, Bindschadler R and Price SF (2002) Changes in west Antarctic ice stream velocities: Observation and analysis. *Journal of Geophysical Research: Solid Earth* **107**(B11), EPM 3–1–EPM 3–22. doi:10.1029/2001JB001029

Journaux B and 6 others (2019) Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear. *The Cryosphere* **13**(5), 1495–1511. doi:10.5194/tc-13-1495-2019

Khazendar A and Jenkins A (2003) A model of marine ice formation within Antarctic ice shelf rifts. *Journal of Geophysical Research: Oceans* **108**(C7), 3235. doi:10.1029/2002JC001673

Kulessa B and 10 others (2019) Seawater softening of suture zones inhibits fracture propagation in Antarctic ice shelves. *Nature Communications* **10**(1), 5491. doi:10.1038/s41467–13539-x

LeDoux CM, Hulbe CL, Forbes MP, Scambos TA and Alley K (2017) Structural provinces of the Ross Ice Shelf, Antarctica. *Annals of Glaciology* **58**(75pt1), 88–98. doi:10.1017/aog.2017.24

Leung J, Hudson TS, Kendall JM and Barcheck G (2025) Evidence that seismic anisotropy captures upstream palaeo-ice fabric: Implications on present-day deformation at Whillans Ice Stream, Antarctica. *Journal of Glaciology* **71**, e39. doi:10.1017/jog.2025.19

Lilien DA and 6 others (2023) Simulating higher-order fabric structure in a coupled, anisotropic ice-flow model: application to Dome C. *Journal of Glaciology* **69**, 2007–2026. doi:10.1017/jog.2023.78

Lilien DA, Rathmann NM, Hvidberg CS and Dahl-Jensen D (2021) Modeling ice-crystal fabric as a proxy for ice-stream stability. *Journal of Geophysical Research: Earth Surface* **126**(9), e2021JF006306. doi:10.1029/2021JF006306

Liu EW, Räss L, Herman F, Podladchikov Y and Suckale J (2024) Spontaneous formation of an internal shear band in ice flowing over topographically variable bedrock. *Journal of Geophysical Research: Earth Surface* **129**(4), e2022JF007040. doi:10.1029/2022JF007040

Logg A and Mardal KA (2012) *Automated Solution Of Differential Equations By The Finite Element Method*. Berlin, Heidelberg: Springer. doi:10.1007/978-3-642-23099-8

Ma Y, Gagliardini O, Ritz C, Gillet-Chaulet F, Durand G and Montagnat M (2010) Enhancement factors for grounded ice and ice shelves inferred from an anisotropic ice-flow model. *Journal of Glaciology* **56**(199), 805–812. doi:10.3189/002214310794457209

MacAyeal DR (1989) Large-scale ice flow over a viscous basal sediment: Theory and application to Ice Stream B, Antarctica. *Journal of Geophysical Research: Solid Earth* **94**(B4), 4071–4087. doi:10.1029/JB094iB04p04071

MacAyeal DR (1993) A tutorial on the use of control methods in ice-sheet mod-
eling. *Journal of Glaciology* **39**(131), 91–98. doi:10.3189/S0022143000015744

Mangeney A, Califano F and Castelnau O (1996) Isothermal flow of an anisotropic ice sheet in the vicinity of an ice divide. *Journal of Geophysical Research: Solid Earth* **101**(B12), 28189–28204. doi:10.1029/96JB01924

Martín C and Gudmundsson GH (2012) Effects of nonlinear rheology, tem-
perature and anisotropy on the relationship between age and depth at ice
divides. *The Cryosphere* **6**(5), 1221–1229. doi:10.5194/tc-6-1221-2012

Martín C, Gudmundsson GH, Pritchard HD and Gagliardini O (2009) On the effects of anisotropic rheology on ice flow, internal structure, and the age-depth relationship at ice divides. *Journal of Geophysical Research: Earth Surface* **114**(F4), 2008JF001204. doi:10.1029/2008JF001204

McCormack FS, Warner RC, Seroussi H, Dow CF, Roberts JL and Treverrow A (2022) Modeling the deformation regime of Thwaites Glacier, West Antarctica, using a simple flow relation for ice anisotropy (ESTAR). *Journal of Geophysical Research: Earth Surface* **127**(3), e2021JF006332. doi:10.1029/2021JF006332

Miles K and 6 others (2025) Influence of the grounding zone on the inter-
nal structure of ice shelves. *Nature Communications* **16**(1), 1–9. doi:10.1038/
s41467-025-58973-2

Minchew BM, Meyer CR, Robel AA, Gudmundsson GH and Simons M (2018) Processes controlling the downstream evolution of ice rheology in glacier shear margins: case study on Rutford Ice Stream, West Antarctica. *Journal of Glaciology* **64**(246), 583–594. doi:10.1017/jog.2018.47

Montagnat M and 9 others (2014) Fabric along the NEEM ice core, Greenland, and its comparison with GRIP and NGRIP ice cores. *The Cryosphere* **8**(4), 1129–1138. doi:10.5194/tc-8-1129-2014

Montagnat M, Chauve T, Barou F, Tommasi A, Beausir B and Fressengeas C (2015) Analysis of dynamic recrystallization of ice from EBSD orientation mapping. *Frontiers in Earth Science* **3**(81). doi:10.3389/feart.2015.00081

Montagnat M and Duval P (2000) Rate controlling processes in the creep of polar ice, influence of grain boundary migration associated with recrystallization. *Earth and Planetary Science Letters* **183**(1), 179–186. doi:10.1016/S0012-821X(00)00262-4

Monz ME and 7 others (2021) Full crystallographic orientation (c and a axes) of warm, coarse-grained ice in a shear-dominated setting: a case study, Storglaciären, Sweden. *The Cryosphere* **15**(1), 303–324. doi:10.5194/tc-15-303-2021

Moore JC, Reid AP and Kipfstuhl J (1994) Microstructure and electrical prop-
erties of marine ice and its relationship to meteoric ice and sea ice. *Journal
of Geophysical Research: Oceans* **99**(C3), 5171–5180. doi:10.1029/93JC0
2832

Morland L (1987) *Unconfined ice-shelf flow*. In *Dynamics of the West Antarctic Ice Sheet*, Springer, pp.99–116. doi:10.1007/978-94-009-3745-1

Naumenko K and Altenbach H (2007) Modeling Of Creep For Structural Analysis, Berlin, Heidelberg: Springer Science & Business Media. doi:10.1007/978-3-540-70839-1

Otosaka IN and 67 others (2023) Mass balance of the Greenland and Antarctic ice sheets from 1992 to 2020. *Earth System Science Data* **15**(4), 1597–1616. doi:10.5194/essd-15-1597-2023

Pettit EC, Thorsteinsson T, Jacobson HP and Waddington ED (2007) The role of creep in flow near an ice divide. *Journal of Glaciology* **53**(181), 277–288. doi:10.3189/172756507782202766

Piazolo S, Montagnat M, Grennerat F, Moulinec H and Wheeler J (2015) Effect of local stress heterogeneities on dislocation fields: Examples from transient creep in polycrystalline ice. *Acta Materialia* **90**, 303–309. doi:10.1016/j.actamat.2015.02.046

Pimienta P, P Duval and VY Lipenkov (1987) *Mechanical behaviour of anisotropic polar ice*. In *The Physical Basis of Ice Sheet Modelling*, IAHS Publication, (170), IAHS Press, Wallingford, UK, pp.57–66.

Placidi L, Greve R, Seddik H and Faria SH (2010) Continuum-mechanical, anisotropic flow model for polar ice masses, based on an anisotropic flow enhancement factor. *Continuum Mechanics and Thermodynamics* **22**(3), 221–237. doi:10.1007/s00161-009-0126-0

Qi C and 8 others (2019) Crystallographic preferred orientations of ice deformed in direct-shear experiments at low temperatures. *The Cryosphere* **13**(1), 351–371. doi:10.5194/tc-13-351-2019
Ranganathan M and Minchew B (2024) A modified viscous flow law for natural glacier ice: Scaling from laboratories to ice sheets. *Proceedings of the National Academy of Sciences* **121**(23), e2309788121. doi:10.1073/pnas. 2309788121

Ranganathan M, Minchew B, Meyer CR and Gudmundsson GH (2020) A new approach to inferring basal drag and ice rheology in ice streams, with applications to West Antarctic Ice Streams. *Journal of Glaciology* **67**, 229–242. doi:10.1017/jog.2020.95

Ranganathan M, Minchew B, Meyer CR and Peč M (2021) Recrystallization of ice enhances the creep and vulnerability to fracture of ice shelves. *Earth and Planetary Science Letters* **576**, 117219. doi:10.1016/j.jpgl.2021.11. 7219.

Rathmann NM and 7 others (2017) Highly temporally resolved response to seasonal surface melt of the Zachariae and 79N outlet glaciers in northeast Greenland. *Geophysical Research Letters* **44**(19), 9805–9814. doi:10.1002/ 2017GL074368

Rathmann NM and 9 others (2022) Elastic wave propagation in anisotropic polycrystals: inferring physical properties of glacier ice. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* **478**(2268), 20220574. doi:10.1098/rspa.2022.0574

Rathmann NM, Hvidberg CS, Grinsted A, Lilien DA and Dahl-Jensen D (2021) Effect of an orientation-dependent non-linear grain fluidity on bulk directional enhancement factors. *Journal of Glaciology* **67**, 569–575. doi:10.1017/jog.2020.117

Rathmann NM and Lilien DA (2021) Inferred basal friction and mass flux affected by crystal-orientation fabrics. *Journal of Glaciology* **68**, 236–252. doi:10.1017/jog.2021.88

Rathmann NM and Lilien DA (2022) On the nonlinear viscosity of the orthotropic bulk rheology. *Journal of Glaciology* **68**, 1243–1248. doi:10.1017/ jog.2022.33

Reese R, Gudmundsson GH, Levermann A and Winkelmann R (2018) The far reach of ice-shelf thinning in Antarctica. *Nature Climate Change* **8**(1), 53–57. doi:10.1038/s41558–0020-x

Retzlaff R and Bentley CR (1993) Timing of stagnation of Ice Stream C, West Antarctica, from short-pulse radar studies of buried surface crevasses. *Journal of Glaciology* **39**(133), 553–561. doi:10.3189/S0022143000016440

Richards DH, Pegler SS, Piazolo S and Harlen OG (2021) The evolution of ice fabrics: A continuum modelling approach validated against laboratory experiments. *Earth and Planetary Science Letters* **556**, 116718. doi:10.1016/j. jpgl.2020.116718.

Richards DH, Pegler SS, Piazolo S, Stoll N and Weikusat I (2023) Bridging the gap between experimental and natural fabrics: modeling ice stream fabric evolution and its comparison with ice-core data. *Journal of Geophysical Research: Solid Earth* **128**(11), e2023JB027245. doi:10.1029/2023JB02 7245

Rignot E, Casassa G, Gogineni P, Krabill W, Rivera A and Thomas R (2004) Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf. *Geophysical Research Letters* **31**(18), 2004GL020697. doi:10.1029/2004GL020697

Rignot E, Mouginot J and Scheuchl B (2017) MEaSUREs InSAR-Based Antarctica Ice Velocity Map, Version 2, Boulder, Colorado, USA. doi:10.5067/ D7GK8F5J8M8R

Samyn D, Svensson A and Fitzsimons SJ (2008) Dynamic implications of discontinuous recrystallization in cold basal ice: Taylor Glacier, Antarctica.

*Journal of Geophysical Research: Earth Surface* **113**(F3), 2006JF000600. doi:10.1029/2006JF000600

Seroussi H and 48 others (2023) Insights into the vulnerability of Antarctic glaciers from the ISMIP6 ice sheet model ensemble and associated uncertainty. *The Cryosphere* **17**(12), 5197–5217. doi:10.5194/tc-17-5197-2023

Shepherd A and 9 others (2019) Trends in Antarctic ice sheet elevation and mass. *Geophysical Research Letters* **46**(14), 8174–8183. doi:10.1029/ 2019GL082182

Shoji H and CC Langway (1985) *Mechanical Properties of Fresh Ice Core from Dye 3. Greenland*, American Geophysical Union (AGU), pp.39–48. doi:10.1029/GM033p0039.

Staroszczyk R and Gagliardini O (1999) Two orthotropic models for strain-induced anisotropy of polar ice. *Journal of Glaciology* **45**(151), 485–494. doi:10.3189/S0022143000001349

Stoll N and 15 others (2024) *EastGRIP ice core reveals the exceptional evolution of crystallographic preferred orientation throughout the Northeast Greenland Ice Stream*. *EGUsphere*, **2024**, 1–34. doi:10.5194/egusphere-2024-2653

Sun S, Cornford SL, Moore JC, Gladstone R and Zhao L (2017) Ice shelf fracture parameterization in an ice sheet model. *The Cryosphere* **11**(6), 2543–2554. doi:10.5194/tc-11-2543-2017

Sun S and Gudmundsson GH (2023) The speedup of Pine Island Ice Shelf between 2017 and 2020: reevaluating the importance of ice damage. *Journal of Glaciology* **69**, 1983–1991. doi:10.1017/jog.2023.76

Svendsen B and Hutter K (1996) A continuum approach for modelling induced anisotropy in glaciers and ice sheets. *Annals of Glaciology* **23**, 262–269. doi:10.3189/S0260305500013525

Thomas RE and 11 others (2021) Microstructure and crystallographic preferred orientations of an azimuthally oriented ice core from a lateral shear margin: Priestley Glacier, Antarctica. *Frontiers in Earth Science* **9**, 702213. doi:10.3389/feart.2021.702213

Treverrow A, Warner RC, Budd WF and Craven M (2010) Meteoric and marine ice crystal orientation fabrics from the Amery Ice Shelf, East Antarctica. *Journal of Glaciology* **56**(199), 877–890. doi:10.3189/ 002214310794457353

Treverrow A, Warner RC, Budd WF, Jacka T and Roberts JL (2015) Modelled stress distributions at the Dome Summit South borehole, Law Dome, East Antarctica: a comparison of anisotropic ice flow relations. *Journal of Glaciology* **61**(229), 987–1004. doi:10.3189/2015Jg0G14J198

Wang L and 10 others (2025a) Estimating marine ice thickness beneath the Amery ice shelf from Airborne radio-echo sounding. *Journal of Glaciology* **71**, e75. doi:10.1017/jog.2025.42

Wang Q, Fan S, Richards DH, Worthington R, Prior DJ and Qi C (2025b) Evolution of crystallographic preferred orientations of ice sheared to high strains by equal-channel angular pressing. *The Cryosphere* **19**(2), 827–848. doi:10.5194/tc-19-827-2025

Wang Y, Lai CY, Prior DJ and Cowen-Breen C (2025c) Deep learning the flow law of Antarctic ice shelves. *Science* **387**(6739), 1219–1224. doi:10.1126/ science.adp3300

Zeising O and 6 others (2023) Improved estimation of the bulk ice crystal fabric asymmetry from polarimetric phase co-registration. *The Cryosphere* **17**(3), 1097–1105. doi:10.5194/tc-17-1097-2023

Zhang Y and 6 others (2024) Formation mechanisms of large-scale folding in Greenland's ice sheet. *Geophysical Research Letters* **51**(16), e2024GL109492. doi:10.1029/2024GL109492