1. Introduction
Ice varies from isotropic to anisotropic, with crystals developing a preferred orientation fabric due to ice flow dynamics. Ice crystals tend to align with the direction of an applied force, and shearing of ice is enhanced along basal planes perpendicular to the applied force (Alley, Reference Alley1988). The orientation of the basal planes is described by a vector referred to as the c-axis. The snow grain orientation is initially isotropic, and anisotropic ice fabric develops at a rate dependent on a variety of physical factors including the temperature, pressure and strain rate of the ice. These dependencies allow the investigation of both the steady-state and transient behaviour of ice through the analysis of observational fabric anisotropy data and ice flow models.
Ice cores have long been used to quantify the direction of ice crystals and the degree of fabric anisotropy with depth (Durand and others, Reference Durand, Gagliardini, de Fleurian, Zwinger and Le Meur2009; Montagnat and others, Reference Montagnat2014; Weikusat and others, Reference Weikusat2017). However, the lack of horizontal spatial information hinders the study of the relationship between ice flow and fabric anisotropy. Apart from ice cores sites, anisotropic information is sparse and relies on geophysical methods such as seismics, both active (Diez and others, Reference Diez2014; Brisbourne and others, Reference Brisbourne, Martín, Smith, Baird, Kendall and Kingslake2019) and passive (Smith and others, Reference Smith2017; Kufner and others, Reference Kufner2023) or radar (Fujita and others, Reference Fujita, Maeno and Matsuoka2006; Drews and others, Reference Drews, Eisen, Steinhage, Weikusat, Kipfstuhl and Wilhelms2012; Matsuoka and others, Reference Matsuoka, Power, Fujita and Raymond2012). The main limitation is that, generally, seismic methods only provide depth-averaged fabric anisotropy and radar only provides fabric anisotropy information in the horizontal plane, perpendicular to the direction of wave propagation. Recent advances in phase coherent radar systems and data processing for use in inferring ice fabric anisotropy have resulted in an increase in the acquisition of observational fabric data (Dall, Reference Dall2010; Young and others, Reference Young, Martín, Christoffersen, Schroeder, Tulaczyk and Dawson2021a; Reference Young2021b; Ershadi and others, Reference Ershadi2022; Jordan and others, Reference Jordan, Martín, Brisbourne, Schroeder and Smith2022). These surveys provide the spatial variability of fabric anisotropy in more extensive areas and aid understanding of the relation between fabric anisotropy and ice flow. Prior to this, the lack of observational data to compare with and constrain models has long hindered the progression of ice fabric modelling.
Although we concentrate only on the coupling of the fabric anisotropy field to the velocity and stress fields in the present study, we note that advances in a full coupling of crystal orientation tensor evolution with viscosity are ongoing with challenges remaining regarding numerical instability (Gerber and others, Reference Gerber2023). Several modelling studies investigating the effects of an anisotropic rheology on ice flow dynamics at ice divides in two dimensions (Martin and others, Reference Martín, Hindmarsh and Navarro2009; Martin and Gudmundsson, Reference Martín and Gudmundsson2012; Lilien and others, Reference Lilien2023) and in ice streams (Gerber and others, Reference Gerber2023; Richards and others, Reference Richards, Pegler, Piazolo, Stoll and Weikusat2023) have shown the importance of including fabric anisotropy or viscous anisotropy in ice flow models. In idealised, two-dimensional simulations, Richards and others Reference Richards, Pegler and Piazolo(2022) investigate the influence of the strain rate and spin on the fabric type. Building on the work of Thorsteinsson and others (Reference Thorsteinsson, Waddington and Fletcher2003); Gillet-Chaulet and others (Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006) and Pettit and others (Reference Pettit, Thorsteinsson, Jacobson and Waddington2007); Rathmann and others (Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen2021) investigated the coupling between fabric anisotropy and viscosity through enhancement factors which result in greater shearing perpendicular to the predominant c-axis direction in two-dimensional numerical ice-flow simulations. Using a computationally efficient anisotropic flow model with a simplified representation of fabric anisotropy coupled to the viscosity via enhancement factors, McCormack and others Reference McCormack, Warner, Seroussi, Dow, Roberts and Treverrow(2022) studied the effect of fabric anisotropy in larger-scale three-dimensional models.
The coupling of fabric anisotropy evolution in ice flow models is not without complication due to challenges with numerical stability, the parameter choice within the fabric anisotropy evolution equation and the coupling with viscosity in the flow law. Furthermore, because an isotropic assumption has traditionally been employed in ice sheet models, model parameters have been optimised to fit isotropic rather than anisotropic ice flow dynamics. For these reasons, a thorough investigation of each step in the process of modelling ice fabric anisotropy is needed.
In this paper, we present the first model of the three-dimensional, diagnostic, fabric anisotropy field of an ice rise using the finite element model Elmer/Ice (Gagliardini and others, Reference Gagliardini2013) across 100 partitions and applied to Derwael Ice Rise (DIR). We investigate the dependence of the fabric anisotropy field on the strain rate and deviatoric stress fields without coupling with the viscosity and without re-crystallisation terms, which allows analysis in the absence of the additional feedback complexity. Our simulations provide novel ice fabric predictions for a number of three-dimensional ice-flow settings including areas of vertical-shear-dominated flank flow, grounding zones and shear zones. We highlight areas where geophysical measurements of ice-fabric types would be most informative to constrain relevant model parameters. Furthermore, the direction of the c-axis across the ice rise is investigated to identify regions where a vertical c-axis assumption is valid in radar measurements of fabric anisotropy, providing a novel method for determining the error in eigenvalues with such an assumption.
1.1. Motivation
Our study is motivated by a lack of progress in large-scale ice-sheet models considering crystal orientation fabric. This is despite observations of strong crystal orientation fabric in polar ice (Alley, Reference Alley1988), knowledge of how anisotropic polar ice is (Duval and others, Reference Duval, Ashby and Anderman1983) and field evidence of the effect crystal orientation fabric has on ice flow (Gerber and others, Reference Gerber2023). After the introduction of the crystal orientation tensor (Gödert, Reference Gödert2003), the flow-induced crystal orientation fabric evolution can be considered by large-scale ice-sheet models (Gagliardini and others, Reference Gagliardini2013).
We believe that the main reasons for the lack of progress are (a) the numerical implementation of crystal orientation fabric evolution is challenging (Seddik and others, Reference Seddik, Greve, Zwinger and Placidi2011), (b) the interpretation of model output and comparison with non-comprehensive observations is complex (Jordan and others, Reference Jordan, Martín, Brisbourne, Schroeder and Smith2022) and (c) essential model parameters within the theory are not yet constrained by observations (Ma and others, Reference Ma, Gagliardini, Ritz, Gillet-Chaulet, Durand and Montagnat2010).
Numerically, the main issue is that numerical dispersion tends to break down the orientation tensor properties. Here, we present a numerically robust, three-dimensiontal model, discuss a framework to interpret model output for comparison with observations and highlight areas where observations could better constrain model parameters.
1.2. DIR
Ice rises are an ideal study location for understanding ice-flow processes, because transitions between different flow regimes occur over comparatively short spatial scales. Furthermore, ice rises regulate flow from the Antarctic Ice Sheet towards the open ocean, controlling ice shelf velocities and the continental grounding line position (Favier and Pattyn, Reference Favier and Pattyn2015; Schannwell and others, Reference Schannwell, Drews, Ehlers, Eisen, Mayer and Gillet-Chaulet2019; Reference Schannwell2020; Henry and others, Reference Henry, Drews, Schannwell and Višnjević2022). Formation, evolution and disintegration of ice rises occur over glacial-interglacial timescales, meaning that remnants of ice properties such as temperature and fabric anisotropy from previous flow regimes may become stored in the slow flowing ice of an ice rise.
Ice rises typically have clear ice divides transitioning into a flank-flow regime on all sides with little to no basal sliding (Matsuoka and others, Reference Matsuoka2015). The grounding line is typically located a few tens of kilometres away from the divide, and the flow field transitions to the surrounding ice shelves through narrow shear zones with large horizontal shear strain rates. We chose DIR (Fig. 1) in East Antarctica as our study site for two reasons. First, we can rely on a predicted three-dimensional steady-state velocity field developed in a previous study (Henry and others, Reference Henry, Schannwell, Višnjević, Millstein, Bons and Eisen2023). This velocity field is based on a transient, thermomechanically coupled full Stokes model with an isotropic rheology. The model was forced with an observationally constrained surface mass balance field, and predictions of the model output were validated with extensive radar observations which are available for this Ice Rise (Koch and others, Reference Koch2024). Second, previous studies suggest that DIR is likely close to steady state (Callens and others, Reference Callens, Drews, Witrant, Philippe and Pattyn2016) possibly with a minor amount of thinning (Drews and others, Reference Drews, Matsuoka, Martín, Callens, Bergeot and Pattyn2015). This justifies the steady-state assumption applied here. DIR is a marine-based ice rise and has a flow divide in the form of a curved arc extended south-west to north-east. The maximum ice thickness is roughly 630 m, and flow velocities in the flank are typically slower than 10 m a−1.
In (a), an overview is shown of DIR, with the surrounding ice shelves and the grounding line (Jezek and others, Reference Jezek, Curlander, Carsey, Wales and Barry2013; Rignot and Scheuchl., Reference Rignot and Scheuchl2017). The velocity field is shown in colour. In (b), the simulated upper surface velocity field of DIR is shown, based on simulations in Henry and others Reference Henry, Schannwell, Višnjević, Millstein, Bons and Eisen(2023). Boxes A and B show the areas referred to as the areas of high horizontal divergence at the tails of the flow divide.
2. Methods
The simulations use output from Henry and others Reference Henry, Schannwell, Višnjević, Millstein, Bons and Eisen(2023) as a starting point in which the finite element model Elmer/Ice (Gagliardini and others, Reference Gagliardini2013) was applied to DIR. Most importantly, we use the predicted three-dimensional, steady-state velocity field to make predictions of the fabric anisotropy resulting from a one-way coupling with the crystal orientation evolution equations detailed below.
2.1. Governing equations
The equations of motion for Stokes flow are written as
\begin{equation}
\boldsymbol{\nabla}\cdot (\boldsymbol{\tau} - P \textbf{I}) + \rho_i \textbf{g} = 0,
\end{equation} where τ is the deviatoric stress tensor, P is the pressure, I is the identity matrix, ρi is the ice density and 
$\textbf{g} = g \hat{\mathbf{e}}_z$ is the gravitational acceleration. The ice is subject to an incompressibility condition,
\begin{equation}
\boldsymbol{\nabla}\cdot \textbf{u} = 0,
\end{equation}and Glen’s flow law,
\begin{equation}
\boldsymbol{\tau} = 2 \eta \dot{{\varepsilon}},
\end{equation} when assuming that ice viscosity is isotropic. This relation describes the nonlinear dependence between the strain rate tensor, 
$\dot{{\varepsilon}}$, and the deviatoric stress tensor. The effective viscosity, η, is
\begin{equation}
\eta = \frac{1}{2} EA(T')^{-1/n} \dot{\varepsilon}_e^{(1-n)/n},
\end{equation} where E is an enhancement factor which is spatially and temporally constant here and 
$A(T')$ is the ice fluidity which is dependent on the ice temperature relative to the pressure melting point, Tʹ, which is solved using a temperature evolution equation as described in Henry and others Reference Henry, Schannwell, Višnjević, Millstein, Bons and Eisen(2023). Here, n is Glen’s flow law exponent, and 
$\dot{\varepsilon}_e = \sqrt{\mathrm{tr} \dot{{\varepsilon}}^2/2}$ is the effective strain rate.
When assuming a direction-dependent viscosity, an alternative flow law is defined as
\begin{equation}
\tau = 2 \hat{\boldsymbol{\eta}} \dot{\varepsilon},
\end{equation} where 
$\hat{\boldsymbol{\eta}}$ is an anisotropic enhancement tensor. In this study, this equation is calculated using the GOLF model (Castelnau and others, Reference Castelnau, Duval, Lebensohn and Canova1996; Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005; Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006; Gagliardini and others, Reference Gagliardini2013) although other methods exist (e.g. Rathmann and others Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen(2021)). Such formulations, coupled with an equation describing the evolution of the crystal orientation, can be unstable in fast-flowing areas (Gerber and others, Reference Gerber2023) such as in the ice shelf in the model domain presented here. We therefore use Eq. (3) rather than Eq. (5) as the constitutive law in our model setup. However, a direction-dependent viscosity is used for calculating the crystal orientation evolution, as outlined below.
Keeping the velocity and stress fields constant in time, a semi-Lagrangian fabric anisotropy evolution equation is coupled and simulations are performed for 20000 years. To initialise the fabric anisotropy simulation, the simulation named n3E0.5dsdt50 in Henry and others Reference Henry, Schannwell, Višnjević, Millstein, Bons and Eisen(2023) is used. The simulation time of 20000 years was deemed appropriate given that in Henry and others Reference Henry, Schannwell, Višnjević, Millstein, Bons and Eisen(2023), it was predicted that the ice in DIR is roughly 8000 years old at a depth of 95 %.
The distribution of crystal orientation in a volume of ice can be described by crystal orientation tensors of varying complexity (Advani and Tucker, Reference Advani and Tucker1987). Higher-order tensors allow ever-finer details of the crystal orientation distribution to be captured. Such complexity is computationally infeasible, and thus, the description of the crystal orientation tensor is reduced to the second-order crystal orientation tensor, 
$\mathbf{a}^{(2)}$, and the fourth-order tensor, 
$\boldsymbol{a}^{(4)}$, thereby neglecting higher-order tensors (Cintra and Tucker, Reference Cintra and Tucker1995; Chung and Kwon, Reference Chung and Kwon2002). Here, c is the c-axis unit vector and the operator 
$\langle \rangle$ denotes the average over all the grains that compose the ice polycrystal. In reality, 
$\boldsymbol{a}^{(2)}$ is dependent on 
$\boldsymbol{a}^{(4)}$, which is dependent on 
$\boldsymbol{a}^{(6)}$ and so on.
It is useful to compute the eigenvalues and eigenvectors of the second-order crystal orientation tensor, 
$\boldsymbol{a}^{(2)}$, in order to understand the fabric anisotropy type. The eigenvalues of 
$\boldsymbol{a}^{(2)}$, namely, λ 1, λ 2 and λ 3, and the corresponding eigenvectors, v1, v2 and v3, of the crystal orientation tensor 
$\mathbf{a}^{(2)}$ satisfy
\begin{equation}
\mathbf{a}^{(2)} \mathbf{v} = \lambda \mathbf{v},
\end{equation} defining each eigenvalue by its subscript such that 
$\lambda_1 \leqslant \lambda_2 \leqslant \lambda_3$. In the case of randomly oriented c-axes or isotropic ice, all three eigenvalues are similar in size such that 
$\lambda_1 \approx \lambda_2 \approx \lambda_3 \approx 1/3$. In areas where 
$\lambda_3 \geqslant\lambda_1 \approx \lambda_2$, ice fabric is said to have a single maximum fabric with the majority of c-axes pointing in a single direction. Where 
$\lambda_2 \approx \lambda_3 \geqslant\lambda_2$, ice is said to have a girdle fabric, with c-axes following an arc or circle. See the Appendix for an explanation of the connection to observational data.
The evolution of the second-order orientation tensor can be described by
\begin{equation}
\begin{aligned}
\bigg(\frac{\partial}{\partial t} + \mathbf{u} \cdot \nabla \bigg) \mathbf{a}^{(2)} & = \mathbf{W}\mathbf{a}^{(2)} - \mathbf{a}^{(2)}\mathbf{W} \\
& \phantom{=} - \iota (\mathbf{C}\mathbf{a}^{(2)} + \mathbf{a}^{(2)}\mathbf{C} - 2 \mathbf{a}^{(4)} : \mathbf{C}),
\end{aligned}
\end{equation}where the tensor C is defined by and implemented in our model as
\begin{equation}
\mathbf{C} = (1 - \alpha) \dot{\boldsymbol{\varepsilon}} + 2 \alpha \boldsymbol{\hat{\eta}} \boldsymbol{\dot{\varepsilon}}.
\end{equation} Here, α is the so-called interaction parameter, which describes the relative influence of the strain rate and stress tensors, W is the spin of the ice and 
$\boldsymbol{\hat{\eta}}$ is the anisotropic enhancement tensor calculated using GOLF (Castelnau and others, Reference Castelnau, Duval, Lebensohn and Canova1996; Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005; Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006; Gagliardini and others, Reference Gagliardini2013). The parameter ι determines the rate at which the crystal orientation tensor influenced by the weighted combination of strain rate and deviatoric stress tensors. The last term, 
$2 \alpha \boldsymbol{\hat{\eta}} \boldsymbol{\dot{\varepsilon}}$, in Eq. (8) is equivalent to 
$\alpha k_s A \tau_e^{n-1} \boldsymbol{\tau}$ in Gagliardini and others Reference Gagliardini(2013), so that
\begin{equation}
\mathbf{C} = (1 - \alpha) \dot{\boldsymbol{\varepsilon}} + \alpha k_s A \tau_e^{n-1} \boldsymbol{\tau}.
\end{equation} We assume that 
$k_s=1$, as no given value was published in the studies which we base our parameter choices on (Martin and others, Reference Martín, Hindmarsh and Navarro2009; Seddik and others, Reference Seddik, Greve, Zwinger and Placidi2011; Martin and Gudmundsson, Reference Martín and Gudmundsson2012; Gagliardini and others, Reference Gagliardini2013). Values in other literature have ranged from 
$k_s \lt 2.5$ (Gagliardini and Meyssonnier, Reference Gagliardini and Meyssonnier2002) to 
$k_s=10$ (Ma and others, Reference Ma, Gagliardini, Ritz, Gillet-Chaulet, Durand and Montagnat2010). For each time step, iterations are performed by repeatedly calculating the fabric anisotropy and the direction-dependent viscosity which would result in that fabric anisotropy. This results in a fabric anisotropy which is dependent on a direction-dependent viscosity through the last term in Eq. (8).
Early work recognised the influence of the cumulative strain and stress on ice crystal c-axis orientation (Alley, Reference Alley1988), which led to the development of a crystal orientation tensor evolution equation (Eq. (7)) dependent on the velocity gradient through the spin and strain rate tensors and the deviatoric stress tensor (Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006). Although it is known that the velocity and stress fields have an influence on the fabric anisotropy field of ice, it remains unclear what the relative influence is. The fabric anisotropy evolution equation (Eq. (7)) is made up of various terms, with the spin tensor, W, acting to rotate the crystal orientation tensor to follow the spin of the ice. Acting opposite to the spin tensor is the tensor, C, which is a weighted combination of the strain rate tensor, 
$\boldsymbol{\dot{\varepsilon}}$, and the deviatoric stress tensor, τ. This combination comes from the fact that the behaviour of macroscopic materials can be limited by two extreme approximations: uniform stress, where the stress in the crystals is assumed to be identical to the macroscopic stress, and Taylor or uniform strain rate, where the strain rate in the crystals is assumed to be identical to the macroscopic strain rate (Gagliardini and others, Reference Gagliardini, Gillet-Chaulet and Montagnat2009). The choice of α and ι in previous studies has been motivated by assumptions of the relative influence of stress and strain rate on crystal orientation evolution. The last term in the equation describes the influence of the higher-order crystal orientation tensor, 
$\boldsymbol{a}^{(4)}$ on the 
$3 \times 3$ crystal orientation tensor, 
$\boldsymbol{a}^{(2)}$.
To solve Eq. (7) we also require a relation between 
$\mathbf{a}^{(2)}$ and 
$\mathbf{a}^{(4)}$, a closure approximation, and we use an invariant-based optimal fitting closure approximation (Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006). The distribution of c-axes can be described by tensors of ever-increasing order, but a compromise is found by using the closure approximation (Advani and Tucker, Reference Advani and Tucker1987; Cintra and Tucker, Reference Cintra and Tucker1995; Chung and Kwon, Reference Chung and Kwon2002) in order to reduce computation time. The effect of the choice of α and ι on the fabric anisotropy field is investigated using combinations of parameters from the previous literature presented in Table 1.
Parameter combinations
The value of α determines the relative influence of the strain rate and deviatoric stress tensors on fabric anisotropy evolution. For example, a value of α = 0 means that there is dependence on the strain rate tensor but no dependence on the deviatoric stress. Alternatively, a value of α = 1 provides dependence on the deviatoric stress tensor but no dependence on the strain rate tensor. The parameter ι adjusts the rate at which the fabric anisotropy field develops in response to the weighted combination of the strain rate and deviatoric stress fields and takes a value between ι = 0 and ι = 1.
Initially, all ice in the model domain is isotropic, described by the tensor 
$\mathbf{a}^{(2)} = \frac{1}{3} \textbf{I}$, where I is the 
$3\times3$ identity matrix. During transient simulation of the fabric anisotropy field dependent on the velocity field, the upper surface is assigned an isotropic boundary conditions on all ice entering the domain due to accumulation.
The fabric anisotropic evolution model described in Eq. (7) together with its boundary and initial conditions is solved using a semi-Lagrangian method as described in Martin and others Reference Martín, Hindmarsh and Navarro(2009). The determinant of 
$\mathbf{a}^{(2)}$ gives information about the degree of crystal orientation of the fabric (Advani and Tucker, Reference Advani and Tucker1987). We define here the degree of crystal orientation with the scalar value,
\begin{equation}
d=1-3^3 \det (\mathbf{a}^{(2)} ).
\end{equation}The degree of crystal orientation varies from zero for isotropic ice to one for single-maximum fabric and is independent of the frame of reference as the determinant is an invariant. We constrain Equation (7) with the condition,
\begin{equation}
\bigg(\frac{\partial}{\partial t} + \mathbf{u} \cdot \nabla \bigg) d \geqslant0,
\end{equation}meaning that ice is constrained to increase in degree of crystal orientation over time. This constraint on the degree of crystal orientation stops numerical dispersion breaking the simulation in areas with high strain rates. Our method is also straightforward to implement in large-scale models that run in parallel environments.
3. Results
3.1. Analysis of simulated fabric anisotropy field
We display predicted eigenvalues for the α = 0, ι = 1 simulation at an elevation of z = 0 which encompasses an ice-depth range of close to 0 m in the ice shelf to over 300 m at the ice rise flow divide where stresses are dominated by vertical compression and lateral extension (Fig. 2). The main characteristic, which is evident in all simulations, is that the ice fabric evolves from an isotropic material at the surface and, under deformation, develops into an anisotropic fabric varying from a single maximum to a girdle fabric (Fig. 3). In the ice rise interior, the largest eigenvalue, λ 3, is much larger than λ 1 and λ 2, particularly at and surrounding the flow divide (Figs. 1b and 2). This results in greater fabric anisotropy here. Further south in the ice rise, where ice flow is hindered by the convergence of flow between the ice rise and the ice shelf, differences are not as substantial. In the ice rise interior away from the grounding line, differences between the smaller two eigenvalues, λ 1 and λ 2 are mostly small, increasing slightly on the stoss side of the ice rise. The stoss side is defined as the side of the ice rise with a flow direction opposing the upstream ice shelf flow direction (Figs. 1b and 2). At the grounding zone on the stoss side of the ice rise, where horizontal convergence of flow occurs, the three eigenvalues are similar in magnitude, whereas at the transition from grounded to floating ice on the lee side of the ice rise, where flow is dominated be extension, a slightly larger difference is seen between the smaller two eigenvalues, λ 1 and λ 2. In the ice shelves away from the grounding zone, there are differing patterns in the relative magnitude of the three eigenvalues. In the ice shelf west of the ice rise where velocities are higher than east of the ice rise and extension dominates, the largest eigenvalue is significantly larger in magnitude compared with the two smaller eigenvalues. In the ice shelf to the east of the ice rise, where extension of flow is an active process, differences in magnitude are generally much lower. The ice shelf south of the ice rise, on the stoss side, shows one eigenvalue much larger than the other two, similar to the magnitudes in the ice shelf on the western side of the ice rise.
The eigenvalues, λ 1, λ 2 and λ 3, of the crystal orientation tensor in the α = 0, ι = 1 simulation at an elevation of z = 0, corresponding to the sea level. The solid lines, black in the plots showing λ 1 and λ 2, and white in the plot showing λ 3, are contours of depth below the upper ice surface, and the dashed lines show the grounding line. The dotted line in the λ 3 figure shows where the cross-section in Figure 3 is taken.
Cross-sections through the flow divide as shown in Figure 2 showing λ 3, the largest eigenvalue for (a) the α = 0, ι = 1 simulation, (b) the α = 0, ι = 0.6 simulation, (c) the α = 1, ι = 1 simulation and (d) the α = 0.06, ι = 1 simulation. The solid lines show isochrones.
3.2. Comparison of simulations with differing parameter choice
In order to compare simulation results across the four combinations of the parameters α and ι in Eq. (7) used in the previous studies stated in Table 1, we use a number of metrics to understand the fabric types evolving in the various flow regimes at DIR. First, we describe the differences in relative magnitude of the eigenvalues of the crystal orientation tensor across the four simulations. If the logarithmic ratio (Woodcock, Reference Woodcock1977), 
$\ln (\lambda_3 / \lambda_2)$, is large, then the c-axes are more concentrated in a single direction. The direction of the largest horizontal crystal orientation eigenvector, 
$\boldsymbol{v}_{2, H}$, shows the predominant horizontal c-axis direction.
The simulations with parameter choices of α = 0, ι = 1 (Fig. 4a) and α = 0.06, ι = 1 (Fig. 4d) have comparable results, with a much greater 
$\ln (\lambda_3 / \lambda_2)$ than predicted by the other two simulations. The simulation results differ, however, in the spatial variation of 
$\ln (\lambda_3 / \lambda_2)$. The α = 0, ι = 1 simulation shows larger values of 
$\ln (\lambda_3 / \lambda_2)$ at the tails of the ice rise flow divide (Boxes A and B in Figure 1b), whereas the α = 0.06, ι = 1 simulation shows slightly greater relative values along the flow divide. Although not negligible, the simulations with parameter choices of α = 0, ι = 0.6 (Fig. 4b) and α = 1, ι = 1 (Fig. 4c) show much lower ratios between the largest and second largest eigenvalues, with the α = 0, ι = 0.6 simulation showing a slightly higher value at the centre of the flow divide. In the α = 1, ι = 1 simulation, a differing pattern of horizontal eigenvectors originating at a point source at the flow divide, whereas the other simulations show alignment of the eigenvector direction along the flow divide. Of note is also the differing eigenvector directions at the tails of the flow divide, with differing patterns of vector divergence and convergence.
The ratio of the two larger eigenvalues of the 
$3 \times 3$ crystal orientation tensor at z = 0 for (a) the α = 0, ι = 1 simulation, (b) the α = 0, ι = 0.6 simulation, (c) the α = 1, ι = 1 simulation and (d) the α = 0.06, ι = 1 simulation. The dashes show the maximum horizontal fabric anisotropy direction. The dashed line shows the grounding line, and the solid line shows contours of the depth below the upper ice surface.
The two smaller eigenvalues, λ 1 and λ 2, are investigated using the metric 
$\ln (\lambda_2 / \lambda_1)$ (Fig. 5). Large ratios between the two smaller eigenvalues indicate that c-axis directions for a given volume of ice are concentrated along an arc. In all simulations, there are high values of 
$\ln (\lambda_2 / \lambda_1)$ at the flow divide and generally low values in the vicinity of the grounding zone on the stoss side of the ice rise, with values differing between simulations elsewhere. The α = 0, ι = 0.6 (Fig. 5b) simulation shows the lowest values overall. The simulations with parameter choices of α = 0, ι = 1 and α = 0.06, ι = 1 show relatively similar results, with the highest eigenvalue ratio at the flow divide and in the areas of the ice rise perpendicular to the flow divide. The α = 1, ι = 1 simulation (Fig. 5c) has areas of a large ratio between λ 2 and λ 1, much larger than any other simulation, and highlights the effect of a strong dependence on the strain-rate tensor on the crystal orientation tensor evolution. Rather than having high values perpendicular to the flow divide, the highest values are along a small band following the flow divide and in the areas of high horizontal divergence at the tails of the flow divide, extending in some areas as far as the grounding line. There are no large differences between λ 2 and λ 1 in the ice shelves, with the largest values being concentrated north-west of the ice rise in all simulations.
The ratio of the two smaller eigenvalues of the 
$3 \times 3$ crystal orientation tensor at z = 0 for (a) the α = 0, ι = 1 simulation, (b) the α = 0, ι = 0.6 simulation, (c) the α = 1, ι = 1 simulation and (d) the α = 0.06, ι = 1 simulation. The dashes show the maximum horizontal fabric anisotropy direction. The dashed line shows the grounding line, and the solid line shows contours of the depth below the upper ice surface.
The Woodcock k value (Woodcock, Reference Woodcock1977), defined as
\begin{equation}
k = \frac{\ln (\lambda_3 / \lambda_2)}{\ln (\lambda_2 / \lambda_1)},
\end{equation} provides a metric by which to investigate which areas are characterised by a single maximum fabric and which are characterised by girdle fabric. Note that k can have values between k = 0 and 
$k = \infty$. In order to better investigate small values of k, we plot 
$\ln(k)$ for each simulation in Table 1. If 
$\ln(k) \gt 0$, then the ice is defined as having a single maximum fabric, and if 
$\ln(k) \lt 0$, the ice is defined as having a girdle fabric. Furthermore, if 
$\ln(k) \gt \gt 0$ and 
$\ln (\lambda_3 / \lambda_2) \gt \gt 0$, then the ice has a strong single maximum. If, on the other hand, 
$\ln(k) \lt \lt 0$ and 
$\ln (\lambda_2 / \lambda_1) \gt \gt 0$, then ice has a strong girdle fabric.
In the α = 0, ι = 1 simulation, there are relatively high 
$\ln(k)$ values at the tails of the flow divide (Fig. 6a), extending in some areas almost to the grounding line. Elsewhere on the ice rise, values of 
$\ln(k)$ are generally close to zero. The α = 0.06, ι = 1 simulation (Fig. 6d) shows similar results except at the north-east of the ice rise, where high values of 
$\ln(k)$ are concentrated closer to the grounding line. The α = 0, ι = 0.6 simulation (Fig. 6b) shows higher values of 
$\ln(k)$ even further from the flow divide and a small area with negative 
$\ln(k)$ values at the north-eastern end of the flow divide. The α = 1, ι = 1 simulation (Fig. 6c) shows, by far, the most negative values of 
$\ln(k)$, concentrated at the two tails of the flow divide and small areas with positive values of 
$\ln(k)$ perpendicular to the flow divide. All four simulations show an almost-continuous band of negative values of 
$\ln(k)$ at the grounding line or a small distance away from the grounding line in the ice shelf. Moving away from the grounding line, a general increase in values of 
$\ln(k)$ is seen, with some exceptions.
The logarithm of the Woodcock k value at an elevation of z = 0 for (a) the α = 0, ι = 1 simulation, (b) the α = 0, ι = 0.6 simulation, (c) the α = 1, ι = 1 simulation and (d) the α = 0.06, ι = 1 simulation. The dashes show the maximum horizontal fabric anisotropy direction. The dashed line shows the grounding line, and the solid line shows contours of the depth below the upper ice surface.
3.3. Metrics for comparison with radar data
In quad-polarimetric radar processing, it is often assumed that because of the dominance of vertical compression, one eigenvector aligns with the vertical direction. In areas where this assumption holds, signal processing can be simplified (Ershadi and others, Reference Ershadi2022; Jordan and others, Reference Jordan, Martín, Brisbourne, Schroeder and Smith2022). Obliquely oriented fabric types can, in theory, also be detected (Matsuoka and others, Reference Matsuoka, Wilen, Hurley and Raymond2009; Rathmann and others, Reference Rathmann, Lilien, Grinsted, Gerber, Young and Dahl-Jensen2022), but thus far this has not been done for observations which are typically only collected in a nadir-viewing geometry. Here, we present results evaluating the applicability of the assumption of one vertical eigenvector across all four simulations.
In general, the predicted tilt angle is small in the grounded area across all simulations (Fig. 7). Note that the grounded area is the area within the dashed line marking the grounding line. The angle remains small, with similar spatial patterns in the simulations with parameter choices of α = 0, ι = 1 and α = 0.06, ι = 1 (Fig. 7a and d, respectively). In the α = 0, ι = 0.6 simulation (Fig. 7b), differences between the z-direction and v3 are generally small, with the exception of slightly higher values on either side of the flow divide. Spatial patterns on the stoss side of the ice rise are similar to those in the above-mentioned α = 0, ι = 1 and α = 0.06, ι = 1 simulations. The simulation with the largest tilt angles is the α = 1, ι = 1 simulation (Fig. 7c), with differences of 
$\pi/8$ to 
$\pi/4$ radians at and perpendicular to the flow divide. In all simulations, the tilt angle in the ice shelf is larger a short distance away from the grounding line, except for small, isolated areas.
The angle between the eigenvector corresponding to the largest eigenvalue and the vertical direction at an elevation of z = 0 corresponding to the sea level for (a) the α = 0, ι = 1 simulation, (b) the α = 0, ι = 0.6 simulation, (c) the α = 1, ι = 1 simulation and (d) the α = 0.06, ι = 1 simulation. The dashed line shows the grounding line, and the solid lines show contours for the depth below the surface in metres.
We furthermore calculate the eigenvectors and eigenvalues of the horizontal crystal orientation tensor, i.e. the upper left 
$2 \times 2$ part of the 
$3 \times 3$ tensor 
$\mathbf{a}^{(2)}$. The reason for this is that if the 
$3 \times 3$ crystal orientation tensor has one strictly vertical eigenvector, the eigenvalues and eigenvectors of the 
$2 \times 2$ tensor correspond to the horizontal eigenvalues and eigenvectors of the 
$3 \times 3$ tensor. We denote the eigenvalues of the 
$2 \times 2$ tensor by 
$\lambda_{1,H}$ and 
$\lambda_{2,H}$, and the corresponding eigenvectors by 
$\mathbf{v}_{1, H}$ and 
$\mathbf{v}_{2,H}$, respectively. By comparing the eigenvalues of the 
$2 \times 2$, horizontal crystal orientation tensor with the eigenvalues of the 
$3 \times 3$ crystal orientation tensor, an error estimate can be found for assuming that the eigenvector corresponding to the largest eigenvalue is aligned with the z-axis.
We present percentage differences between λ 1 and 
$\lambda_{1,H}$ (Fig. 8), as well as λ 2 and 
$\lambda_{2,H}$ (Fig. 9). We find that if a difference exists between λ 1 and 
$\lambda_{1,H}$, or λ 2 and 
$\lambda_{2,H}$, then the 
$2 \times 2$ tensor eigenvalues always underestimate the smaller two 
$3 \times 3$ tensor eigenvalues. Furthermore, percentage differences tend to be higher for the smaller eigenvalues λ 1 and 
$\lambda_{1,H}$ than for λ 2 and 
$\lambda_{2,H}$, but exceptions to this are seen. The α = 0, ι = 1 simulation shows a negligible percentage difference between λ 1 and 
$\lambda_{1,H}$ (Fig. 8a) and small percentage differences for λ 2 and 
$\lambda_{2,H}$ of up to 
$~10$ % (Fig. 9a). A similar spatial distribution of percentage differences appears for λ 2 and 
$\lambda_{2,H}$ in the α = 0.06, ι = 1 simulation (Fig. 9d). Most notable is the percentage difference between λ 1 and 
$\lambda_{1,H}$ in the α = 1, ι = 1 simulation (Fig. 8c), where values are above 20 % in a large area at and perpendicular to the flow divide. In this simulation (α = 1, ι = 1), percentage differences between λ 2 and 
$\lambda_{2,H}$ (Figs. 9c) are high compared to other simulations, but the high values localised at and in the area perpendicular to the flow divide occupy a smaller area than the percentage differences between λ 1 and 
$\lambda_{1,H}$ (Fig. 8c). In the α = 0, ι = 0.6 and the α = 0.06, ι = 1 simulations, percentage differences between λ 1 and 
$\lambda_{1,H}$ are negligible across a large portion of the ice rise, in particular on the stoss side (Fig. 8b, d). The highest percentage differences are either side of the flow divide with values up to and over 20 %, reducing further away from the flow divide. Percentage differences between λ 2 and 
$\lambda_{2,H}$ in the α = 0, ι = 0.6 simulation are negligible except for small areas at the flow divide (Fig. 9b). As the assumption that the eigenvector corresponding with the largest eigenvalue is made for grounded ice, we do not analyse the differences between eigenvalues of the 
$2 \times 2$ and the 
$3 \times 3$ crystal orientation tensors in the ice shelves.
The percentage difference between 
$\lambda_{1,H}$ and λ 1 at an elevation of z = 0 for (a) the α = 0, ι = 1 simulation, (b) the α = 0, ι = 0.6 simulation, (c) the α = 1, ι = 1 simulation and (d) the α = 0.06, ι = 1 simulation. The solid line contours show depth below surface, and the dashed line is the grounding line.
The percentage difference between 
$\lambda_{2,H}$ and λ 2 at an elevation of z = 0 for (a) the α = 0, ι = 1 simulation, (b) the α = 0, ι = 0.6 simulation, (c) the α = 1, ι = 1 simulation and (d) the α = 0.06, ι = 1 simulation. The solid line contours show depth below surface, and the dashed line is the grounding line.
Next, we present results for the difference between the eigenvalues, 
$\lambda_{2,H} - \lambda_{1,H}$, of the 
$2 \times 2$ horizontal part of the 
$3 \times 3$ crystal orientation tensor. This metric corresponds directly with results from radar data where a vertical eigenvector assumption is made (Fig. 10). The α = 0, ι = 1 (Fig. 10a) and the α = 0.06, ι = 1 (Fig. 10d) simulations show similar results both spatially and in terms of the magnitude of differences, with differences staying below 0.1 at and in the vicinity of the flow divide. The largest differences in the eigenvalues occur on the eastern side of the ice rise, close to the grounding line as well as in small isolated areas elsewhere on the stoss side of the ice rise in the vicinity of the grounding line. Although eigenvalue differences differ substantially in magnitude between the two simulations, the α = 0, ι = 0.6 (Fig. 10b) and the α = 1, ι = 1 (Fig. 10c) simulations show similar spatial patterns. In the α = 0, ι = 0.6 simulation (Fig. 10b), the highest differences between 
$\lambda_{2,H}$ and 
$\lambda_{1,H}$ area seen at the tails of the flow divide, extending towards the grounding line north-east of the ice rise and in isolated areas on the east of the ice rise. Elsewhere, differences between 
$\lambda_{2,H}$ and 
$\lambda_{1,H}$ remain below 0.1. The α = 1, ι = 1 simulation (Fig. 10c) shows large differences between 
$\lambda_{2,H}$ and 
$\lambda_{1,H}$ of over 0.3 at the flow divide and extending to the grounding line, particularly on the lee side of the ice rise. As in all other simulations, differences of over 0.3 are seen on the eastern side of the ice rise close to the grounding line.
The difference between the two eigenvalues of the 
$2 \times 2$ horizontal crystal orientation tensor for (a) the α = 0, ι = 1 simulation, (b) the α = 0, ι = 0.6 simulation, (c) the α = 1, ι = 1 simulation and (d) the α = 0.06, ι = 1 simulation. The solid lines show contours for the depth below the surface in metres, and the dashed line is the grounding line.
Finally, we investigate the dependence of the fabric anisotropy field on the strain rate field and, in particular, the relationship between the largest horizontal fabric anisotropy direction and the largest horizontal strain rate direction. Given that we expect c-axes to point perpendicular to the direction of maximum stretching, we calculate the largest horizontal strain rate direction by calculating the eigenvalues and eigenvectors of the upper left 
$2 \times 2$ tensor of the 
$3 \times 3$ strain rate tensor, 
$\dot{\boldsymbol{\varepsilon}}$. We denote the eigenvectors of the horizontal strain rate tensor by 
$\bf{w}_{1,H}$ and 
$\boldsymbol{w}_{2,H}$, where 
$\boldsymbol{w}_{1,H} \leqslant \boldsymbol{w}_{2,H}$. We then calculate the angle between the maximum horizontal strain rate direction, 
$\boldsymbol{w}_{2,H}$, and the maximum horizontal fabric anisotropy direction, 
$\boldsymbol{v}_{2, H}$ (Fig. 11). In the α = 0, ι = 1 simulation (Fig. 11a) and the α = 0.06, ι = 1 simulation (Fig. 11d), the angle between 
$\boldsymbol{w}_{2,H}$ and 
$\boldsymbol{v}_{2, H}$ remains close to 
$\pi / 2$, meaning that the two vectors are mostly close to perpendicular. Exceptions are in small areas of divergence at the two tails of the flow divide. In the α = 0, ι = 0.6 simulation (Fig. 11b), the angle between 
$\boldsymbol{w}_{2,H}$ and 
$\boldsymbol{v}_{2, H}$ is mainly high everywhere on the ice rise except at the two tails of the flow divide where values are consistently closer to zero. Lastly, the α = 1, ι = 1 simulation (Fig. 11c) shows a large area on the stoss side of the ice rise, away from the flow divide, where 
$\boldsymbol{w}_{2,H}$ and 
$\boldsymbol{v}_{2, H}$ are mostly close to perpendicular. At the flow divide, there is a thin line where the two vectors are perpendicular and in the areas perpendicular to the flow divide, the angle between the two vectors is generally 
$ \gt \!\pi / 4$. This simulation shows the largest areas where the angle between 
$\boldsymbol{w}_{2,H}$ and 
$\boldsymbol{v}_{2, H}$ is close to zero from the extremes of the flow divide as far or almost as far as the grounding line.
The angle between the eigenvectors corresponding to the larger horizontal crystal orientation eigenvalue and the larger horizontal strain rate eigenvalue for (a) the α = 0, ι = 1 simulation, (b) the α = 0, ι = 0.6 simulation, (c) the α = 1, ι = 1 simulation and (d) the α = 0.06, ι = 1 simulation.
4. Discussion
We present the results of three-dimensional simulations of the full-tensor fabric anisotropy field of DIR via a coupling with the velocity and stress fields. The simulated domain includes the surrounding ice shelf, allowing analysis of the fabric anisotropy field across the grounding zone. Based on four previous studies (Martin and others, Reference Martín, Hindmarsh and Navarro2009; Seddik and others, Reference Seddik, Greve, Zwinger and Placidi2011; Martin and Gudmundsson, Reference Martín and Gudmundsson2012; Gagliardini and others, Reference Gagliardini2013), we simulate the fabric anisotropy field varying parameter choices controlling the relative influence of the strain rate and deviatoric stress tensors. By decomposition of the simulated crystal orientation tensor into eigenvalues and eigenvectors, our results provide a necessary step towards the comparison of observed ice fabric with three-dimensional modelled results for a variety of flow regimes.
4.1. Dependence of the fabric anisotropy field on the velocity and stress fields
The various flow regimes in the ice rise influence the fabric anisotropy field in different ways, with past flow regimes of a volume of ice also having an effect on the present fabric anisotropy, which is captured due to the semi-Lagrangian implementation of the crystal orientation tensor evolution in the simulations in this study.
Based on previous choices (Martin and others, Reference Martín, Hindmarsh and Navarro2009; Seddik and others, Reference Seddik, Greve, Zwinger and Placidi2011; Martin and Gudmundsson, Reference Martín and Gudmundsson2012; Gagliardini and others, Reference Gagliardini2013), we adjust the parameters α and ι in Eq. (7). The parameter α, taking values 
$\alpha \in [0, 1]$, controls the relative influence of the strain rate and deviatoric stress tensors, with a value of α = 0 meaning that the deviatoric stress tensor is ignored and a value of α = 1 meaning that the strain rate tensor is ignored. The parameter ι, taking values of 
$\iota \in [0, 1]$, controls the rate of change the crystal orientation tensor, 
$\bf{a}^{(2)}$, undergoes in response to the combination of the strain rate and deviatoric stress tensor. For example, for a value of ι close to zero, apart from rotational changes due to the spin tensor, the crystal orientation tensor changes minimally in response to the combination of the strain rate and deviatoric stress tensors. Because the strain rate and deviatoric stress tensor are co-dependent, differentiation between the effect of the various combinations of α and ι on the crystal orientation tensor is not trivial. Broadly speaking, although this is not always the case, the behaviour of the fabric anisotropy field in the α = 0, ι = 0.6 simulation shows similar results to the α = 1, ι = 1 simulation, despite the relatively large differences in parameter choice. Conversely, the α = 0, ι = 1 and the α = 0.06, ι = 1 simulations show only slightly differing results, likely explained by the similarity in values of α and ι.
While certain anisotropic fabric types can be expected under particular flow and stress regimes, the length of time the ice has undergone deformation in this regime also has an influence on the reflected fabric anisotropy field. In Figure 4, it can be seen that at the elevation z = 0, the ice has developed into a relatively strong single maximum fabric in the α = 0, ι = 1 simulation (Fig. 4a) and the α = 0.06, ι = 1 simulation (Fig. 4d) compared with the other two simulations. On the other hand, the α = 1, ι = 1 simulation (Fig. 4c) does not show a strong single maximum. This indicates that the strain rate term of Eq. (9) has a stronger influence on the development of a single maximum under vertical compression at a flow divide than the deviatoric stress tensor term of Eq. (9). Interestingly, in the α = 0, ι = 1 simulation (Fig. 4a) and the α = 0.06, ι = 1 simulation (Fig. 4d), the areas with a higher single maximum are located at the two tails of the flow divide. Although not explicitly simulated in this study, we expect that the fabric type in the area of high horizontal divergence at the tails of the flow divide would show a similar fabric to that forming at a dome flow divide, defined by a flow divide with horizontal flowlines originating at a point source. In the vicinity of the flow divide, the α = 1, ι = 1 simulation shows a vastly different horizontal fabric anisotropy direction than the other simulations (Fig. 5c). The eigenvector corresponding to the largest horizontal eigenvalue originates at a point source as opposed to being parallel to the flow divide in the other simulations. Values of 
$\ln (\lambda_2 / \lambda_1)$ are high along the flow divide and in relatively large areas where there is high horizontal divergence at the tails of the flow divide, meaning that the ice has a girdle fabric. Of note is that Figure 11c shows that in the areas of girdle fabric, the maximum horizontal strain rate direction is not perpendicular to the maximum horizontal fabric anisotropy direction. The reason for the differences between the α = 1, ι = 1 and the other three simulations is due to the largely different α value, which in this case results in ignoring the strain rate tensor completely.
Moving into the flanks of the ice rise, where vertical shear dominates, although the flow of ice in the grounded areas is horizontally divergent, the fabric types vary substantially. We identify a number of factors as having an influence on the ice fabric type: (a) whether flank flow is perpendicular to the flow divide or at the tails of the flow divide, (b) whether the ice is on the stoss or lee side of the ice rise and (c) how close ice is to the grounding line and what type of flow regime is active at the grounding zone. In flow regimes dominated by vertical shear, as is the case in the areas perpendicular to the flow divide, it is expected that fabric shows a single maximum with a slight offset to the vertical (Llorens and others, Reference Llorens2022). However, the α = 0, ι = 1 (Fig. 4a) and the α = 0.06, ι = 1 (Fig. 4d) simulations show the strongest single maxima perpendicular to the flow divide, but Figure 7a and d show that these simulations have the smallest tilt angle compared to the other simulations.
Across the grounding line, the predicted ln(k) values allow differentiation between single maximum fabric 
$(\ln (k) \gt 0)$ and girdle fabrics 
$(\ln (k) \lt 0$, Fig. 6). It must be noted, however, that a large positive or negative value of 
$\ln (k)$ does not necessarily imply that the ice has a strong single maximum or strong girdle fabric, respectively. Any analysis of the 
$\ln (k)$ field must be done in conjunction with plots of 
$\ln (\lambda_2 / \lambda_1)$ (Fig. 5) and 
$\ln (\lambda_3 / \lambda_2)$ (Fig. 4) (Woodcock, Reference Woodcock1977). In all simulations, (
$\ln (\lambda_3 / \lambda_2)$) is the lowest at the grounding zone and in the northern and eastern sides of the ice shelf, meaning that there is no single maximum in these areas. Contrary to this, the metric for a girdle fabric (
$\ln (\lambda_2 / \lambda_1)$) shows a small band at the grounding zone with values of 
$\ln (\lambda_2 / \lambda_1) = 0.5$ to 1. This weak girdle fabric may be due to the transition from a fabric with a vertical single maximum, to a single maximum aligned with the horizontal plane, via a transition across the grounding zone.
4.2. Implications for viscosity coupling
One ultimate goal in modelling the fabric anisotropy field of ice is to be able to fully couple the velocity and stress fields to the fabric anisotropy field and vice versa (note, we will refer to these simulations as fully coupled hereafter). In this study, we have explored the influence of the velocity and stress fields on the fabric anisotropy field using three-dimensional simulations, but have not investigated the coupling of the fabric anisotropy field to the velocity field via a direction-dependent viscosity. Here, we discuss areas of DIR where we expect the fabric type to influence viscosity as well as barriers to fully coupled, three-dimensional anisotropic simulations.
Generally speaking, ice crystal c-axes rotate in a direction which allows for optimal shearing, with studies showing that if ice is highly anisotropic, deformation in the plane perpendicular to the c-axis is significantly enhanced compared with deformation in other directions (Duval and others, Reference Duval, Ashby and Anderman1983). Areas where ice has a single maximum will experience ease of deformation parallel to basal planes and stiffening in response to deformation parallel to the c-axes. In the simulations of DIR, a single maximum fabric is most apparent in the α = 0, ι = 1 and the α = 0.06, ι = 1 simulations (Fig. 4). If enhancement factors were used (Gillet-Chaulet and others, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and Castelnau2005; Ma and others, Reference Ma, Gagliardini, Ritz, Gillet-Chaulet, Durand and Montagnat2010) in a coupling with the ice viscosity, this would allow more ease of vertical shear in these areas which is important in the evolution of Raymond arches at flow divides, as shown in two-dimensional, anisotropic simulations (Martin and others, Reference Martín, Hindmarsh and Navarro2009; Martin and Gudmundsson, Reference Martín and Gudmundsson2012; Drews and others, Reference Drews, Matsuoka, Martín, Callens, Bergeot and Pattyn2015).
A number of hurdles hinder the three-dimensional simulation of fully coupled ice anisotropy including computational expense, numerical instabilities, a lack of comparison with observations, parameter uncertainty and challenges associated with parallelisation. In a study by Gerber and others Reference Gerber(2023) in which flow tube simulations of the Elmer/Ice model of fully coupled ice anisotropy were performed, it was noted that numerical instabilities significantly hindered the simulation of areas with high velocities. Furthermore, computational expense and an increased memory load means that a directional-dependent viscosity creates a further hindrance. Despite this, simulations such as those presented in this paper allow ease of comparison with observations and should be included in comparisons between a hierarchy of models with varying complexity.
4.3. Model representation of anisotropy
Ice rises are an ideal location for testing models of ice anisotropy models as they contain various flow regimes which are relevant for other locations in the Antarctic Ice Sheet. Ice rises contain many of the same features as entire ice sheets, and results can inform expectations on studies in larger scales (Matsuoka and others, Reference Matsuoka2003). The shear zones on the western and eastern sides of DIR are representative of behaviour typical at lateral shear zones in ice shelves or ice streams (Smith and others, Reference Smith2017). Similarly, flow divides, flank flow and ice shelf flow are typical in many areas of Antarctica. More studies are needed to examine whether an anisotropic flow law is needed in order to make accurate projections of sea level rise. This is especially important as sudden speed up of ice, for example, in the initiation of an ice stream, is associated with anisotropic behaviour (Thorsteinsson and others, Reference Thorsteinsson, Waddington and Fletcher2003; Lilien and others, Reference Lilien, Rathmann, Hvidberg and Dahl-Jensen2021; Gerber and others, Reference Gerber2023). This necessitates further studies for comparison between various hierarchies of isotropic and anisotropic ice flow models with varying model complexity as well as comparison with observations. In our approach, fabric anisotropy is modelled using a crystal orientation tensor, which can describe many fabric types, but does not capture cone fabric shapes as modelled by Pettit and others Reference Pettit, Thorsteinsson, Jacobson and Waddington(2007) or multiple single maxima. An example of a representation of crystal orientation which can capture complex fabrics has been developed by Rathmann and others Reference Rathmann, Hvidberg, Grinsted, Lilien and Dahl-Jensen(2021). Furthermore, the introduction of additional, higher-order crystal orientation tensors in Eq. (7) would allow a more accurate representation of crystal orientation distribution, yet would introduce additional computational expensive.
In addition to difficulties in choosing a correct mathematical representation of an fabric anisotropy, Eq. (7) is challenging to solve numerically as it is an advection equation, with no diffusive term. We find that the use of semi-Lagrangian methods reduces dispersion, as they introduce diffusion through interpolation (Advani and Tucker, Reference Advani and Tucker1987). However, in areas with strong gradients in the velocity field, we nonetheless find that numerical dispersion deteriorates the solution with time. The alternative approach to this numerical problem has been to either explicitly incorporate numerical diffusion in Eq. (7) as in Seddik and others Reference Seddik, Greve, Zwinger and Placidi(2011) or to include re-crystallisation terms as in Gagliardini and others Reference Gagliardini(2013). We find it to be more effective to constrain the degree of orientation, but our approach can easily be adapted to include a re-crystallisation term in Eq. (7).
4.4. Framework for comparison with observations
A crucial component of accurate modelling of the fabric anisotropy field of ice is the comparison with observations, which, until recently, have been possible only by comparison with a limited number of ice cores. From ice core data, it is possible to construct a crystal orientation tensor which is directly comparable to simulations such as those presented in this work. Recently, advances in the measurement of anisotropic fabric using radar have allowed the acquisition of a significant amount of data with much greater ease than expensive ice coring projects (Young and others, Reference Young, Martín, Christoffersen, Schroeder, Tulaczyk and Dawson2021a; Reference Young2021b; Ershadi and others, Reference Ershadi2022; Reference Ershadi2025). However, assumptions made in the processing of data mean that the methods currently used are only valid where the dominant c-axis direction is exactly vertical (Rathmann and others, Reference Rathmann, Lilien, Grinsted, Gerber, Young and Dahl-Jensen2022).
Thus far, observational fabric anisotropy data have been used for inferring ice flow history including flow re-organisation (Durand and others, Reference Durand2007; Matsuoka and others, Reference Matsuoka, Power, Fujita and Raymond2012; Brisbourne and others, Reference Brisbourne, Martín, Smith, Baird, Kendall and Kingslake2019; Reference Brisbourne, Kendall, Kufner, Hudson and Smith2021) as well as for investigation of the state of the ice fabric in ice conditions such as ice streams (Smith and others, Reference Smith2017; Kufner and others, Reference Kufner2023), in layers of enhanced shear (Pettit and others, Reference Pettit2011) and at ice rises (Drews and others, Reference Drews, Matsuoka, Martín, Callens, Bergeot and Pattyn2015; Ershadi and others, Reference Ershadi2025). In terms of making comparisons between observations and anisotropic models, matching of isochrones and velocities has been performed (Drews and others, Reference Drews, Matsuoka, Martín, Callens, Bergeot and Pattyn2015; McCormack and others, Reference McCormack, Warner, Seroussi, Dow, Roberts and Treverrow2022), but few studies have compared directly the observed and modelled fabric anisotropy fields (Lilien and others, Reference Lilien2023).
Although we do not explicitly make comparisons between the observed and modelled fabric anisotropy fields of DIR, we have described steps necessary for such a comparison. Our results show very differing three-dimensional fabric anisotropy fields depending on the chosen influence of the strain rate and deviatoric stress tensors in Eq. (7) used in previous studies Martin and others (Reference Martín, Hindmarsh and Navarro2009); Seddik and others (Reference Seddik, Greve, Zwinger and Placidi2011); Martin and Gudmundsson (Reference Martín and Gudmundsson2012); Gagliardini and others (Reference Gagliardini2013) and draw attention to the need to constrain equation parameters. We have provided a method to compare the difference between the two horizontal eigenvalues (Fig. 10), sometimes referred to as the horizontal anisotropy, which is directly comparable to radar observations. We recommend radar surveys in the areas of high horizontal flow divergence at the tails of flow divides, as these are areas where differing values of α and ι show significantly differing fabric anisotropy fields.
5. Conclusions
In this study, we have modelled the three-dimensional, diagnostic, fabric anisotropy field of DIR including domain partitioning using a crystal orientation tensor evolution equation describing the spatial distribution of ice crystal c-axes in a given ice volume. Simulations are performed varying the relative influence of the strain rate and deviatoric stress tensors in the crystal orientation tensor evolution equation as used in four previous studies (Martin and others, Reference Martín, Hindmarsh and Navarro2009; Seddik and others, Reference Seddik, Greve, Zwinger and Placidi2011; Martin and Gudmundsson, Reference Martín and Gudmundsson2012; Gagliardini and others, Reference Gagliardini2013), with results showing significantly differing ice fabrics across the flow divide and the flanks of DIR. Lastly, we provide a modelling framework for comparison with radar observations, outlining areas where the assumption of one vertical eigenvector may not hold and resulting errors in horizontal eigenvalues.
Data availability statement
The code needed to perform the simulations and create the presented in the figures here are available in a publicly accessible, permanent repository (Henry, Reference Henry2024d). The accompanying model input and output data as well as the processed data are also available in a publicly accessible, permanent repository (Henry, Reference Henry2024c). The simulations presented here are dependent on the results of simulations presented in a previous study (Henry and others, Reference Henry, Schannwell, Višnjević, Millstein, Bons and Eisen2023), for which code and data are available (Henry, Reference Henry2024a, Reference Henryb).
Acknowledgements
C.H. was supported by the Deutsche Forschungsgemeinschaft in the framework of the priority programme 1158 ‘Antarctic Research with comparative investigations in Arctic ice areas’ by a grant SCHA 2139/1-1. R.D. was partially supported by an Emmy Noether Grant of the Deutsche Forschungsgemeinschaft (DR 822/3-1). This work used resources of the Deutsches Klimarechenzentrum granted by its Scientific Steering Committee. We would like to thank the reviewers and editor for comments which helped to improve the clarity of the text. We would furthermore like to thank Luca Schmidt for an internal review of the manuscript.
Author contribution
A.C.J.H. and C.M. designed the study, A.C.J.H. performed the simulations with input from C.M., A.C.J.H. performed the analysis with input from C.M. and A.C.J.H. wrote the manuscript which was edited by C.M. and R.D.
Competing interests
We have no competing interests to declare.
Appendix
The representation of the fabric anisotropy field using a crystal orientation tensor,
\begin{equation}
\mathbf{a}^{(2)} = \begin{pmatrix}
a^{(2)}_{xx} & a^{(2)}_{xy} & a^{(2)}_{xz} \\
a^{(2)}_{yx} & a^{(2)}_{yy} & a^{(2)}_{yz} \\
a^{(2)}_{zx} & a^{(2)}_{zy} & a^{(2)}_{zz} \end{pmatrix},
\end{equation}describing the c-axis distribution is ideal for comparison with ice core or radar observations in that an equivalent crystal orientation tensor can be constructed from observational data by determining the cross-products of cosines of each c-axis and summing up such that
\begin{equation}
\mathbf{a}^{(2)} = \begin{pmatrix}
a^{(2)}_{xx} & a^{(2)}_{xy} & a^{(2)}_{xz} \\
a^{(2)}_{yx} & a^{(2)}_{yy} & a^{(2)}_{yz} \\
a^{(2)}_{zx} & a^{(2)}_{zy} & a^{(2)}_{zz} \\
\end{pmatrix}
= \frac{1}{N} \begin{pmatrix}
\sum l_i^2 & \sum l_i m_i & \sum l_i n_i \\
\sum m_i l_i & \sum m_i^2 & \sum m_i n_i \\
\sum n_i l_i & \sum n_i m_i & \sum n_i^2 \\
\end{pmatrix},
\end{equation} where N is the number of c-axes being summed over. The cosines of the c-axis vectors are defined by taking the cosine of the angle between the c-axis vector, ci, and each positive coordinate axis, 
$\boldsymbol{\hat{e}}_x$, 
$\boldsymbol{\hat{e}}_y$ and 
$\boldsymbol{\hat{e}}_z$, so that
\begin{equation}
\begin{aligned}
l_i & = c_{i,x}, \\
m_i & = c_{i,y}, \\
n_i & = c_{i,z}.
\end{aligned}
\end{equation}From the orientation tensor, a number of metrics can be calculated to investigate the degree of crystal orientation, the predominant c-axes directions and the fabric type.
Open access
