2.1. Momentum balance

The SSA momentum balance is (Morland, Reference Morland1987; MacAyeal, Reference MacAyeal1989)

(1)\begin{align} \boldsymbol{\nabla} \cdot (H \mathbf{R}) = \rho g H \boldsymbol{\nabla} S , \end{align}

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

(2)\begin{align} \mathbf{R} &= \begin{bmatrix} 2\tau_{xx}+\tau_{yy} & \tau_{xy} \\ \tau_{xy} & \tau_{xx}+2\tau_{yy} \\ \end{bmatrix}. \end{align}

The stress deviator components τ_xx, τ_yy and τ_xy follow from the power-law rheology of ice assuming dislocation creep is the dominant deformation mechanism,

(3)\begin{align} \boldsymbol{\tau} &= A^{-1/n} \dot{\epsilon}_\mathrm{e}^{(1-n)/n} \mathbf{C}, \end{align}

(4)\begin{align} \dot{\epsilon}_\mathrm{e} &= \sqrt{\frac{\mathbf{C}:\dot{\boldsymbol{\epsilon}}}{2}}, \end{align}

where A is the depth-averaged flow rate factor, $\dot{\boldsymbol{\epsilon}}$ is the strain rate tensor, $\dot{\epsilon}_\mathrm{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

(5)\begin{align} \dot{\boldsymbol{\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}, \end{align}

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

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 $\sim 100$.

2.2. Bulk rheology

In the case of the Glen isotropic rheology, the tensorial part is $\mathbf{C} = \dot{\boldsymbol{\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, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005) has three mutually perpendicular planes of reflection symmetry with normals m₁, m₂, and m₃, 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 Reference Lilien2023, the tensorial part of the orthotropic rheology is

(6)\begin{align} \mathbf{C} = \sum_{i=1}^3 \big( \eta_i \mathbf{P}_i \mathbf{P}_i + \eta_{i+3} \mathbf{P}_{i+3}\mathbf{P}_{i+3} \big): \dot{\boldsymbol{\epsilon}} , \end{align}

where $\mathbf{P}_i = \mathbf{I}/3-\mathbf{m}_i \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{\boldsymbol{\epsilon}}$ and $\mathbf{P}_{i+3} : \dot{\boldsymbol{\epsilon}}$ give the longitudinal and shear strain-rate components in the rheological symmetry frame 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 m_i. These can be conveniently written in terms of the fabric-induced directional enhancement factors E_ij (Rathmann and Lilien, Reference Rathmann and Lilien2022)

(7)\begin{align} \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 \left( 2 E_{j_i j_i}^{2/(n+1)} E_{k_i k_i}^{2/(n+1)} - E_{ii}^{4/(n+1)}\right) } ,\quad \eta_{i+3} = 2E_{j_i k_i}^{-2/(n+1)} . \end{align}

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)

(8)\begin{equation}E_{ij}=\frac{{\mathbf m}_i\cdot\dot{\boldsymbol\epsilon}(\hat{\boldsymbol\tau})\cdot{\mathbf m}_j}{{\mathbf m}_i\cdot{\dot{\boldsymbol\epsilon}}_\text{iso}(\hat{\boldsymbol\tau})\cdot{\mathbf m}_j},\end{equation}

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

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

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

We mention in passing that modeling ice flow using an orthotropic rheology is not new (Staroszczyk and Gagliardini, Reference Staroszczyk and Gagliardini1999; Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006; Ma and others, Reference Ma, Gagliardini, Ritz, Gillet-Chaulet, Durand and Montagnat2010; Lilien and others, Reference Lilien, Rathmann, Hvidberg and Dahl-Jensen2021; Reference Lilien2023) and is capable of (i) reproducing the viscous anisotropy of Dye 3 ice core samples (Rathmann and Lilien Reference Rathmann and Lilien2022; revisited in Supplementary section B), (ii) modeling the fabric profile of the EPICA Dome C ice core (Durand and others, Reference Durand2007), and (iii) explaining the formation of Raymond bumps with double peaks (Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009).

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 Reference Montagnat(2014) or Faria and others (Reference Faria, Weikusat and Azuma2014a, Reference Faria, Weikusat and Azuma2014b) for an introduction to the physics and observations of crystal processes.

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).

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 $\sim{-15}^{\circ}\mathrm{C}$ (De La Chapelle and others, Reference De La Chapelle, Castelnau, Lipenkov and Duval1998; Samyn and others, Reference Samyn, Svensson and Fitzsimons2008; Qi and others, Reference Qi2019).

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 (subdominant to lattice rotation), although it may impact the reduction of stored strain energy (Montagnat and Duval, Reference Montagnat and Duval2000). For this reason, and because some work (Lilien and others, Reference Lilien2023; Richards and others, Reference Richards, Pegler, Piazolo, Stoll and Weikusat2023) 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 Reference Castelnau, Duval, Lebensohn and Canova(1996) and Svendsen and Hutter Reference Svendsen and Hutter(1996), lattice rotation is modeled as an advection process on the surface of the unit sphere (S ²), which is supposed to capture the strain-induced rotation of crystal lattices (c axes). In this model, the c-axis velocity field on S ² depends solely on the bulk stretching $\dot{\boldsymbol{\epsilon}}$ and spin ω tensors and is given by $\dot{\mathbf{c}} = (\boldsymbol{\omega} + \boldsymbol{\omega}_{\mathrm{p}})\cdot \mathbf{c}$, where c is an arbitrary c-axis, $\boldsymbol{\omega}_{\mathrm{p}}=\iota (\mathbf{c}^2\cdot\dot{\boldsymbol{\epsilon}} - \dot{\boldsymbol{\epsilon}}\cdot\mathbf{c}^2)$ is the plastic spin, and c² is the outer product of c with itself. We take ι = 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 ι could take another value to mimic the activity of non-basal slip systems (Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006; Richards and others, Reference Richards, Pegler, Piazolo and Harlen2021), 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.

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 initially-isotropic fabric (d–f).

2.3.2. DDRX

Following Placidi and others Reference Placidi, Greve, Seddik and Faria(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{\boldsymbol{\epsilon}}:\dot{\boldsymbol{\epsilon}}/2} A_{\Gamma}\exp(-Q_{\Gamma}/RT)$ accounts for the preferential (Arrhenius) activation at warm temperatures and the effect of strain-rate magnitude (Richards and others, Reference Richards, Pegler, Piazolo and Harlen2021; Lilien and others, Reference Lilien2023). 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}\mathrm{J\,mol}^{-1}$ according to a recent calibration by Lilien and others Reference Lilien2023. The deformability D is the square of the basal-plane resolved shear stress, $(\boldsymbol{\tau}\cdot\boldsymbol{\tau}):\mathbf{c}^2 - \boldsymbol{\tau}:\mathbf{c}^4:\boldsymbol{\tau}$, relative to the average value if the fabric were isotropic

(10)\begin{align} D = \frac{(\boldsymbol{\tau}\cdot\boldsymbol{\tau}):\mathbf{c}^2 - \boldsymbol{\tau}:\mathbf{c}^4:\boldsymbol{\tau}}{\boldsymbol{\tau}:\boldsymbol{\tau}/5}, \end{align}

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 \lt \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 Γ 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, Reference De La Chapelle, Castelnau, Lipenkov and Duval1998; Hamann and others, Reference Hamann, Weikusat, Azuma and Kipfstuhl2007). 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 initially-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 Γ₀ 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, Reference Faria2001) (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, Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen2021; Richards and others, Reference Richards, Pegler, Piazolo and Harlen2021)

(11)\begin{align} \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), \end{align}

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 orientation space gives the usual mass density of glacier ice $\rho=\int_{S^2} \psi \mathrm{d}{\Omega}$, where $\mathrm{d}{\Omega}=\sin(\theta)\mathrm{d}{\theta}\mathrm{d}{\phi}$ is the infinitesimal solid angle. Note that $\psi_l^{-m} = (-1)^m (\psi_l^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, Reference Svendsen and Hutter1996; Richards and others, Reference Richards, Pegler, Piazolo and Harlen2021)

(12)\begin{align} \frac{\mathrm{D} \psi}{\mathrm{D}t} = -\boldsymbol{\nabla}_{S^2}\cdot(\psi \dot{\mathbf{c}}) + \Gamma \psi, \end{align}

where $\nabla_{S^2}$ is the gradient on S ² and ${\mathrm{D}}/{\mathrm{D}t}$ is the material derivative. The harmonic expansion (11) transforms the fabric model (12) into a high-dimensional matrix problem (Rathmann and Lilien, Reference Rathmann and Lilien2021; Rathmann and others, Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen2021),

(13)\begin{align} \frac{\mathrm{D} \mathbf{s}}{\mathrm{D}t} = (\mathbf{M}_\mathrm{LROT}+\mathbf{M}_\mathrm{DDRX})\cdot\mathbf{s}, \end{align}

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

(14)\begin{align} \mathbf{s} = [\psi_0^0,\psi_2^{0},\psi_2^{1},\psi_2^{2},\psi_4^{0},\cdots,\psi_4^{4},\cdots,\psi_L^{0},\cdots,\psi_L^{L}] \in \mathbb{C}^{(L + 2)^2/4} \end{align}

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 $\overline{\mathbf{s}}(x,y,t)$, and implicitly assuming that the velocity field and gradients refer to their horizontal xy parts, $\overline{\mathbf{s}}$ is governed by

(15)\begin{align} \frac{\text{D}{\overline{\mathbf{s}}}}{\text{D}t} = (\mathbf{M}_\mathrm{LROT}+\mathbf{M}_\mathrm{DDRX})\cdot\overline{\mathbf{s}} + \frac{a_{\mathrm{sfc}}}{H}(\mathbf{s}_{\mathrm{sfc}} - \overline{\mathbf{s}}) + \frac{a_{\mathrm{sub}}}{H}(\mathbf{s}_{\mathrm{sub}} - \overline{\mathbf{s}}), \end{align}

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 $\overline{\mathbf{s}}$ to tend towards the characteristic fabric states of ice that accumulates on the surface $\mathbf{s}_{\mathrm{sfc}}$ or subglacially $\mathbf{s}_{\mathrm{sub}}$ (likely isotropic), depending on the positively-defined accumulation rates $a_{\mathrm{sfc}}$ and $a_{\mathrm{sub}}$.

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.”.

2.4. Viscous anisotropy

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

We follow Rathmann and Lilien Reference Rathmann and Lilien2021 and Lilien and others Reference Lilien2023 who calculated E_ij by substituting $\dot{\boldsymbol{\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:

(16)\begin{align} E_{ij} = (1-\alpha) \frac{ \mathbf{m}_{i}\cdot\dot{\boldsymbol{\epsilon}}^{\mathrm{Sachs}} (\hat{\boldsymbol{\tau}})\cdot \mathbf{m}_{j} }{ \mathbf{m}_{i}\cdot\dot{\boldsymbol{\epsilon}}^{\mathrm{Sachs}}_{\mathrm{iso}}(\hat{\boldsymbol{\tau}})\cdot \mathbf{m}_{j} } + \alpha \frac{ \mathbf{m}_{i}\cdot\dot{\boldsymbol{\epsilon}}^{\mathrm{Taylor}} (\hat{\boldsymbol{\tau}})\cdot \mathbf{m}_{j} }{ \mathbf{m}_{i}\cdot\dot{\boldsymbol{\epsilon}}^{\mathrm{Taylor}}_{\mathrm{iso}}(\hat{\boldsymbol{\tau}})\cdot \mathbf{m}_{j} } , \end{align}

where $\alpha\in[0;1]$. The Sachs and Taylor contributions are $\dot{\boldsymbol{\epsilon}}^{\mathrm{Sachs}} = \langle{\dot{\boldsymbol{\epsilon}}'(\boldsymbol{\tau})}\rangle$ and $\dot{\boldsymbol{\epsilon}}^{\mathrm{Taylor}} = \langle{\boldsymbol{\tau}'(\dot{\boldsymbol{\epsilon}})}\rangle^{-1}$, where $\boldsymbol{\tau}'$ and $\dot{\boldsymbol{\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

(17)\begin{align} \langle{\cdot}\rangle = \frac{1}{\rho}\int_{S^2} (\cdot)\, \psi\ \mathrm{d}{\Omega} . \end{align}

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

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

(18)\begin{align} \dot{\boldsymbol{\epsilon}}' &= A' \Bigg( \boldsymbol{\tau}' - \frac{E'_{cc}-1}{2}(\boldsymbol{\tau}' : \mathbf{c}^2)\mathbf{I} +\frac{3(E'_{cc}-1)-4(E'_{ca}-1)}{2} \boldsymbol{\tau}' : \nonumber\\ & \quad\langle{\mathbf{c}^4}\rangle +(E'_{ca}-1)(\boldsymbol{\tau}'\cdot\mathbf{c}^2 + \mathbf{c}^2\cdot\boldsymbol{\tau}') \Bigg) , \end{align}

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

The free parameters of our micromechanical model (Aʹ, $E'_{cc}$, $E'_{ca}$, α) 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{\boldsymbol{\epsilon}}^{\mathrm{Sachs}}$ and $\dot{\boldsymbol{\epsilon}}^{\mathrm{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, Reference Rathmann2022), implying that choosing a small truncation L is sufficient as long as regularization does not affect harmonics of degree $l\leq 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, Reference Duval, Ashby and Anderman1983). 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\geq 10$ and therefore $\dim(\mathbf{s})\geq 36$ (Rathmann and others, Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen2021).

3.1. Biases

The hypothesis that viscous anisotropy can be represented by (reduced to) a scalar enhancement E implies that 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

(20)\begin{align} \dot{\epsilon}_\mathrm{e,iso} = \sqrt{\frac{\dot{\boldsymbol{\epsilon}}:\dot{\boldsymbol{\epsilon}}}{2}}, \end{align}

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.

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 ϑ between the fabric principal direction m₁ 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 ₁₂ 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 noncoaxial behavior is permitted by the orthotropic rheology, unlike the Glen isotropic rheology where τ and $\dot{\boldsymbol{\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. ϵ_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 τ 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 compression/tension 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}_\mathrm{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 τ (or $\dot{\boldsymbol{\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 (Reference Placidi, Greve, Seddik and Faria2010) 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 (Reference Graham, Morlighem, Warner and Treverrow2018); McCormack and others (Reference McCormack, Warner, Seroussi, Dow, Roberts and Treverrow2022)) 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 (Reference Placidi, Greve, Seddik and Faria2010) 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_{\mathrm{c}}$ and $E_{\mathrm{s}}$ to these minimally and maximally favorable states, respectively, whereas for intermediate states a semi-empirical scaling is proposed based on Azuma (Reference Azuma1995):

(21)\begin{align} E = \begin{cases} (1-E_{\mathrm{c}})\langle{D}\rangle^{8/21(E_{\mathrm{s}}-1)/(1-E_{\mathrm{c}})} + E_{\mathrm{c}}\ \text{for}\ 0 \leq \langle{D}\rangle \leq 1, \\ \dfrac{4\langle{D}\rangle^2(E_{\mathrm{s}}-1) + 25 -4E_{\mathrm{s}}}{21}\ \text{for}\ 1 \lt \langle{D}\rangle \leq 5/2. \end{cases} \end{align}

Here, we follow Placidi and others Reference Placidi, Greve, Seddik and Faria(2010) by setting $E_{\mathrm{c}}=0.1$ and $E_{\mathrm{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, τ can be exchanged for $\dot{\boldsymbol{\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}_\mathrm{e}$ in (4) is a useful quantity to consider, which fulfills the effective Norton–Bailey creep law (Naumenko and Altenbach, Reference Naumenko and Altenbach2007, p. 21) for both isotropic and orthotropic rheologies:

(22)\begin{align} \dot{\epsilon}_\mathrm{e}=A\tau_{\mathrm{e}}^n , \end{align}

where $\tau_{\mathrm{e}}$ is the effective deviatoric stress. The motivation for considering $\dot{\epsilon}_\mathrm{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}_\mathrm{e}$ satisfies the requirement that if deformation is (say) simple shear aligned with ${\mathbf m}_2$ and ${\mathbf m}_3$ so that $\dot{\boldsymbol{\epsilon}} \propto \mathbf{P}_4$, then $\dot{\epsilon}_\mathrm{e}$ is dominated by the term involving the corresponding shear enhancement E ₂₃ since $\mathbf{P}_4:\dot{\boldsymbol{\epsilon}}$ is the only nonzero strain-rate projection in (6).

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

(23)\begin{align} \frac{\dot{\epsilon}_\mathrm{e}}{\dot{\epsilon}_\mathrm{e,iso}} = \sqrt{\frac{\eta_4(\mathbf{P}_4:\dot{\boldsymbol{\epsilon}})^2}{\dot{\boldsymbol{\epsilon}}:\dot{\boldsymbol{\epsilon}}}} = E^{-1/(n+1)}_{23} . \end{align}

Hence $\dot{\epsilon}_\mathrm{e}/\dot{\epsilon}_\mathrm{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)

(24)\begin{align} E = \left(\frac{\dot{\epsilon}_\mathrm{e}}{\dot{\epsilon}_\mathrm{e,iso}}\right)^{q}, \end{align}

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_{\mathrm{e}} \lt \tau_{\mathrm{e,iso}} \Rightarrow \dot{\epsilon}_\mathrm{e} \lt \dot{\epsilon}_\mathrm{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}_\mathrm{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 for several deformation kinematics 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}\mathrm{C})$, and $\Gamma_0=\Gamma_0({-5}^{\circ}\mathrm{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_{\mathrm{c}}=0.2$ and $E_{\mathrm{s}}=3.2$ enables CAFFE to approximately reproduce the behavior of EIE, henceforth referred to as CAFFE $^\dagger$ (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).

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_ij (black line). MODF insets show the modeled fabrics at selected strains.

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, Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen2021). 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.

Regarding the biases mentioned above, caused by fabric incompatibility 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 enhancement (solid lines approximately agree in Figure 6), but significantly overestimate longitudinal strain rates when the stress state is a superposition of shear and tension (dashed purple and green lines disagree with the dashed black line in Figure 6b).

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}\,\mathrm{km}$, length $L={100}\,\mathrm{km}$, and uniform initial thickness $H_0={500}\,\mathrm{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,

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

Figure 8. Schematic of the idealized half-width ice shelf model.

where 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 ₀. 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, Reference Huth, Duddu and Smith2021) 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}\mathrm{C}$ to make the effect of fabric on modeled viscosity unambiguous. The canonical rate factor A(T) by Cuffey and Paterson (Reference Cuffey and Paterson2010, p. 72–74) is used, but we recognize that recent work (Fan and others, Reference Fan, Wang, Prior, Breithaupt, Hager and Wallis2025) 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 (Reference Cuffey and Paterson2010, 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, Reference Behn, Goldsby and Hirth2021) despite recent work suggesting that $n\simeq 4$ (Bons and others, Reference Bons2018; Ranganathan and Minchew, Reference Ranganathan and Minchew2024; Fan and others, Reference Fan, Wang, Prior, Breithaupt, Hager and Wallis2025) or n < 3 (Wang and others, Reference Wang, Lai, Prior and Cowen-Breen2025c) 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 \overline{\mathbf{s}}$ is added to (15) for numerical stability following Rathmann and Lilien Reference Rathmann and Lilien2021 and Lilien and others Reference Lilien2023 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.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}\mathrm{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.

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 $^\dagger$ 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).

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 \mathrm{d}{\Omega}$ (Ismail and Mainprice, Reference Ismail and Mainprice1998) (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 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.

Figure 10. Orthotropic model results in steady state when DDRX is the dominant crystal process. Caption same as in Figure 9.

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 horizontal (Cartesian) components of E_ij expected to be most important for the simulated flows (Figure 8). In the model without DDRX, both EIE and CAFFE $^\dagger$ 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 fabric-induced hardening is somewhat underestimated and that CAFFE $^\dagger$ 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 $^\dagger$ underestimating E for DDRX fabrics in simple shear (Figure 7d,g). Whether there are more optimal values of $E_{\mathrm{c}}$ and $E_{\mathrm{s}}$ that improve the CAFFE $^\dagger$ misfit for DDRX fabrics was not explored.

Note that we considered CAFFE $^\dagger$, rather than CAFFE (no dagger), for these experiments, since comparing misfits is meaningful only if $E_{\mathrm{c}}$ and $E_{\mathrm{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).