## Modelling the three-dimensional, diagnostic anisotropy field of an ice rise

A. Clara J. Henry 1, 2, * , Carlos Mart´ ın , and Reinhard Drews 3 2

1 Max Planck Institute for Meteorology, Hamburg, Germany 2 Department of Geosciences, University of T¨bingen, T¨bingen, u u Germany

3 British Antarctic Survey, Natural Environmental Research Council, Cambridge, UK

* Now at: Department of Mathematics, Stockholm University, Sweden

## Abstract

Polar ice develops anisotropic crystal orientation fabrics under deformation, yet ice is most often modelled as an isotropic fluid. We present threedimensional simulations of the crystal orientation fabric of Derwael Ice Rise including the surrounding ice shelf using a crystal orientation tensor evolution equation corresponding to a fixed velocity field. We use a semiLagrangian numerical method that constrains the degree of crystal orientation evolution to solve the equations in complex flow areas. We perform four simulations based on previous studies, altering the rate of evolution of the crystal anisotropy and its dependence on a combination of the strain rate and deviatoric stress tensors. We provide a framework for comparison with radar observations of the anisotropy field, outlining areas where the assumption of one vertical eigenvector may not hold and provide resulting errors in measured eigenvalues. We recognise the areas of high horizontal divergence at the ends of the flow divide as important areas to make comparisons with observations. Here, poorly constrained model parameters result in the largest difference in fabric type. These results are important in the planning of future campaigns for gathering data to constrain model parameters and as a link between observations and computationally-efficient, simplified models of anisotropy.

## 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, 1988). The orientation of

Correspondence: clara.henry@math.su.se

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 anisotropy data and ice flow models.

Ice cores have long been used to quantify the the direction of ice crystals and the degree of anisotropy with depth (Durand et al., 2009; Montagnat et al., 2014; Weikusat et al., 2017). However, the lack of horizontal spatial information hinders the study of the relationship between ice flow and anisotropy. Apart from ice cores sites, anisotropic information is sparse and relies on geophysical methods such as seismics, both active (Diez et al., 2014; Brisbourne et al., 2019) and passive (Smith et al., 2017; Kufner et al., 2023) or radar (Fujita et al., 2006; Matsuoka et al., 2012; Drews et al., 2012). The main limitation is that, generally, seismic methods only provide depthaveraged anisotropy and radar only provide 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 anisotropy have resulted in an increase in the acquisition of observational fabric data (Dall, 2010; Young et al., 2021a,b; Jordan et al., 2022; Ershadi et al., 2022). These surveys provide the spatial variability of fabric in more extensive areas and aid understanding of the relation between ice 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 anisotropy field to the velocity and stress fields in the presented 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 et al., 2023). Several modelling studies investigating the effects of an anisotropic rheology on ice flow dynamics at ice divides in two dimensions (Mart´ ın et al., 2009; Mart´ ın and Gudmundsson, 2012; Lilien et al., 2023) and in ice streams (Gerber et al., 2023; Richards et al., 2023) have shown the importance of including anisotropy evolution in ice flow models. In idealised, two-dimensional simulations, Richards et al. (2022) investigate the influence of the strain-rate and spin on the fabric type. Building on the work of Thorsteinsson et al. (2003), Gillet-Chaulet et al. (2006), and Pettit et al. (2007), Rathmann et al. (2021) investigated the coupling between 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 anisotropy coupled to the viscosity via enhancement factors, McCormack et al. (2022) studied the effect of anisotropy on larger-scale three-dimensional models.

The coupling of anisotropy evolution in ice flow models is not without complication due to challenges with numerical stability, the parameter choice within the anisotropy evolution equation, and the coupling with viscosity in the Glen's flow law equation. 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 anisotropy is needed.

In this paper, we present the first model of the three-dimensional anisotropy field of an ice rise using the finite element model Elmer/Ice (Gagliardini et al., 2013) across 100 partitions and applied to Derwael Ice Rise. We investigate the dependence of the anisotropy field on the strain rate and deviatoric stress fields without coupling with the viscosity and without recrystallisation 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 verticalshear-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 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, 1988), knowledge of how anisotropic polar ice is (Duval et al., 1983) and field evidence of the effect crystal orientation fabric has on ice flow (Gerber et al., 2023). After the introduction of the crystal orientation tensor (G¨dert, 2003), the flow-induced o crystal orientation fabric evolution can be considered by large-scale ice-sheet models (Gagliardini et al., 2013).

We believe that the main reasons for the lack of progress are: (a) the numerical implementation of crystal orientation fabric evolution is challenging (Seddik et al., 2011), (b) the interpretation of model output and comparison with non-comprehensive observations is complex (Jordan et al., 2022), and (c) essential model parameters within the theory are not yet constrained by observations (Ma et al., 2010).

Numerically, the main issue is that numerical dispersion tends to break down the orientation tensor properties. Here, we present a numericallyrobust, 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 Derwael Ice Rise (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 con-

Figure 1: In (a), an overview is shown of Derwael Ice Rise, with the surrounding ice shelves and the grounding line (Jezek et al., 2013; Rignot and Scheuchl., 2017) The velocity field is shown in colour. In (b), the simulated upper surface velocity field of Derwael Ice Rise is shown, based on simulations in Henry et al. (2023). The boxes A and B show the areas referred to as the areas of high horizontal divergence at the tails of the flow divide.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-c2f829cd-f91e-4e68-a7ed-2ed06d3d12e7.jpg)

tinental grounding line position (Favier and Pattyn, 2015; Schannwell et al., 2019, 2020; Henry et al., 2022). Formation, evolution and disintegration of ice rises occur over glacial-interglacial timescales, meaning that remnants of ice properties such as temperature and 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 et al., 2015). 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 Derwael Ice Rise (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 et al., 2023). This velocity field is based on a transient, thermomechanically-coupled full Stokes model with an isotropic rheology. The model was forced with an observationallyconstrained surface mass balance field and predictions of the model output were validated with extensive radar observations which are available for this Ice Rise (Koch et al., 2023). Second, previous studies suggest that DIR is likely close to steady state (Callens et al., 2016) possibly with a minor amount of thinning (Drews et al., 2015). 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 .

## 2 Methods

The simulations use output from Henry et al. (2023) as a starting point in which the finite element model Elmer/Ice (Gagliardini et al., 2013) was applied to DIR. Most importantly, we use the predicted three-dimensional, steady-state velocity field to make predictions of the anisotropic ice properties 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

$$\nabla \cdot ( \tau - P \mathbf I ) + \rho _ { i } \mathbf g = 0,$$

where τ is the deviatoric stress tensor, P is the pressure, I is the identity matrix, ρ i is the ice density and g = g ˆ e z is the gravitational acceleration. The ice is subject to an incompressibility condition,

$$\nabla \cdot \mathbf u = 0,$$

and the Glen's flow law,

$$\tau = 2 \eta \dot { \varepsilon },$$

which describes the nonlinear dependence between the strain rate tensor, ˙, ε and the deviatoric stress tensor. The effective viscosity, η , is

$$\eta = \frac { 1 } { 2 } E A ( T ^ { \prime } ) ^ { - 1 / n } \dot { \varepsilon } _ { e } ^ { ( 1 - n ) / n }.$$

where E is an enhancement factor which is spatially and temporally constant here, 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 et al. (2023). Here, n is the Glen's flow law exponent, and ˙ ε e = √ tr ˙ ε 2 / 2 is the effective strain rate.

Keeping the velocity and stress fields constant in time, a semi-Lagrangian anisotropy evolution equation is coupled and simulations are performed for 20000 years. To initialise the anisotropy simulation, the simulation named n3E0.5dsdt50 in Henry et al. (2023) is used. The simulation time of 20000 years was deemed appropriate given that in Henry et al. (2023), it was predicted that the ice in DIR is roughly 8000 years old at a depth of 95 %.

The anisotropy of ice is described by the second and fourth order orientation tensors, a (2) and a (4) , respectively defined by

$$a _ { i j } ^ { ( 2 ) } = \langle c _ { i } c _ { j } \rangle$$

and

$$a _ { i j k l } ^ { ( 4 ) } = \langle c _ { i } c _ { j } c _ { k } c _ { l } \rangle,$$

where c is the c -axis unit vector and the operator ⟨⟩ denotes the average over all the grains that compose the ice polycrystal.

The representation of the anisotropy field using a crystal orientation tensor,

$$\mathbf a ^ { ( 2 ) } = \begin{pmatrix} a _ { x x } ^ { ( 2 ) } & a _ { x y } ^ { ( 2 ) } & a _ { x z } ^ { ( 2 ) } \\ a _ { y x } ^ { ( 2 ) } & a _ { y y } ^ { ( 2 ) } & a _ { y z } ^ { ( 2 ) } \\ a _ { z x } ^ { ( 2 ) } & a _ { z y } ^ { ( 2 ) } & a _ { z z } ^ { ( 2 ) } \end{pmatrix},$$

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

$$\mathbf a ^ { ( 2 ) } = \begin{pmatrix} a _ { x x } ^ { ( 2 ) } & a _ { x y } ^ { ( 2 ) } & a _ { x z } ^ { ( 2 ) } \\ a _ { y x } ^ { ( 2 ) } & a _ { y y } ^ { ( 2 ) } & a _ { y z } ^ { ( 2 ) } \\ a _ { z x } ^ { ( 2 ) } & a _ { z y } ^ { ( 2 ) } & a _ { z z } ^ { ( 2 ) } \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}.,$$

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, c i , and each positive coordinate axis, ˆ e x , ˆ e y and ˆ e z , so that

$$\Psi ^ { ( 1 ) } _ { i } = \frac { c _ { i } \cdot \hat { e } _ { x } } { | c _ { i } | } \\ m _ { i } = \frac { c _ { i } \cdot \hat { e } _ { y } } { | c _ { i } | } \\ n _ { i } = \frac { c _ { i } \cdot \hat { e } _ { x } } { | c _ { i } | }. \\ \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { \Phi } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } & \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 2 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { (( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 3 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { - ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1) } _ { i + 1 } = \Phi ^ { \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1)) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } | \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1,) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } { \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 2 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ { ( 1 ) } _ { i + 1 } = \Phi ^ {$$

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. The eigenvalues, namely λ 1 , λ 2 and λ 3 , and the corresponding eigenvectors, v 1 , v 2 and v 3 , of the crystal orientation tensor a (2) satisfy

$$\mathbf a ^ { ( 2 ) } \mathbf v = \lambda \mathbf v$$

defining each eigenvalue by its subscript such that λ 1 ≤ λ 2 ≤ λ 3 . In the case of randomly oriented c -axes or isotropic ice, all three eigenvalues are similar in size such that λ 1 ≈ λ 2 ≈ λ 3 ≈ 1 3. / In areas where λ 3 ≥ λ 1 ≈ λ 2 , ice fabric is said to have a single maximum fabric with the majority of c -axes pointing in a single direction. Where λ 2 ≈ λ 3 ≥ λ 2 , ice is said to have a girdle fabric, with c -axes following an arc or circle.

The evolution of the second-order orientation tensor can be described by

$$\left ( \frac { \partial } { \partial t } + \mathfrak { u } \cdot \mathbf \nabla \right ) \mathbf a ^ { ( 2 ) } & = \mathbf W \mathbf a ^ { ( 2 ) } - \mathbf a ^ { ( 2 ) } \mathbf W \\ & \quad - \iota ( \mathbf C \mathbf a ^ { ( 2 ) } + \mathbf a ^ { ( 2 ) } \mathbf C - 2 \mathbf a ^ { ( 4 ) } \colon \mathbf C ),$$

where the tensor C is defined by

$$\mathbf C = ( 1 - \alpha ) \dot { \varepsilon } + \alpha \frac { 1 } { 2 \eta _ { 0 } } \tau,$$

and α is the so-called interaction parameter, which describes the relative influence of the strain rate and stress tensors. The parameter ι determines the rate at which the crystal orientation tensor is influences by the weighted combination of strain-rate and deviatoric stress tensors. The last term, α 1 2 η 0 τ , in Eq. (10) is represented as αk Aτ s n -1 e τ in Gagliardini et al. (2013), but the simulation calculations are equivalent. The spin of the ice is denoted by the spin tensor, W , and

$$\eta _ { 0 } = \frac { 1 } { 2 } ( E A ) ^ { - 1 / n } \left ( \frac { 1 } { 2 } \text{tr} ( \dot { \varepsilon } ^ { 2 } ) ^ { \frac { 1 - n } { 2 n } } \right ).$$

Early work recognised the influence of the cumulative strain and stress on ice crystal c -axis orientation (Alley, 1988), which led to the development of a crystal orientation tensor evolution equation (Eq. (9)) dependent on the velocity gradient through the spin and strain rate tensors, and the deviatoric stress tensor (Gillet-Chaulet et al., 2006). Although it is known that the velocity and stress fields have an influence on the anisotropy field of ice, it remains unclear what the relative influence is. The anisotropy evolution equation (Eq. (9)) 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, ˙ ε , 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 et al., 2009). The choice of α and ι in previous studies have 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, a (4) on the 3 × 3 crystal orientation tensor, a (2) . In reality, a (2) is dependent on further, higher-order tensors, but it has been shown that this dependence is negligible (Chung and Kwon, 2002).

To solve Eq. (9) we also require a relation between a (2) and a (4) , a closure approximation, and we use an invariant-based optimal fitting closure approximation (Gillet-Chaulet et al., 2006). The distribution of c -axes can be described by tensors of ever-increasing order, but a compromise is found by using the closure approximation in order to reduce computation time. The effect of the choice of α and ι on the anisotropy field is investigated using combinations of parameters from previous literature presented in Table 1.

The value of α determines the relative influence of the strain rate and deviatoric stress tensors on 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 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.

Table 1: Parameter combinations

|     | α      | ι     | Source                          |
|-----|--------|-------|---------------------------------|
| (a) | 0      | 1     | Mart´n ı et al. (2009)          |
| (b) | 0      | 0 . 6 | Seddik et al. (2011)            |
| (c) | 1      | 1     | Mart´n ı and Gudmundsson (2012) |
| (d) | 0 . 06 | 1     | Gagliardini et al. (2013)       |

Initially, all ice in the model domain is isotropic, described by the tensor a (2) = 1 3 I , where I is the 3 × 3 identity matrix. During transient simulation of the 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 anisotropic evolution model described in Eq. (9) together with its boundary and initial conditions is solved using a semi-Lagrangian method as described in Mart´ ın et al. (2009). The determinant of a (2) gives information about the degree of crystal orientation of the fabric (Advani and Tucker, 1987). We define here the degree of crystal orientation with the scalar value,

$$d = 1 - 3 ^ { 3 } \det \left ( \mathbf a ^ { ( 2 ) } \right ).$$

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 (9) with the condition,

$$\left ( \frac { \partial } { \partial t } + \mathfrak { u } \cdot \nabla \right ) d \geq 0,$$

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

Figure 2: The eigenvalues, λ 1 , λ 2 and λ 3 , of the crystal orientation tensor in the α = 0, ι = 1 simulation at an elevation of z = 0, corresponding to 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 Fig. 3 is taken.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-0271e074-fe67-4934-853a-1b1537a4879d.jpg)

Figure 3: Cross-sections through the flow divide as shown in the Fig. 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-01df0e07-b5d6-4647-9f9e-20e707052052.jpg)

Figure 4: The ratio of the two larger eigenvalues of the 3 × 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 anisotropy direction. The dashed line shows the grounding line and the solid line show contours of the depth below the upper ice surface.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-d1c636f0-88fd-4d04-9e42-eeb768165e9a.jpg)

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.

## 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. (9) 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 Derwael Ice Rise. 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, 1977), ln( λ /λ 3 2 ), is large, then the c -axes are more concentrated in a single direction. The direction of the largest horizontal crystal orientation eigenvector, 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( λ /λ 3 2 ) than predicted by the other two simulations. The simulation results differ, however, in the spatial variation of ln( λ /λ 3 2 ). The α = 0, ι = 1 simulation shows larger values of ln( λ /λ 3 2 ) at the tails of the ice rise flow divide (Boxes A and B in Fig. 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 two smaller eigenvalues, λ 1 and λ 2 are investigated using the metric ln( λ /λ 2 1 ) (Fig. 5). Large ratios between the two smaller eigenvalues indicates that c -axis directions for a given volume of ice are concentrated along an arc. In all simulations, there are high values of ln( λ /λ 2 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 show 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.

Figure 5: The ratio of the two smaller eigenvalues of the 3 × 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 anisotropy direction. The dashed line shows the grounding line and the solid line show contours of the depth below the upper ice surface.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-82753d3b-e5ae-4100-b0d3-e25f86a1b977.jpg)

Figure 6: 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 anisotropy direction. The dashed line shows the grounding line and the solid line show contours of the depth below the upper ice surface.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-968e8423-d0b9-49e7-8229-70965be50fbb.jpg)

The Woodcock k value (Woodcock, 1977), defined as

$$k = \frac { \ln ( \lambda _ { 3 } / \lambda _ { 2 } ) } { \ln ( \lambda _ { 2 } / \lambda _ { 1 } ) },$$

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 = ∞ . 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( λ /λ 3 2 ) &gt;&gt; 0, then the ice has a strong single maximum. If, on the other hand, ln( k ) &lt;&lt; 0 and ln( λ /λ 2 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

Figure 7: The angle between the eigenvector corresponding to the largest eigenvalue and the vertical direction at an elevation of z = 0 corresponding to 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-5c96a0cf-dd17-4849-9fc0-4d60a62e4bd3.jpg)

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.

## 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 (Jordan et al., 2022; Ershadi et al., 2022). Obliquely oriented fabric types can, in theory, also be detected (Matsuoka et al., 2009; Rathmann et al., 2022), 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 (Figs. 7a and 7d, respectively). In the α = 0, ι = 0 6 simulation (Fig. . 7b), differences between the z -direction and v 3 are

Figure 8: The percentage difference between λ 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-0fe2b6d4-4e33-4d22-ba03-45276ca4b619.jpg)

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 π/ 8 to π/ 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.

We furthermore calculate the eigenvectors and eigenvalues of the horizontal crystal orientation tensor, i.e. the upper left 2 × 2 part of the 3 × 3 tensor a (2) . The reason for this is that if the 3 × 3 crystal orientation tensor has one strictly vertical eigenvector, the eigenvalues and eigenvectors of the 2 × 2 tensor correspond to the horizontal eigenvalues and eigenvectors of the 3 × 3 tensor. We denote the eigenvalues of the 2 × 2 tensor by λ 1 ,H and λ 2 ,H , and the corresponding eigenvectors by v 1 ,H and v 2 ,H , respectively. By comparing the eigenvalues of the 2 × 2, horizontal crystal orientation tensor with the eigenvalues of the 3 × 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 λ 1 ,H (Fig. 8), as well as λ 2 and λ 2 ,H (Fig. 9). We find that if a difference exists between λ 1 and λ 1 ,H , or λ 2 and λ 2 ,H , then the 2 × 2 tensor eigenvalues always underestimate the smaller two 3 × 3 tensor eigenvalues. Furthermore, percentage differences tend to be higher for the smaller eigenvalues λ 1 and λ 1 ,H than for λ 2 and

λ 2 ,H , but exceptions to this are seen. The α = 0, ι = 1 simulation shows a negligible percentage difference between λ 1 and λ 1 ,H (Fig. 8a) and small percentage differences for λ 2 and λ 2 ,H of up to 10 % (Fig. 9a). A similar spatial distribution of percentage differences appears for λ 2 and λ 2 ,H in the α = 0 06, . ι = 1 simulation (Fig. 9d). Most notable is the percentage difference between λ 1 and λ 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 λ 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 λ 1 ,H (Fig. 8c). In the α = 0, ι = 0 6 and the . α = 0 06, . ι = 1 simulations, percentage differences between λ 1 and λ 1 ,H are negligible across a large portion of the ice rise, in particular on the stoss side (Figs. 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 λ 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 × 2 and the 3 × 3 crystal orientation tensors in the ice shelves.

Next, we present results for the difference between the eigenvalues, λ 2 ,H -λ 1 ,H , of the 2 × 2 horizontal part of the 3 × 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 λ 2 ,H and λ 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 λ 2 ,H and λ 1 ,H remain below 0 1. . The α = 1, ι = 1 simulation (Fig. 10c) shows large differences between λ 2 ,H and λ 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.

Finally, we investigate the dependence of the anisotropy field on the strain rate field, and in particular the relationship between the largest horizontal 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 × 2 tensor

Figure 9: The percentage difference between λ 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-bdf11f62-e5a3-4460-9f64-a2d0d92d624b.jpg)

of the 3 × 3 strain rate tensor, ˙. ε We denote the eigenvectors of the horizontal strain rate tensor by w 1 ,H and w 2 ,H , where w 1 ,H ≤ w 2 ,H . We then calculate the angle between the maximum horizontal strain rate direction, w 2 ,H , and the maximum horizontal anisotropy direction, v 2 ,H (Fig. 11). In the α = 0, ι = 1 simulation (Fig. 11a) and the α = 0 06, . ι = 1 simulation (Fig. 11d), the angle between w 2 ,H and v 2 ,H remains close to π/ 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 w 2 ,H and 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 w 2 ,H and 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; π/ 4. This simulation shows the largest areas where the angle between w 2 ,H and v 2 ,H is close to zero from the extremes of the flow divide as far, or almost as far as the grounding line.

## 4 Discussion

We present the results of three-dimensional simulations of the full-tensor anisotropy field of Derwael Ice Rise via a coupling with the velocity and

Figure 10: The difference between the two eigenvalues of the 2 × 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-453bf0d1-9371-4311-b85f-bc856c0b622e.jpg)

stress fields. The simulated domain includes the surrounding ice shelf, allowing analysis of the anisotropy field across the grounding zone. Based on four previous studies (Mart´ ın et al., 2009; Seddik et al., 2011; Mart´ ın and Gudmundsson, 2012; Gagliardini et al., 2013), we simulate the 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 threedimensional modelled results for a variety of flow regimes.

## 4.1 Dependence of the anisotropy field on the velocity and stress fields

The various flow regimes in the ice rise influence the anisotropy field in different ways, with past flow regimes of a volume of ice also having an effect on the present 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 (Mart´ ın et al., 2009; Seddik et al., 2011; Mart´ ın and Gudmundsson, 2012; Gagliardini et al., 2013), we adjust the parameters α and ι in Eq. (9). The parameter α , taking values α ∈ [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 ι ,

Figure 11: 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01069v1-figure-c1c7d9b1-9580-46a0-9db7-de8bdd1480ab.jpg)

taking values of ι ∈ [0 , 1], controls the rate of change the crystal orientation tensor, 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 via Glen's flow law (Eq. (3)), differentiation between the various combinations of α and ι is not trivial. Broadly speaking, though this is not always the case, the behaviour of the 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 anisotropy field. In Fig. 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. (10) 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. (10). 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 anisotropy direction than the other simulations (Fig. 5c). The eigenvector corresponding to the largest horizontal eigenvalue originate at a point source as opposed to being parallel to the flow divide in the other simulations. Values of ln( λ /λ 2 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 Fig. 11c shows that in the areas of girdle fabric, the maximum horizontal strain rate direction is not perpendicular to the maximum horizontal 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 et al., 2022). 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 Figs. 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( λ /λ 2 1 ) (Fig. 5) and ln( λ /λ 3 2 ) (Fig. 4) (Woodcock, 1977). In all simulations, (ln( λ /λ 3 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( λ /λ 2 1 )) shows a small band at the grounding zone with values of ln( λ /λ 2 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 anisotropy field of ice is to be able to fully couple the velocity and stress fields to the 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 anisotropy field using three-dimensional simulations, but have not investigated the coupling of the anisotropy field to the velocity field via a direction-dependent viscosity. Here, we discuss areas of Derwael Ice Rise 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 et al., 1983). 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 Derwael Ice Rise, 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 et al., 2005; Ma et al., 2010) 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 (Mart´ ın et al., 2009; Mart´ ın and Gudmundsson, 2012; Drews et al., 2015).

A number of hurdles hinder the three-dimensional simulation of a fullycoupled anisotropy field including computational expense, numerical instabilities, a lack of comparison with observations, parameter uncertainty and challenges associated with parallelisation. In a study by Gerber et al. (2023) in which flow tube simulations the Elmer/Ice model of a fully-coupled anisotropy field 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 direction-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 anisotropy 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 et al., 2003). The shear zones on the western and eastern sides of Derwael Ice Rise are representative of behaviour typical at lateral shear zones in ice shelves or ice streams (Smith et al., 2017). 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 et al., 2003; Lilien et al., 2021; Gerber et al., 2023). 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, 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 et al. (2007) or multiple single maxima. An example of a representation of crystal orientation which can capture more complex fabrics has been developed by Rathmann et al. (2021).

In addition to difficulties in choosing a correct mathematical representation of an anisotropic fabric, Eq. (9) is challenging to solve numerically as it is an advection equation, with no diffusive term. We find that the use of semi-Lagrangian methods reduce dispersion, as they introduce diffusion through interpolation (Advani and Tucker, 1987). 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. (9) as in Seddik et al. (2011), or to include re-crystallisation terms as in Gagliardini et al. (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. (9).

## 4.4 Framework for comparison with observations

A crucial component of accurate modelling of the 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 has allowed the acquisition of a significant amount of data with much greater ease than expensive ice coring projects (Young et al., 2021a,b; Ershadi et al., 2022). 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 et al., 2022).

Thus far, observational anisotropy data has been used for inferring ice flow history including flow re-organisation (Durand et al., 2007; Matsuoka et al., 2012; Brisbourne et al., 2019, 2021) as well as for investigation of the state of the ice fabric in ice conditions such as ice streams (Smith et al., 2017; Kufner et al., 2023), in layers of enhanced shear (Pettit et al., 2011) and at ice rises (Drews et al., 2015). In terms of making comparisons between observations and anisotropic models, matching of isochrones and velocities has been performed (Drews et al., 2015; McCormack et al., 2022), but few studies have compared directly the observed and modelled anisotropy fields (Lilien et al., 2023).

Although we do not explicitly make comparisons between the observed

and modelled anisotropy fields of Derwael Ice Rise, we have described steps necessary for such a comparison. Our results show very differing threedimensional anisotropy fields depending on the chosen influence of the strain rate and deviatoric stress tensors in Eq. (9) used in previous studies (Mart´ ın et al., 2009; Seddik et al., 2011; Mart´ ın and Gudmundsson, 2012; Gagliardini et al., 2013) 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 anisotropy fields.

## 5 Conclusions

In this study, we have modelled the three-dimensional anisotropy field of Derwael Ice Rise (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 (Mart´ ın et al., 2009; Seddik et al., 2011; Mart´ ın and Gudmundsson, 2012; Gagliardini et al., 2013), 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.

## 6 Acknowledgements

C. Henry was supported by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the priority programme 1158 'Antarctic Research with comparative investigations in Arctic ice areas' by a grant SCHA 2139/11. R. Drews was partially supported by an Emmy Noether Grant of the Deutsche Forschungsgemeinschaft (DR 822/3-1). This work used resources of the Deutsches Klimarechenzentrum (DKRZ) granted by its Scientific Steering Committee (WLA). We would like to thank Luca Schmidt for comments which improved the manuscript.

## 7 Author contribution

ACJH and CM designed the study, ACJH performed the simulations with input from CM, ACJH performed the analysis with input from CM, ACJH wrote the manuscript which was edited by CM and RD.

## 8 Code and data availability

The code needed to perform the simulations and create the presented in the figures here are available in a publicly-accessible, permanent repository (Henry, 2024a). The accompanying model input and output data, as well as the processed data are also available in a publicly-accessible, permanent repository (Henry, 2024c). The simulations presented here are dependent on the results of simulations presented in a previous study (Henry et al., 2023), for which code and data are available (Henry, 2024b,d).

## References

- S. G. Advani and I. Tucker, Charles L. The Use of Tensors to Describe and Predict Fiber Orientation in Short Fiber Composites. Journal of Rheology , 31(8):751-784, 11 1987. doi: 10.1122/1.549945.
- R. B. Alley. Fabrics in polar ice sheets: development and prediction. Science , 240(4851):493-495, 1988.
- A. M. Brisbourne, C. Mart´ ın, A. M. Smith, A. F. Baird, J. M. Kendall, and J. Kingslake. Constraining Recent Ice Flow History at Korff Ice Rise, West Antarctica, Using Radar and Seismic Measurements of Ice Fabric. Journal of Geophysical Research: Earth Surface , 124(1):175-194, 2019. doi: 10.1029/2018JF004776.
- A. M. Brisbourne, M. Kendall, S.-K. Kufner, T. S. Hudson, and A. M. Smith. Downhole distributed acoustic seismic profiling at skytrain ice rise, west antarctica. The Cryosphere , 15(7):3443-3458, 2021. doi: 10. 5194/tc-15-3443-2021.
- D. Callens, R. Drews, E. Witrant, M. Philippe, and F. Pattyn. Temporally stable surface mass balance asymmetry across an ice rise derived from radar internal reflection horizons through inverse modeling. Journal of Glaciology , 62(233):525-534, 2016. doi: 10.1017/jog.2016.41.
- D. H. Chung and T. H. Kwon. Invariant-based optimal fitting closure approximation for the numerical prediction of flow-induced fiber orientation. Journal of Rheology , 46(1):169-194, 2002. doi: 10.1122/1.1423312.
- J. Dall. Ice sheet anisotropy measured with polarimetric ice sounding radar. In 2010 IEEE International Geoscience and Remote Sensing Symposium , pages 2507-2510, 2010. doi: 10.1109/IGARSS.2010.5653528.
- A. Diez, O. Eisen, I. Weikusat, J. Eichler, C. Hofstede, P. Bohleber, T. Bohlen, and U. Polom. Influence of ice crystal anisotropy on seismic velocity analysis. Annals of Glaciology , 55(67):97-106, 2014. doi: 10.3189/2014AoG67A002.
- R. Drews, O. Eisen, D. Steinhage, I. Weikusat, S. Kipfstuhl, and F. Wilhelms. Potential mechanisms for anisotropy in ice-penetrating radar

data. Journal of Glaciology , 58(209):613-624, 2012. doi: 10.3189/ 2012JoG11J114.

- R. Drews, K. Matsuoka, C. Mart´ ın, D. Callens, N. Bergeot, and F. Pattyn. Evolution of Derwael Ice Rise in Dronning Maud Land, Antarctica, over the last millennia. Journal of Geophysical Research: Earth Surface , 120 (3):564-579, 2015. doi: 10.1002/2014JF003246. URL http://doi.wiley. com/10.1002/2014JF003246 .
- G. Durand, F. Gillet-Chaulet, A. Svensson, O. Gagliardini, S. Kipfstuhl, J. Meyssonnier, F. Parrenin, P. Duval, and D. Dahl-Jensen. Change in ice rheology during climate variations-implications for ice flow modelling and dating of the epica dome c core. Climate of the Past , 3(1):155-167, 2007. doi: 10.5194/cp-3-155-2007.
- G. Durand, O. Gagliardini, B. de Fleurian, T. Zwinger, and E. Le Meur. Marine ice sheet dynamics: Hysteresis and neutral equilibrium. Journal of Geophysical Research , 114(F3):F03009, Aug. 2009. doi: 10.1029/ 2008JF001170.
- P. Duval, M. Ashby, and I. Anderman. Rate-controlling processes in the creep of polycrystalline ice. The Journal of Physical Chemistry , 87(21): 4066-4074, 1983. doi: 10.1021/j100244a014.
- M. R. Ershadi, R. Drews, C. Mart´ ın, O. Eisen, C. Ritz, H. Corr, J. Christmann, O. Zeising, A. Humbert, and R. Mulvaney. Polarimetric radar reveals the spatial distribution of ice fabric at domes and divides in east antarctica. The Cryosphere , 16(5):1719-1739, 2022. doi: 10.5194/ tc-16-1719-2022.
- L. Favier and F. Pattyn. Antarctic ice rise formation, evolution, and stability. Geophysical Research Letters , 42(11):4456-4463, June 2015. ISSN 00948276, 1944-8007. doi: 10.1002/2015GL064195.
- S. Fujita, H. Maeno, and K. Matsuoka. Radio-wave depolarization and scattering within ice sheets: a matrix-based model to link radar and ice-core measurements and its application. Journal of Glaciology , 52(178):407-424, 2006. doi: 10.3189/172756506781828548.
- O. Gagliardini, F. Gillet-Chaulet, and M. Montagnat. A review of anisotropic polar ice models: from crystal to ice-sheet flow models. Physics of Ice Core Records II , 68, 2009.
- O. Gagliardini, T. Zwinger, F. Gillet-Chaulet, G. Durand, L. Favier, B. de Fleurian, R. Greve, M. Malinen, C. Mart´ ın, P. R˚back, J. Ruokoa lainen, M. Sacchettini, M. Sch¨ afer, H. Seddik, and J. Thies. Capabilities and performance of Elmer/Ice, a new-generation ice sheet model. Geoscientific Model Development , 6(4):1299-1318, 2013. doi: 10.5194/ gmd-6-1299-2013.

- T. A. Gerber, D. A. Lilien, N. M. Rathmann, S. Franke, T. J. Young, F. Valero-Delgado, M. R. Ershadi, R. Drews, O. Zeising, A. Humbert, et al. Crystal orientation fabric anisotropy causes directional hardening of the northeast greenland ice stream. Nature Communications , 14(1):2653, 2023. doi: 10.1038/s41467-023-38139-8.
- F. Gillet-Chaulet, O. Gagliardini, J. Meyssonnier, M. Montagnat, and O. Castelnau. A user-friendly anisotropic flow law for ice-sheet modeling. Journal of glaciology , 51(172):3-14, 2005. doi: 10.3189/ 172756505781829584.
- F. Gillet-Chaulet, O. Gagliardini, J. Meyssonnier, T. Zwinger, and J. Ruokolainen. Flow-induced anisotropy in polar ice and related ice-sheet flow modelling. Journal of non-newtonian fluid mechanics , 134(1-3):33-43, 2006. doi: 10.1016/j.jnnfm.2005.11.005.
- G. G¨ odert. A mesoscopic approach for modelling texture evolution of polar ice including recrystallization phenomena. Annals of Glaciology , 37:23-28, 2003. doi: 10.3189/172756403781815375.
- A. C. J. Henry. Simulation and post-processing code presented in 'modelling the three-dimensional, diagnostic anisotropy field of an ice rise' (Version v1). Zenodo. https://doi.org/10.5281/zenodo.13143326 , 2024a. Accessed: 2024-07-31.
- A. C. J. Henry. Code for the publication 'Predicting the three-dimensional stratigraphy of an ice rise' (Version v1). Zenodo. https://doi.org/10. 5281/zenodo.12180129 , 2024b. Accessed: 2024-06-20.
- A. C. J. Henry. Model input data, simulation output data and processed data presented in 'modelling the three-dimensional, diagnostic anisotropy field of an ice rise' (Version v1). Zenodo. https://doi.org/10.5281/ zenodo.13143591 , 2024c. Accessed: 2024-07-31.
- A. C. J. Henry. Data for the publication 'predicting the three-dimensional stratigraphy of an ice rise' (Version v1). Zenodo. https://doi.org/10. 5281/zenodo.12183293 , 2024d. Accessed: 2024-06-20.
- A. C. J. Henry, R. Drews, C. Schannwell, and V. Viˇ snjevi´. c Hysteretic evolution of ice rises and ice rumples in response to variations in sea level. The Cryosphere , 16(9):3889-3905, 2022. doi: 10.5194/tc-16-3889-2022.
- A. C. J. Henry, C. Schannwell, V. Viˇnjevi´, J. Millstein, P. Bons, O. Eisen, s c and D. R. Predicting the three-dimensional age-depth field of an ice rise. 2023. doi: 10.22541/essoar.169230234.44865946/v1.
- K. C. Jezek, J. C. Curlander, F. Carsey, C. Wales, and R. G. Barry. Ramp amm-1 sar image mosaic of antarctica, version 2, 2013.
- T. M. Jordan, C. Mart´ ın, A. M. Brisbourne, D. M. Schroeder, and A. M. Smith. Radar characterization of ice crystal orientation fabric and

- anisotropic viscosity within an antarctic ice stream. Journal of Geophysical Research: Earth Surface , 127(6):e2022JF006673, 2022. doi: https://doi.org/10.1029/2022JF006673.
- I. Koch, R. Drews, S. Franke, D. Jansen, F. M. Oraschewski, L. S. Muhle, V. Viˇ snjevi´, c K. Matsuoka, F. Pattyn, and O. Eisen. Radar internal reflection horizons from multisystem data reflect ice dynamic and surface accumulation history along the princess ragnhild coast, dronning maud land, east antarctica. Journal of Glaciology , page 1-19, 2023. doi: 10.1017/jog.2023.93.
- S. Kufner, J. Wookey, A. M. Brisbourne, C. Mart´ ın, T. S. Hudson, J. M. Kendall, and A. M. Smith. Strongly depth-dependent ice fabric in a fastflowing antarctic ice stream revealed with icequake observations. Journal of Geophysical Research: Earth Surface , 128(3):e2022JF006853, 2023. doi: 10.1029/2022JF006853.
- D. A. Lilien, N. M. Rathmann, C. S. Hvidberg, and D. Dahl-Jensen. Modeling ice-crystal fabric as a proxy for ice-stream stability. Journal of Geophysical Research: Earth Surface , 126(9), 2021. doi: 10.1029/ 2021JF006306.
- D. A. Lilien, N. M. Rathmann, C. S. Hvidberg, A. Grinsted, M. R. Ershadi, R. Drews, and D. Dahl-Jensen. Simulating higher-order fabric structure in a coupled, anisotropic ice-flow model: application to dome c. Journal of Glaciology , page 1-20, 2023. doi: 10.1017/jog.2023.78.
- M.-G. Llorens, A. Griera, P. D. Bons, I. Weikusat, D. J. Prior, E. GomezRivas, T. De Riese, I. Jimenez-Munt, D. Garc´ ıa-Castellanos, and R. A. Lebensohn. Can changes in deformation regimes be inferred from crystallographic preferred orientations in polar ice? The Cryosphere , 16(5): 2009-2024, 2022. doi: 10.5194/tc-16-2009-2022.
- Y. Ma, O. Gagliardini, C. Ritz, F. Gillet-Chaulet, G. Durand, and M. Montagnat. Enhancement factors for grounded ice and ice shelves inferred from an anisotropic ice-flow model. Journal of Glaciology , 56(199):805-812, 2010. doi: 10.3189/002214310794457209.
- C. Mart´ ın and G. H. Gudmundsson. Effects of nonlinear rheology, temperature and anisotropy on the relationship between age and depth at ice divides. The Cryosphere , 6(5):1221-1229, 2012. doi: 10.5194/ tc-6-1221-2012.
- C. Mart´ ın, R. C. A. Hindmarsh, and F. J. Navarro. On the effects of divide migration, along-ridge flow, and basal sliding on isochrones near an ice divide. Journal of Geophysical Research , 114(F2):F02006, Apr. 2009. ISSN 0148-0227. doi: 10.1029/2008JF001025. URL http://doi.wiley.com/ 10.1029/2008JF001025 .
- K. Matsuoka, T. Furukawa, S. Fujita, H. Maeno, S. Uratsuka, R. Naruse, and O. Watanabe. Crystal orientation fabrics within the antarctic ice sheet

- revealed by a multipolarization plane and dual-frequency radar survey. Journal of Geophysical Research: Solid Earth , 108(B10):2003JB002425, 2003. doi: 10.1029/2003JB002425.
- K. Matsuoka, L. Wilen, S. Hurley, and C. Raymond. Effects of birefringence within ice sheets on obliquely propagating radio waves. IEEE Transactions on Geoscience and Remote Sensing , 47(5):1429-1443, 2009. doi: 10.1109/ TGRS.2008.2005201.
- K. Matsuoka, D. Power, S. Fujita, and C. F. Raymond. Rapid development of anisotropic ice-crystal-alignment fabrics inferred from englacial radar polarimetry, central west antarctica. Journal of Geophysical Research: Earth Surface , 117(F3), 2012. doi: https://doi.org/10.1029/2012JF002440.
- K. Matsuoka, R. C. Hindmarsh, G. Moholdt, M. J. Bentley, H. D. Pritchard,
- J. Brown, H. Conway, R. Drews, G. Durand, D. Goldberg, T. Hattermann,
- J. Kingslake, J. T. Lenaerts, C. Mart´ ın, R. Mulvaney, K. W. Nicholls,
- F. Pattyn, N. Ross, T. Scambos, and P. L. Whitehouse. Antarctic ice rises and rumples: Their properties and significance for ice-sheet dynamics and evolution. Earth-Science Reviews , 150:724-745, 2015. doi: 10.1016/j. earscirev.2015.09.004.
- F. S. McCormack, R. C. Warner, H. Seroussi, C. F. Dow, J. L. Roberts, and A. Treverrow. Modeling the deformation regime of thwaites glacier, west antarctica, using a simple flow relation for ice anisotropy (estar). Journal of Geophysical Research: Earth Surface , 127(3):e2021JF006332, 2022. doi: 10.1029/2021JF006332.
- M. Montagnat, N. Azuma, D. Dahl-Jensen, J. Eichler, S. Fujita, F. GilletChaulet, S. Kipfstuhl, D. Samyn, A. Svensson, and I. Weikusat. Fabric along the neem ice core, greenland, and its comparison with grip and ngrip ice cores. The Cryosphere , 8(4):1129-1138, 2014. doi: 10.5194/ tc-8-1129-2014.
- E. C. Pettit, T. Thorsteinsson, H. P. Jacobson, and E. D. Waddington. The role of crystal fabric in flow near an ice divide. Journal of Glaciology , 53 (181):277-288, 2007. doi: 10.3189/172756507782202766.
- E. C. Pettit, E. D. Waddington, W. D. Harrison, T. Thorsteinsson, D. Elsberg, J. Morack, and M. A. Zumberge. The crossover stress, anisotropy and the ice flow law at siple dome, west antarctica. Journal of Glaciology , 57(201):39-52, 2011. doi: 10.3189/002214311795306619.
- N. M. Rathmann, C. S. Hvidberg, A. Grinsted, D. A. Lilien, and D. DahlJensen. Effect of an orientation-dependent non-linear grain fluidity on bulk directional enhancement factors. Journal of Glaciology , 67(263):569-575, 2021. doi: 10.1017/jog.2020.117.
- N. M. Rathmann, D. A. Lilien, A. Grinsted, T. A. Gerber, T. J. Young, and D. Dahl-Jensen. On the limitations of using polarimetric radar sounding to infer the crystal orientation fabric of ice masses. Geophysical Research Letters , 49(1):e2021GL096244, 2022. doi: 10.1029/2021GL096244.

- D. H. Richards, S. S. Pegler, and S. Piazolo. Ice fabrics in two-dimensional flows: beyond pure and simple shear. The Cryosphere , 16(10):4571-4592, 2022. doi: 10.5194/tc-16-4571-2022.
- D. H. Richards, S. S. Pegler, S. Piazolo, N. Stoll, and I. Weikusat. Bridging the gap between experimental and natural fabrics: Modeling ice stream fabric evolution and its comparison with ice-core data. Journal of Geophysical Research: Solid Earth , 128(11):e2023JB027245, 2023. doi: https://doi.org/10.1029/2023JB027245.
- J. M. Rignot, E. and B. Scheuchl. Measures insar-based antarctica ice velocity map, version 2, 2017. URL https://nsidc.org/data/NSIDC-0484/ versions/2 .
- C. Schannwell, R. Drews, T. A. Ehlers, O. Eisen, C. Mayer, and F. GilletChaulet. Kinematic response of ice-rise divides to changes in ocean and atmosphere forcing. The Cryosphere , 13(10):2673-2691, 2019. doi: 10. 5194/tc-13-2673-2019.
- C. Schannwell, R. Drews, T. A. Ehlers, O. Eisen, C. Mayer, M. Malinen, E. C. Smith, and H. Eisermann. Quantifying the effect of ocean bed properties on ice sheet geometry over 40 000 years with a full-Stokes model. The Cryosphere , 14(11):3917-3934, 2020. doi: 10.5194/tc-14-3917-2020.
- H. Seddik, R. Greve, T. Zwinger, and L. Placidi. A full stokes ice flow model for the vicinity of dome fuji, antarctica, with induced anisotropy and fabric evolution. The Cryosphere , 5(2):495-508, 2011. doi: 10.5194/ tc-5-495-2011.
- E. C. Smith, A. F. Baird, J. M. Kendall, C. Mart´ ın, R. S. White, A. M. Brisbourne, and A. M. Smith. Ice fabric in an antarctic ice stream interpreted from seismic anisotropy. Geophysical Research Letters , 44(8):3710-3718, 2017. doi: 10.1002/2016GL072093.
- T. Thorsteinsson, E. D. Waddington, and R. C. Fletcher. Spatial and temporal scales of anisotropic effects in ice-sheet flow. Annals of Glaciology , 37:40-48, 2003. doi: 10.3189/172756403781815429.
- I. Weikusat, D. Jansen, T. Binder, J. Eichler, S. H. Faria, F. Wilhelms, S. Kipfstuhl, S. Sheldon, H. Miller, D. Dahl-Jensen, and T. Kleiner. Physical analysis of an antarctic ice core-towards an integration of microand macrodynamics of polar ice. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences , 375(2086): 20150347, 2017. doi: 10.1098/rsta.2015.0347.
- N. H. Woodcock. Specification of fabric shapes using an eigenvalue method. Geological Society of America Bulletin , 88(9):1231, 1977. doi: 10.1130/ 0016-7606(1977)88 1231:SOFSUA 2.0.CO;2. ⟨ ⟩
- T. J. Young, C. Mart´ ın, P. Christoffersen, D. M. Schroeder, S. M. Tulaczyk, and E. J. Dawson. Rapid and accurate polarimetric radar measurements

of ice crystal fabric orientation at the western antarctic ice sheet (wais) divide ice core site. The Cryosphere , 15(8):4117-4133, 2021a.

- T. J. Young, D. M. Schroeder, T. M. Jordan, P. Christoffersen, S. M. Tulaczyk, R. Culberg, and N. L. Bienert. Inferring ice fabric from birefringence loss in airborne radargrams: Application to the eastern shear margin of thwaites glacier, west antarctica. Journal of Geophysical Research: Earth Surface , 126(5):e2020JF006023, 2021b. doi: 10.1029/2020JF006023.
A. Clara J. Henry Carlos Martín Reinhard Drews
Modelling the three-dimensional, diagnostic anisotropy field of an ice rise
Polar ice develops anisotropic crystal orientation fabrics under deformation,
yet ice is most often modelled as an isotropic fluid. We present
three-dimensional simulations of the crystal orientation fabric of Derwael Ice
Rise including the surrounding ice shelf using a crystal orientation tensor
evolution equation corresponding to a fixed velocity field. We use a
semi-Lagrangian numerical method that constrains the degree of crystal
orientation evolution to solve the equations in complex flow areas. We perform
four simulations based on previous studies, altering the rate of evolution of
the crystal anisotropy and its dependence on a combination of the strain rate
and deviatoric stress tensors. We provide a framework for comparison with radar
observations of the anisotropy field, outlining areas where the assumption of
one vertical eigenvector may not hold and provide resulting errors in measured
eigenvalues. We recognise the areas of high horizontal divergence at the ends
of the flow divide as important areas to make comparisons with observations.
Here, poorly constrained model parameters result in the largest difference in
fabric type. These results are important in the planning of future campaigns
for gathering data to constrain model parameters and as a link between
observations and computationally-efficient, simplified models of anisotropy.