Simulating higher-order fabric structure in a coupled, anisotropic ice-flow model: application to Dome C

Abstract

Ice-crystal fabric can induce mechanical anisotropy that significantly affects flow, but ice-flow models generally do not include fabric development or its effect upon flow. Here, we incorporate a new spectral expansion of fabric, and more complete description of its evolution, into the ice-flow model Elmer/Ice. This approach allows us to model the effect of both lattice rotation and migration recrystallization on large-scale ice flow. The fabric evolution is coupled to flow using an unapproximated non-linear orthotropic rheology that better describes deformation when the stress and fabric states are misaligned. These improvements are most relevant for simulating dynamically interesting areas, where recrystallization can be important, tuning data are scarce and rapid flow can lead to misalignment between stress and fabric. We validate the model by comparing simulated fabric to ice-core and phase-sensitive radar measurements on a transect across Dome C, East Antarctica. With appropriately tuned rates for recrystallization, the model is able to reproduce observations of fabric. However, these tuned rates differ from those previously derived from laboratory experiments, suggesting a need to better understand how recrystallization acts differently in the laboratory compared to natural settings.

Introduction and background

The crystal orientation fabric of ice (or simply fabric) refers to the ensemble of grain c-axis orientations (basal plane normal directions) making up the polycrystal. Ice displays a strong mechanical anisotropy; for an individual grain, shear is ~10⁴ times easier perpendicular to the c-axis than parallel to it (Duval and others, 1983). At the polycrystal scale, grain interactions reduce this effect, but single-maximum fabrics (i.e. polycrystals with a strong alignment of c-axes in one direction) have been found to shear more easily by an order of magnitude compared to isotropic fabrics (Pimienta and Duval, 1987). This induced anisotropy can have significant implications for ice flow by localizing shear near the bed (Thorsteinsson, 2001; Rathmann and Lilien, 2022a) or in the margins of ice streams (Minchew and others, 2018; Grinsted and others, 2022). Despite the strong effect that fabric can have upon flow, large-scale models of ice flow generally neglect the fabric-induced mechanical anisotropy. In part, the effects are neglected due to unknown fabrics in the dynamic regions where fabric may most affect flow, since challenging conditions have limited direct measurements of fabric in such areas (Jackson and Kamb, 1997; Thomas and others, 2021). Active and passive seismics (e.g. Bentley, 1971; Smith and others, 2017; Lutz and others, 2022) and phase-sensitive radar (pRES) (e.g. Brisbourne and others, 2019; Jordan and others, 2019) provide some constraints if surface access is possible, though depth resolution and precision are limited. In the fastest flowing regions, new methods have only recently allowed inference of fabric from airborne radar data (Young and others, 2021), thus allowing the first continuous, spatially extensive inferences of fabric. Models of fabric development can potentially be used to estimate the fabric between sparse measurements, particularly in dynamically interesting areas such as ice streams.

Modeling crystal processes

Fabric-evolution models range from small-scale (sub-meter) models that track individual ice grains to large-scale (kilometer or greater) models that rely on statistical descriptions of fabric (e.g. Gillet-Chaulet and others, 2006; Llorens and others, 2016). These models broadly include two processes: deformation that mechanically affects the distribution of grain orientations, a process termed lattice rotation, and mass movement between grains, termed recrystallization.

The effect of lattice rotation (the apparent rotation of c-axes) has been suggested to depend linearly on the deformation gradient (e.g. Svendsen and Hutter, 1996), assuming most deformation occurs as slip on the basal planes. However, this simple form does not apply at the polycrystal scale. The bulk fabric evolution depends on how the bulk stress and strain transfer to the smaller grain scale; though simple assumptions about uniform stress and strain are often used, the reality is more complex than can currently be represented by large-scale models (Castelnau and others, 1996; Rathmann and Lilien, 2022a). Nevertheless, models of fabric evolution that account for lattice rotation alone, under the assumption of homogeneous strain or stress over the polycrystal scale, are able to do a good job reproducing ice-core measurements in some locations (Gillet-Chaulet and others, 2006; Rathmann and Lilien, 2022a), suggesting that this process is well-described at least near ice divides and domes where such ice cores are drilled.

Recrystallization dominates fabric development in some regimes (e.g. Alley, 1988), but is generally neglected in large-scale fabric models. Recrystallization encompasses three processes: normal grain growth, rotation recrystallization (also called continuous dynamic recrystallization, or CDRX), and migration recrystallization (also called discontinuous dynamic recrystallization, or DDRX). Normal grain growth refers to the orientation-independent enlargement of grains near the surface of the ice sheet, and does not strongly affect fabric (Montagnat and Duval, 2000). Rotation recrystallization refers to the division of grains into subgrains during deformation, and has been found to occur in relatively shallow parts of ice sheets (e.g. Weikusat and others, 2009). However, since the new grains are split from their parents gradually, they do not lead to a significantly different preferred orientation, but rather cause a slowing of fabric development (Montagnat and Duval, 2000), which is sometimes modeled as a diffusive process on the orientation of the grains (Gödert, 2003; Placidi and others, 2010; Richards and others, 2021). Migration recrystallization refers to the movement of mass to (possibly new) grains well-oriented for the in situ stress (Chapelle and others, 1998). Though early evidence (Alley, 1988; Duval and Castelnau, 1995) suggested that migration recrystallization only affects warm ($\ge -10^{\circ}$C), deep ice, some ice-core evidence has indicated that the process occurs shallower (Diprinzio and others, 2005; Kipfstuhl and others, 2006, 2009). These observations led to a proposal of a new paradigm for the activation of migration recrystallization at sufficient stress, even if ice is cold (Faria and others, 2014). The bulk effect of migration recrystallization on fabric can be simulated as a production/decay process that leads to grains well-oriented for the in situ stress or strain rate (Gödert and Hutter, 1998; Faria and others, 2003; Placidi and others, 2010) and the necessary model parameters have recently been tuned to experimental data (Richards and others, 2021).

Current anisotropic ice-flow models

Models of fabric development are, with few exceptions, separate from ice-flow models. Only a few large-scale ice-flow models consider the effect of fabric on flow, let alone track its development explicitly. To the extent that the effect of fabric is included, it is often as a strain-rate enhancement factor, which can be isotropic or anisotropic (i.e. the material response may or may not be invariant to rotations from a reference state) depending on how the scalar factor is calculated (Seddik and others, 2008; Graham and others, 2018). Since, in practice, much ice-flow modeling infers a scalar enhancement factor prior to prognostic simulations (e.g. Larour and others, 2005; Arthern and Gudmundsson, 2010), the influence of fabric may be subsumed into the bulk viscosity along with other factors causing enhancements, e.g. variations in mechanical damage or temperature. Such scalar compensation may successfully account for the instantaneous effect of anisotropy on strain in the dominant direction (Rathmann and Lilien, 2022a), but by definition cannot accurately describe the viscosity in all directions; if conditions change or if multiple deformation modes are active, scalar compensation may lead to incorrect transient results.

Other models account for anisotropy by assuming that fabric adjusts instantaneously to the steady state fabric pattern consistent with in situ stress and strain (Graham and others, 2018), and thus obviate the need for a time-evolving fabric model. However, this assumption is poorly supported for many parts of ice-sheet interiors, where fabric is thought to develop slowly by lattice rotation (Alley, 1988). For example, the stress state is unconfined compression at all depths beneath a perfect dome, but fabric varies from isotropic at the surface to a single maximum at depth (e.g. Durand and others, 2009). To determine fabric accurately everywhere, lattice rotation and recrystallization must both be accounted for, which necessitates a model of fabric development that includes transient and advective effects (e.g. Thorsteinsson and others, 2003).

Finally, one large-scale model, Elmer/Ice, is capable of simulating fabric evolution explicitly and couples that evolution to flow (Gillet-Chaulet and others, 2006; Seddik and others, 2008). In that module, fabric is assumed to develop solely due to lattice rotation, and the closure approximation, which relieves the need to track infinitely detailed information about fabric orientation, contains implicit site-specific tuning; the effect of this closure approximation is essentially untested due to the lack of closure-approximation-free model for comparison. The fabric evolution module is not commonly used despite the widespread applications of Elmer/Ice. This limited use stems from the added computational cost of simulating the fabric, more free model parameters typically unconstrained by observations, and violation of the lattice-rotation-only assumption in fast-flowing areas where the ice-flow model is most commonly applied.

Fabric representation

To model fabric development and its effect on bulk directional ice viscosities, a mathematical description of grain orientations is needed. A full description of the fabric would contain information about the orientation, size, shape and number of grains that makes up the polycrystal, properties that can be simulated explicitly using small (sub-meter scale) microstructural models (e.g. Durand, 2004). In practice, large-scale fabric models generally ignore the size and shape of grains, and explicitly track the orientation alone (Gillet-Chaulet and others, 2006; Seddik and others, 2008; Richards and others, 2021). Large-scale models of fabric development rely on the distribution of preferred grain orientations in order to construct grain-weighted-average quantities such as directional viscosities (e.g. Rathmann and Lilien, 2022b). Different weights can be used to construct average quantities, such as the distribution of number of grains or their mass in orientation space, S ² (surface of the unit sphere).

Orientation distribution functions

In the simplest approach, the orientation distribution gives equal weight to all grains, defined as the orientation distribution function (ODF):

(1)$${\rm ODF} = {\psi( \theta,\; \phi) \over N},\; $$

where ψ is the c-axis orientation density, θ and ϕ are the polar and azimuth angles, respectively, $N = \int _{S^2} \psi \, {\rm d}\Omega$ is the total number of grains and the differential of the solid angle is dΩ = sin (θ) dθ dϕ.

An alternative orientation distribution may be defined by regarding the grain mass density as a function of orientation, ϱ*(θ,  ϕ), termed the orientation mass density, a concept that arises from the theory of mixtures of continuous diversity (Faria, 2001). By analogy to the ODF, we define a mass orientation distribution function (MODF), which is ϱ* normalized:

(2)$${\rm MODF} = {\varrho^\ast ( \theta,\; \phi) \over \rho_{\rm i}},\; $$

where $\rho _{\rm i} = \int _{S^2}\varrho ^\ast \, {\rm d}\Omega$ is simply the density of ice.

Structure tensors and closure

An arbitrary (M)ODF will not necessarily have a closed form, so large-scale models have previously proposed several mathematical representations of fabric, perhaps the simplest of which is the cone- and girdle-angle model (e.g. Thorsteinsson, 2001). More complex models have relied on the fact that any (M)ODF can be expanded in terms of a structure-tensor series (e.g. Svendsen and Hutter, 1996), useful for directly calculating bulk directional viscosities. The structure tensors are the moments of the fabric orientation distribution, defined as the sum of the weighted m-fold outer product of grain c-axes with themselves. For a discrete set of N grains with c-axis denoted c_i, it is defined as

(3)$${\bf a}^{( m) } = \sum_i^N w_i \, ( {\boldsymbol c}_i\otimes) ^m,\; $$

where (c_i ⊗ )^m is the m-repeated outer product, and m is the (even) order of the tensor (Advani and Tucker, 1987). The weights w _i are normalized such that $\sum _i^N w_i = 1$. If no shape/size information is available, then all grains are typically weighted equally (w _i = 1/N) so (3) becomes the arithmetic mean, and the a^(m) are moments of the ODF. If the w _i are taken to be the volume or mass fraction (i.e. the MODF itself), then the a^(m) are the moments of the MODF.

In practice, the structure-tensor series representation of fabric is often truncated at the second-order contribution a⁽²⁾ (i.e. excluding a^(m) for m > 2) (Gillet-Chaulet and others, 2006; Seddik and others, 2008). a⁽²⁾ alone is sufficient to represent simple yet common fabrics such as single-maximum and girdle fabrics found in deep ice cores (Gow and Williamson, 1976; Advani and Tucker, 1987; Weikusat and others, 2017; Bauer and Böhlke, 2021). Such ice-core measurements often discard (or cannot recover) the principal directions and retain only the eigenvalues of a⁽²⁾, here ordered such that λ ₃ ≥ λ ₂ ≥ λ ₁ (though conventions differ).

Complicated fabrics that cannot be expressed in terms of a⁽²⁾ alone have also been observed. For example, ice cores from dynamic areas (Jackson and Kamb, 1997; Gerbi and others, 2021) and laboratory experiments (Journaux and others, 2019; Fan and others, 2020) show multi-maxima and ‘diamond’ fabrics; a⁽²⁾ cannot capture such structure. Tracking higher-order structure tensors would allow more complex fabrics to be modeled, but this quickly becomes challenging for technical reasons, including degeneracy of the tensors (e.g. the information contained in a⁽²⁾ is included in a⁽⁴⁾; Bauer and Böhlke, 2021), complications of tracking moments of the fabric rather than the fabric itself, and the lack of closure for lattice rotation (the evolution of a⁽²⁾ depends on a⁽⁴⁾, that of a⁽⁴⁾ on a⁽⁶⁾, etc.).

Expansion series approach

Some of the issues of the structure-tensor approach can be avoided by describing the fabric as a series of increasingly finer anomalies from isotropy in spectral space. Richards and others (2021) and Rathmann and others (2021) independently proposed using a spherical harmonic expansion of ϱ* or ψ for polycrystalline viscoplastic problems. This approach has long been used for elastic problems (Turner, 1999), but is new to glaciology. In this formulation, an arbitrary ϱ* can be described by a harmonic series expansion

(4)$$\varrho^\ast ( \theta,\; \phi) = \sum_{l = 0}^{L}\sum_{m = -l}^{l} \hat\varrho_{l}^{m} Y_{l}^{m}( \theta,\; \phi) ,\; $$

where L is the order of the approximation (spectral truncation), $\hat \varrho _{l}^{m}$ are the complex expansion coefficients and $Y_{l}^{m}$ the spherical harmonic functions. Since ϱ* is antipodally symmetric, $\hat \varrho _l^m = 0$ if l is odd. This approach describes the MODF (ϱ*/ρ _i) directly, without parameterization and has recently been applied to model individual ice parcels (Rathmann and others, 2021; Richards and others, 2021) and to a simplified, slab model of ice flow (Rathmann and Lilien, 2022a).

Physical and technical advantages are discussed below. Here, it is sufficient to note that, for a given L, expansion (4) contains the same information as structure tensors up to order L (Rathmann and others, 2021).

Modeling the ODF or MODF

Here, like Richards and others (2021), we consider ϱ*, since it provides a framework where recrystallization can naturally be modeled as the transfer of mass between grains with different orientations. Indeed, when recrystallization is active, the total number of grains is not conserved (i.e. ∂N/∂t ≠ 0) unlike total mass (per unit volume). However, Rathmann and Lilien (2022a) proposed using the same DDRX model math for representing the effect on the ODF, in which case grains (orientations) are understood to spontaneously decay and nucleate in equal proportion so that N is conserved. Since, ψ and ϱ* are mathematically identical up to their normalization in these two models, the MODF and ODF are identical, and hence modeled a^(m), too, although by formulating the model in terms of ϱ* it rests on stronger physical grounds.

Motivation: the need for a higher-order model

In this work, we adopt the spherical harmonic expansion of ϱ*, allowing arbitrary grain orientation distributions to be represented. Before proceeding, we demonstrate that such an approach is necessary because the existing structure-tensor-based representations (and closure schemes) are inadequate for modeling migration recrystallization. Although such models are usually formulated in terms of the ODF, the limitations described next would persist even if the models were reinterpreted in terms of the MODF.

When lattice rotation dominates fabric evolution, the evolution of a⁽²⁾ can be written in closed form by parameterizing a⁽⁴⁾ in terms of a⁽²⁾. Figure 1 makes this clear, showing the relationship between the unique lowest-order expansion coefficients of the MODF, $\hat \varrho _2^0$ and $\hat \varrho _4^0$, derived from ice cores and lab experiments (shown by hollow markers), in a frame where fabric is approximately rotationally symmetry around the vertical axis (by first rotating the fabric before comparison, no generality is lost by considering such a reference frame). Since the corresponding unique lowest-order tensorial components are linearly related to $\hat \varrho _2^0$ and $\hat \varrho _4^0$, a tensorial correlation between a⁽⁴⁾ and a⁽²⁾ likewise exists. Gillet-Chaulet and others (2005) noted this by considering an analytical model of fabric development, for multiple types of deformation, leading to the tensorial closure function a⁽⁴⁾(a⁽²⁾) currently used by Elmer/Ice's anisotropic flow solver. Indeed, the colored lines in Figure 1 show that a zero-dimensional model (Rathmann and others, 2021) considering only lattice rotation produces a single consistent relationship between second- and fourth-order terms of the fabric under simple shear, unconfined extension or unconfined compression, confirming that there is a function a⁽⁴⁾(a⁽²⁾) that applies regardless of the type of deformation. Thus, if one is concerned only with the evolution of second-order moments of the fabric under lattice rotation, the closure approximation of Gillet-Chaulet and others (2005) is sufficient and higher-order fabric representations provide little benefit.

Figure 1. Relationship between normalized expansion coefficients $\hat \varrho _2^0/\hat \varrho _0^0$ and $\hat \varrho _4^0/\hat \varrho _0^0$ in zero-dimensional modeling (lines) and from observations (markers). All fabric states are rotated into a (nearly) vertically symmetric reference frame. The two components fully capture the second- and fourth-order fabric strength in the case of vertical symmetry, where the white/colored region represents the space of possible fabric states (MODFs) with a vertical symmetry. Lines show modeled fabric evolution using SpecFab (Rathmann and others, 2021) with lattice rotation alone (brown/green) and with DDRX alone (red). Increasing the strength of recrystallization causes greater deviation from the lattice rotation trajectories, eventually resulting in a steady relationship between $\hat \varrho _2^0$ and $\hat \varrho _4^0$ regardless of deformation amount (solid red circle). Hollow markers show fabrics observed in the laboratory and in ice cores. The marker shape indicates the deformation regime: triangles for extension, diamonds and crosses for compression and squares for simple shear. Color indicates source for each dataset: blue from Thorsteinsson and others (1997), pink from Treverrow and others (2016), purple from Westhoff and others (2021), orange from Voigt (2017), green from Thomas and others (2021) and yellow from Qi and others (2019). Crosses indicate DDRX-affected fabrics during relatively warm deformation tests: red markers from Fan and others (2020), and gray markers from Hunter and others (2023). Ball plots show the corresponding MODFs at different points in the state space. Gray contours show modeled enhancement factors for vertical compression/extension, E _zz. Any vertical spread in observations at a single $\hat \varrho _2^0/\hat \varrho _0^0$ cannot be captured by a traditional closure approximation.

Migration recrystallization makes the correlation function multi-valued, preventing use of a closure approximation. Figure 1 demonstrates this, showing that DDRX-affected fabrics measured in the lab (red and gray markers) and an existing DDRX model (red lines) cause fabric-state trajectories to diverge from the approximate correlation found for lattice-rotation-induced fabrics. Hence, $\hat \varrho _4^0$ is not always uniquely defined by $\hat \varrho _2^0$ (or a⁽⁴⁾ by a⁽²⁾). Thus, since the evolution of a⁽²⁾ depends on a⁽⁴⁾, a⁽⁴⁾ must be allowed to freely evolve, too, when migration recrystallization is active.

To some degree, the spread in fabric observed in laboratory and ice-core data suggests that the effect of recrystallization is always present (Fig. 1). These data thus suggest that migration recrystallization should always be included in fabric-evolution models, even if lattice rotation is the primary control on fabric development in many areas. Including this effect is particularly important for understanding the implications of the fabric upon ice flow, since the directional viscosities are thought to depend on both a⁽²⁾ and a⁽⁴⁾ (e.g. Gillet-Chaulet and others, 2006; Rathmann and Lilien, 2022b). The gray contours in Figure 1 show the modeled compressional/extensional enhancement factor, E _zz, along the z-axis (the viscosity model is introduced below), indicating fabric-induced hardness may be overestimated by $\sim \! 40\%$ if DDRX is active (red and gray markers) but a correlation based on lattice rotation is assumed (green model line). In such a case, using a model of fabric development to constrain directional viscosity is only useful if the model considers migration recrystallization as well as lattice rotation; otherwise, an isotropic ice-flow model has similar error.

Including fourth-order terms, and consequently modeling recrystallization, requires either including higher-order tensors in structure-tensor-based models, and thus a closure approximation of a⁽⁴⁾ in terms of a⁽⁶⁾ or a similar higher-order closure, or a switch to another approach. Here, we incorporate the above spectral description of fabric as a new module in Elmer/Ice (Zwinger and others, 2007; Gagliardini and others, 2013), a state-of-the-art ice-flow model, providing an alternative to the current structure-tensor-based fabric representation. In coupling to flow, we also implement the recently derived, unapproximated nonlinear orthotropic rheology, validated against Dye-3 deformation tests (Rathmann and Lilien, 2022b), which is expected to be important when fabric and stress are misaligned.

We consider a higher-order approximation of the fabric (L = 6, but not limited thereto), and model the (transient) effect of both lattice rotation and migration recrystallization on large-scale fabric development and ice flow, while removing any closure approximations and resulting error in the viscosity. Using rates of rotation and migration recrystallization tuned to laboratory measurements and ice-core data, we apply the new nonlinear, tensorially anisotropic, thermally coupled ice-flow model to a realistic geometry. While the added complexity of this model is most relevant for more dynamic areas, the dearth of data prevents meaningful validation in such regions. Instead, we model a transect across Dome C, East Antarctica and compare to ice-core (Durand and others, 2009) and phase-sensitive radar (Corr and others, 2021; Ershadi and others, 2022) measurements of fabric there.

Notation

Throughout, bold symbols indicate tensors and vectors. The operator $\nabla _{S^2}$ denotes the gradient in orientation space (S ²), while $\nabla$ is the gradient in Cartesian space. The subscript _g denotes a grain parameter. Other operators follow standard conventions; a full list of symbols is given in Table 1.

Table 1. Notation used in the text

Model

We are interested in modeling fabric development and its effect upon ice flow. This is a coupled problem since the fabric evolution depends on the stress and strain rates, and vice versa. Moreover, since recrystallization and the isotropic part of the viscosity are temperature dependent, the temperature field must be simulated or specified. Finally, the ice geometry (i.e. a free surface) must be updated based on mass conservation. We first present the model equations, then describe the application to Dome C.

Fabric model

We account for the combined effects of lattice rotation, migration recrystallization and rotation recrystallization in a single model of fabric development. Specifically, we follow Placidi and others (2010); Rathmann and Lilien (2022a) and Richards and others (2021) by modeling the rate-of-change of ϱ*(x,  t,  θ,  ϕ) as

(5)$${{\rm D}\varrho^\ast \over {\rm D}t} = \underbrace{-\nabla_{S^2}{\boldsymbol \cdot}( \varrho^\ast \dot{{\boldsymbol c}}) }_{{\rm Lattice\ Rotation}} + \underbrace{\Gamma\varrho^\ast }_{{\rm DDRX}} + \underbrace{\Lambda_0\nabla_{S^2}^2\varrho^\ast }_{{\rm CDRX}},\; $$

where

(6)$${{\rm D}\varrho^\ast \over {\rm D}t} = {\partial\varrho^\ast \over \partial t} + {\boldsymbol u}\,{\boldsymbol \cdot}\nabla\varrho^\ast $$

is the material derivative, u(x,  t) is the bulk velocity, $\dot {{\boldsymbol c}}( \theta ,\; \, \phi ,\; \, {\boldsymbol u})$ is the rate-of-change in orientation of an individual c-axis due to lattice rotation (c-axis angular velocity field on S ²), Γ(θ,  ϕ,  τ) is an orientation-dependent mass production/loss rate describing the effect of migration recrystallization given the deviatoric stress tensor τ and Λ₀ is the scalar rate of rotation recrystallization. We address each of these terms in turn.

Lattice rotation

Let ${\dot {\boldsymbol \epsilon }} = ( \nabla {\boldsymbol u} + \nabla {\boldsymbol u}^{\intercal} ) / 2$ be the bulk strain-rate tensor, and ${\boldsymbol W} = ( \nabla {\boldsymbol u} - \nabla {\boldsymbol u}^{\intercal} ) / 2$ the bulk spin tensor. If lattice rotation is regarded as a kinematic process (i.e. grains reorient passively), the apparent rotation of an individual c-axis can be modeled as

(7)$$\dot{{\boldsymbol c}} = {\boldsymbol W} {\boldsymbol \cdot}\, {\boldsymbol c} - \iota\left[{{ \dot{\boldsymbol \epsilon}}}_{\rm g} {\boldsymbol \cdot}\, {\boldsymbol c}- ( {\boldsymbol c}^\intercal{\boldsymbol \cdot}{{\dot{\boldsymbol \epsilon}}}_{\rm g}{\boldsymbol \cdot}\, {\boldsymbol c}) {\boldsymbol c}\right],\; $$

where c = [sin (θ)cos (ϕ),  sin (θ)sin (ϕ),  cos (θ)], ${\dot {\boldsymbol \epsilon }}_{\rm g}$ is the strain rate at the grain scale, and $\iota$ is a parameter determining the ratio of basal slip to rigid body rotation (e.g. Svendsen and Hutter, 1996). As discussed in the section on calibration below, we take $\iota = 1$, equivalent to assuming that the deformation in simple shear behaves like a deck of cards (following, e.g. Gagliardini and Meyssonnier, 1999). In practice, only ${\dot {\boldsymbol \epsilon }}$ is known/computed, while individual grains are subject to a heterogeneous strain-rate field, ${{\dot {\boldsymbol \epsilon }}}_{\rm g}$, where ${{\dot {\boldsymbol \epsilon }}}_{\rm g} = {\dot {\boldsymbol \epsilon }}$ is not guaranteed. This heterogeneity is not easily modeled, so two assumptions are common: the Sachs hypothesis, that stress is homogeneous, and the Taylor hypothesis, that the strain rate is homogeneous. Here, we use an intermediate approximation following Gillet-Chaulet and others (2006), namely

(8)$${\dot{\boldsymbol \epsilon}}_{\rm g} = ( 1-\alpha_{{\rm LR}}) {\dot{\boldsymbol \epsilon}} + \alpha_{{\rm LR}} {1\over \eta_0} {\boldsymbol \tau},\; $$

where α _LR ∈ [0; 1] is an interaction parameter dictating the compromise between the Taylor and Sachs hypotheses and η ₀ is the nonlinear viscosity (defined in Eqn (18)). We take α _LR = 0.06 following previous work (Ma and others, 2010; Lilien and others, 2021), though we expect the results to be insensitive to this choice for α _LR < 0.1 (Martín and others, 2009).

Migration recrystallization

Migration recrystallization is modeled as an orientation-dependent production/decay process proportional to the square of the shear stress resolved on basal planes, similar to Placidi and others (2010). In effect, this parameterization compares the fabric state to the state most favorable for deformation, causing the MODF to decrease in directions that are unfavorable (little basal-plane resolved shear stress) and increase in favorable directions (high basal-plane resolved shear stress). The production/decay rate, Γ, is taken to be

(9)$$\Gamma( \theta,\; \phi,\; {\boldsymbol \tau}) = \Gamma_0\left(D( \theta,\; \phi,\; {\boldsymbol \tau}) -\langle D( \theta,\; \phi,\; {\boldsymbol \tau}) \rangle\right),\; $$

where Γ₀ is a rate factor, and the deformability, D, is given by (Placidi and others, 2010; Richards and others, 2021)

(10)$$D( \theta,\; \phi,\; {\boldsymbol \tau}) = {( {\boldsymbol \tau}{\bi\, \colon \, }{\boldsymbol \tau}) {\bi\, \colon \, }{\boldsymbol c}^{( 2) } - {\boldsymbol \tau}{\bi\, \colon \, }{\boldsymbol c}^{( 4) }{\bi\, \colon \, }{\boldsymbol \tau}\over {\boldsymbol \tau}{\bi\, \colon \, }{\boldsymbol \tau}}.$$

D describes the favorable/unfavorable c-axis directions and is dictated solely by the stress, while 〈D〉 describes how favorably the current fabric configuration is, which depends on the local second- and fourth-order structure tensors and τ.

This formulation differs slightly from previous work (Placidi and others, 2010; Richards and others, 2021) by using the stress tensor, rather than the strain-rate tensor, to determine the deformability (and thus the directions in which grain growth occurs). This is simply an alternative parameterization of the underlying process, which is thought to depend on a number of microphysical parameters including dislocation density in addition to τ and ${\dot {\boldsymbol \epsilon }}$ (Faria, 2006a, 2006b). Since previous work considered a scalar viscosity, τ and ${\dot {\boldsymbol \epsilon }}$ were co-axial, and hence using either in Eqn (10) was equivalent. Because migration recrystallization is thought to depend on stress rather than the strain rate (e.g. Duval and Castelnau, 1995; Chapelle and others, 1998), the distinction is relevant when coaxiality cannot be assumed, such as is the case for the bulk flow law we introduce below. However, we note that while the using τ versus ${\dot {\boldsymbol \epsilon }}$ in Eqn (10) will have an effect on the fabric development, neither choice is superior to the other, and a more complete treatment would likely require a full microstructural model (e.g. Faria, 2006b).

Finally, we are left to determine an overall rate for this process. While the rate as well as the direction in orientation space of migration recrystallization is thought to depend on several microphysical properties (Faria, 2006a, 2006b), we retain the assumption of previous work (Richards and others, 2021) that the total amount of recrystallization is related only to the total strain, not rate or stress at which the strain occurred. We thus parameterize the overall rate factor as

(11)$$\Gamma_0 = {{\dot{\epsilon}}}_{\rm E}A_\Gamma {\rm e}^{-Q_\Gamma/( RT) },\; $$

where $A_\Gamma$ and $Q_\Gamma$ are tunable parameters, ${{\dot {\epsilon }}}_{\rm E}$ is the effective strain rate and R is the universal gas constant. This differs slightly from Richards and others (2021), which used a linear dependence on temperature, based upon laboratory experiments, while we use an Arrhenius relationship based on ice-core data indicating rapid activation above $-10^\circ$C (Alley, 1988; Duval and Castelnau, 1995). The Arrhenius activation is used to make the model applicable closer to the melting point; within the regime of temperatures considered by Richards and others (2021), the Arrhenius relationship is close to linear and the approaches are essentially equivalent.

Rotation recrystallization

Rotation recrystallization is modeled as a diffusive process in orientation space (Placidi and others, 2010) with diffusion coefficient (or rate factor) Λ₀. This is a phenomenological description, consistent with the definition of rotation recrystallization in that new grains are created by splitting off from parent grains with similar orientations. Since ice-core evidence indicates a gradual activation of rotation recrystallization compared to migration recrystallization, we follow Richards and others (2021) exactly by setting

(12)$$\Lambda_0 = {{\dot{\epsilon}}}_{\rm E}\left(m_\Lambda T + b_\Lambda\right),\; $$

where $m_\Lambda$ and $b_\Lambda$ are the slope and intercept of the temperature relationship.

Fabric numerics

Numerically, we represent ϱ* using the spectral expansion in Eqn (4) with L = 6, leading to 28 $\hat \varrho _l^m\in {\opf C}$. The problem is solved on a finite-element mesh, in two or three dimensions. The model mesh is Eulerian, so the fabric evolution has an advective component (last term in Eqn (6)) and a reactive component (right-hand side of Eqn (5)). In our implementation, these two contributions are handled separately. Dϱ*/Dt at each node is determined from the current fabric state, ϱ*, using the SpecFab code of Rathmann and Lilien (2022a). The advective term is handled using the existing code in Elmer/Ice with using residual-free bubble elements for stabilization.

We treat the real and imaginary components of $\hat \varrho _l^m$ separately, leading to 56 values to track per computational node. However, not all coefficients are independent; we leverage the Hermitian (anti)symmetry of the spectral expansion

(13)$$\hat\varrho_{l}^{-m} = ( -1) ^m\left(\hat\varrho_{l}^{m}\right)^\ast ,\; $$

which allows direct calculation of the coefficients m < 0 from those with m > 0. Moreover, the fabric model (Eqn (5)) conserves mass, and we assume incompressibility, so ρ _i is constant and thus $\hat \varrho ^0_0$ is constant as well since $\rho _{\rm i} = \sqrt {4\pi }\hat \varrho ^0_0$. Taking the above considerations into account, 27 unique coefficients are left (though this discussion focuses on L = 6, our Elmer/Ice implementation automatically handles arbitrary L while only computing the unique, independent components). Finally, when ∂u/∂y = 0 (generally equivalent to when the domain is two dimensional (2-D)), a further 12 components are identically zero (although this does not mean that the fabric is 2-D, only that it has useful symmetries, which cause some coefficients to be null). Thus, for L = 6, there are 15 fabric components to consider in 2-D and 27 in three-dimensional (3-D).

Fabric evolution is determined using an implicit timestepping scheme, iterating coefficient-wise. The key difference compared to small-scale spectral models (Rathmann and others, 2021; Richards and others, 2021) is that advection greatly increases the problem dimension. That issue is manifest most strongly in 3-D, where a large number of coefficients, combined with many neighbors for a node in the mesh, means that the matrix problem to be solved has a prohibitively wide bandwidth if all fabric coefficients are solved for simultaneously. In order to set up a scheme that would allow an arbitrary L, we thus found it necessary to use a coefficient-wise solution similar to the existing structure-tensor-based fabric module in Elmer/Ice (Gillet-Chaulet and others, 2006).

We find that artificial diffusion is necessary in both orientation space (when lattice rotation is active) and in Cartesian space (when migration recrystallization is active) to regularize the problem. For stabilization in orientation space, we use hyper diffusion, which reduces noise in higher-order harmonics while having less effect on lower-order coefficients; the diffusion coefficients were tuned for L = 6 so that, in the steady state, $\hat \varrho _2^m$ are unaffected, whereas $\hat \varrho _4^m$ are up to 40% smaller compared to L ≥ 20 when subject to unconfined compression. How this affects the ability to simulate strong single maxima, thus limiting the magnitude of simulated bulk directional enhancement factors, is treated in Appendix A. For the stabilization in Cartesian space, we introduce a diffusion term of the form $\zeta _0\nabla ^2\hat \varrho _l^m$, where ζ ₀ = 5.0 × 10⁻³ a⁻¹ m⁻¹, to the left-hand side of Eqn (5). In theory, this diffusion limits the maximum spatial gradient in fabric, but since fabric transitions slowly this is not a problem in practice; the term only serves to prevent numerical noise in near-zero coefficients from growing during long simulations. Demonstrations of fabric evolution in simple cases, and comparison to existing models, are shown in Appendix A.

Calibration of model parameters

Before applying the model, it is necessary to constrain the rate factors for migration and rotation recrystallization. In order to evaluate the suitability of different rate factors, we compare the fabric predicted by a simplified zero-dimensional model with that observed in the EPICA Dome C (EDC) core (Durand and others, 2009). In this calibration, the core was assumed to be drilled at a perfect dome (undergoing unconfined compression), with constant pure shear from surface to bed (i.e. a Nye model of dome strain; Nye, 1963). We assumed steady state, so that the fabric evolution of a single, zero-dimensional ice parcel at increasing depth gives the modeled fabric profile of the core. At the surface, the downward velocity was set to be 1.53 cm a⁻¹, the mean accumulation throughout the ice core (Bazin and others, 2013), with zero velocity at the bed. Borehole temperature measurements were used to determine T in order to calculate Λ₀ and Γ₀ (Buizert and others, 2021). The fabric was assumed to be a weak single maximum with $a^{( 2) }_{zz} = 0.44$ at 214 m depth, matching the first measurement of EDC (Durand and others, 2009). Below, fabric evolution followed Eqn (5), under the assumption that the stress and strain are coaxial so that migration recrystallization can be modeled without knowing the rheology. For this purpose, we must assume that α _LR = 0 (i.e. the Taylor hypothesis), since τ cannot be calculated in this zero-dimensional model. The model was implemented purely using SpecFab (Rathmann and others, 2021), and solved for all complex fabric coefficients simultaneously using forward-Euler timestepping with 200–3200 year time steps.

For evaluation, the modeled fabric was compared to fabric measured on vertical thin sections cut every 11–50 m along the EDC core (Durand and others, 2009). Although the orientation of the core and therefore the thin sections are unknown, we can still compare the three fabric eigenvalues and identify the vertical component. Moreover, for this simple model, the two horizontal eigenvalues are equal, so there is no need to identify the orientation of the core. For this comparison, we assume that the ODF and MODF are equivalent, i.e. that the grain size distribution is independent of c-axis direction.

We consider two sets of model parameters: one tuned to laboratory data and one to ice-core data (Fig. 2). For the ‘laboratory’ values, we calibrated to the same data as Richards and others (2021), and used their calibrated Λ₀ = 0.001 × T + 0.21, but took $\iota = 1$, and assumed that Γ₀ is defined by Eqn 11. To find $A_\Gamma$ and $Q_\Gamma$, the parameters for this relationship, we used the observations from Richards and others (2021), Table 1, but fit an exponential relationship rather than a linear one. The resulting values are $A_\Gamma = 1.91\times 10^7$ and $Q_\Gamma = 3.36\times 10^{4}$ J mol⁻¹. This set of parameters leads to an RMSE of 0.19 (27% error) in the largest eigenvalue compared to observations (Fig. 2, full blue line). Next, again assuming $\iota = 1$, we found a set of ‘ice-core’-calibrated parameters by brute-force optimizing $b_\Lambda$ and $A_\Gamma$ to minimize the misfit to the EDC data. To do so, we assumed that the temperature dependence of rotation and migration recrystallization are accurately described by the laboratory-calibrated dependencies (i.e. $m_\Lambda = 1.26\times 10^{-3}$ ($^\circ$C)⁻¹ and $Q_\Gamma = 3.36\times 10^{4\circ }$ as above). A 2-D grid of values for $b_\Lambda$ and $A_\Gamma$ was then explored, first varying by orders of magnitude around the laboratory values and subsequently refining once the order of magnitude was determined. We note that, in this scheme, if $b_\Lambda$ was small, Λ₀ could be negative; since this is physically implausible, Λ₀ was assumed to be zero under such conditions. The optimized value for $A_\Gamma$ was $A_\Gamma = 4.3\times 10^7$, which is of the same order as the laboratory-derived value. The optimal value of $b_\Lambda$, however, was 0.02, leading to Λ₀ being 0 in the top 1000 m and negligible below that. Physically, rotation recrystallization is thought to be most important in the upper portions of the ice sheet, so the parameterization makes little sense with this recalibrated value. For the ‘ice-core’-calibrated parameters, we thus took Λ₀ = 0. These parameters led to an RMSE of 0.08 in the largest eigenvalue (11% misfit; dashed blue line in Fig. 2), substantially better than that from the laboratory-calibrated values. Moreover, this misfit was highest near the bottom of the ice core, where substantial spread in the measurements prevents better fitting (Fig. 2).

Figure 2. Rate factor calibration. (a) Modeled eigenvalues, using a zero-dimensional model, resulting from different rate factors. Colors indicate eigenvalue number (blue for λ ₃, orange for λ ₂ and green for λ ₁). Squares show measured fabric (Durand and others, 2009). Lines show modeled fabric, using the zero-dimensional model, with lab-calibrated, ice-core-calibrated and Richards and others (2021) (R2021) parameters indicated by the solid, dashed and dotted lines, respectively. (b) Temperature in the EDC borehole, from Buizert and others (2021).

For completeness, we also modeled the ice-core fabric using the calibration provided by Richards and others (2021) without modification: $\iota = 0.026\times T + 1.95$, Λ₀ = 0.001 × T + 0.21, and Γ₀ = 0.176 × T + 6.09. This leads to an RMSE of 0.18 in the largest eigenvalue (27% error; dotted blue line in Fig. 2), similar to using the ‘laboratory’ calibration. The difference between this and our ‘laboratory’ calibration stems almost exclusively from the difference in $\iota$; the exponential dependence of Γ₀ on temperature has a relatively minor effect. This third set of parameters was not used for large-scale modeling, since it did not produce a good fit to the data and our implementation of the large-scale model assumes $\iota = 1$.

In addition to the calibration discussed above, we tested sensitivity to three alternative calibration schemes: one restricting the depths at which the model was tuned to <3000 m and two using a Dansgaard–Johnsen model of strain (Dansgaard and Johnsen, 1969) rather than a Nye model, also tuned to data <3000 m depth. The restriction in depth avoids possible issues due to a non-dome-like stress state that arguably exists near the bed beneath Dome C (Durand and others, 2007). The Dansgaard–Johnsen model uses a somewhat more realistic description of deformation beneath the divide (specifically that the strain rate goes linearly to zero below a certain depth, dubbed the ‘kink height’). Since the height of the kink is poorly constrained, we tested both 0.2 and 0.4 times the ice thickness. Restricting calibration to shallow depths resulted very similar misfit in those depths to the full-depth inference, and rate factors differed by <$5\%$. With either kink height, the Dansgaard–Johnsen model resulted in a smaller model-data misfit ($\sim \!7\%$) than the Nye model, but the inferred parameters did not differ significantly from those obtained using the Nye model (difference of $< \!5\%$). The results of the additional calibrations are shown in Supplementary Figure S1. While there may be benefits to using one of these alternative calibrations, the full-depth Nye calibration arguably relies on the fewest assumptions, and thus it is what we use in subsequent simulations.

Ice flow

Ice flow is an incompressible Stokes flow

(14)$$\nabla \,{\boldsymbol \cdot}\, {\boldsymbol u} = 0,\; $$

governed by the momentum balance

(15)$$\nabla \,{\boldsymbol \cdot}\, {\boldsymbol \sigma} + \rho_{\rm i}{\boldsymbol g} = 0,\; $$

where u is the bulk velocity, σ is the Cauchy stress tensor, ρ _i the density of ice and g the gravity vector.

For a constitutive relation, we adopt an unapproximated nonlinear extension of the orthotropic flow law (Rathmann and Lilien, 2022b), which, following plastic potential theory, has a nonlinear viscosity that depends on all orthotropic strain-rate tensor invariants:

(16)$$\matrix{ I_1 = {\dot{\boldsymbol \epsilon}}{\bi\, \colon \, }( {\bf m}_2\otimes{\bf m}_2-{\bf m}_3\otimes{\bf m}_3) /2,\; \cr I_2 = {\dot{\boldsymbol \epsilon}}{\bi\, \colon \, }( {\bf m}_3\otimes{\bf m}_3-{\bf m}_1\otimes{\bf m}_1) /2,\; \cr I_3 = {\dot{\boldsymbol \epsilon}}{\bi\, \colon \, }( {\bf m}_1\otimes{\bf m}_1-{\bf m}_2\otimes{\bf m}_2) /2,\; \cr I_4 = {\dot{\boldsymbol \epsilon}}{\bi\, \colon \, }( {\bf m}_2\otimes{\bf m}_3 + {\bf m}_3\otimes{\bf m}_2) /2,\; \cr I_5 = {\dot{\boldsymbol \epsilon}}{\bi\, \colon \, }( {\bf m}_3\otimes{\bf m}_1 + {\bf m}_1\otimes{\bf m}_3) /2,\; \cr I_6 = {\dot{\boldsymbol \epsilon}}{\bi\, \colon \, }( {\bf m}_1\otimes{\bf m}_2 + {\bf m}_2\otimes{\bf m}_1) /2,\; }$$

where the m_i's are the three reflection-symmetry directions that the fabric is presumed to have (taken to be coincident with the eigenvectors of a⁽²⁾).

Posed in inverse form (i.e. the deviatoric stress as a function of strain rate) the constitutive relation is

(17)$${\boldsymbol \tau} = \eta_0\sum_{i = 1}^3\bigg[{\eta_i}( I_{\,j_i}-I_{k_i}) {{\bf I}-3{\bf m}_i\otimes{\bf m}_i\over 2} + {\eta_{i + 3}}I_{i + 3}{{\bf m}_{j_i}\otimes{\bf m}_{k_i} + {\bf m}_{k_i}\otimes{\bf m}_{j_i}\over 2}\bigg],\; $$

where j = (2,  3,  1) and k = (3,  1,  2) are index variables, I is the identity and the nonlinear, isotropic part of the viscosity is

(18)$$\eta_0 = A^{-1/n}\left(\sum_{i = 1}^3\left[{\eta_i}( I_{\,j_i}-I_{k_i}) ^2 + {\eta_{i + 3}}I_{i + 3}^2\right]\right)^{(1-n)/2n},\; $$

where n is the flow law exponent (taken to be the canonical n = 3). Let

(19)$$\gamma = \sum_{i = 1}^3\left[2E_{\,j_i j_i}^{2/( n + 1) }E_{k_i k_i}^{2/( n + 1) }-E_{ii}^{4/( n + 1) }\right],\; $$

then the six dimensionless, relative directional viscosities are (for i = 1,  2,  3)

(20)$$\eta_i = {4\over 3\gamma}\left(E_{\,j_ij_i}^{2/( n + 1) } + E_{k_ik_i}^{2/( n + 1) }-E_{ii}^{2/n + 1) }\right),\; $$

(21)$$\eta_{i + 3} = 2E_{\,j_i k_i}^{-2/( n + 1) },\; $$

which depend on the eigenenhancements, E _ij, defined as the bulk directional enhancement factors (induced by fabric) in the directions of the fabric symmetry axes, m_i.

Directional enhancement factors

We are left to provide a mechanism to calculate the directional enhancement factors induced by an anisotropic fabric, defined as

(22)$$E_{ij}\equiv{{\bf m}_i {\bi\, \cdot \, } {\dot{\boldsymbol \epsilon}}( \hat{{\boldsymbol \tau}} ) {\bi\, \cdot \, }{\bf m}_j\over {\bf m}_i {\bi\, \cdot \, } {\dot{\boldsymbol \epsilon}}_{\rm iso}( \hat{{\boldsymbol \tau}} ) {\bi\, \cdot \, }{\bf m}_j}\ \quad {\rm for}\; i,\; j = 1,\; 2,\; 3,\; $$

where ${\dot {\boldsymbol \epsilon }}$ is the corresponding forward form of Eqn (17), ${\dot {\boldsymbol \epsilon }}_{\rm iso}$ is the strain-rate assuming isotropy (ϱ* = const.), and

(23)$$\hat{{\boldsymbol \tau}} = \tau_0\left\{\matrix{ {\bf I}/3-{\bf m}_i\otimes{\bf m}_i & \rm{if }\quad i = j\cr {\bf m}_i\otimes{\bf m}_j + {\bf m}_j\otimes{\bf m}_i & \rm{if }\quad i \neq j }\right.$$

are idealized longitudinal (i = j) and shear (i ≠ j) stress-tensor states (with some magnitude τ ₀) aligned with the fabric principal directions, m_i.

If we assume that grains do not interact, determining E _ij as a function of the fabric amounts to constructing a suitable grain-averaged rheology, over the grains that compose the polycrystal. Various interactionless averages have been used; end members are the Sachs hypothesis (constant stress) and the Taylor hypothesis (constant strain rate). Here, we follow Rathmann and Lilien (2022a) and take a linear combination of the enhancement factors found using the Sachs and Taylor hypotheses

(24)$$E_{ij}( {\dot{\boldsymbol \epsilon}}) = ( 1-\alpha_{{\rm rheo}}) E_{ij}^{{\rm Sachs}}( {\dot{\boldsymbol \epsilon}}) + \alpha_{{\rm rheo}}E^{{\rm Taylor}}_{ij}( {\dot{\boldsymbol \epsilon}}) ,\; $$

where α _rheo is an interaction parameter controlling the relative weights of the two hypotheses, and $E_{ij}^{{\rm Sachs}}( {\dot {\boldsymbol \epsilon }})$ and $E_{ij}^{{\rm Taylor}}( {\dot {\boldsymbol \epsilon }})$ are the enhancement factors assuming constant stress and strain rates, respectively. We calculate $E_{ij}^{{\rm Sachs}}( {\dot {\boldsymbol \epsilon }})$ and $E_{ij}^{{\rm Taylor}}( {\dot {\boldsymbol \epsilon }})$ assuming a linear, transversely isotropic grain rheology; details of the procedure can be found in Rathmann and Lilien (2022a), but note that the grain rheology depends on two grain parameters, E′_cc and E′_ca, which determine the compressional and shear enhancement of grains, respectively.

For the three required grain parameters, we take α _rheo = 0.0125, E′_cc = 1, and E′_ca = 10³ following Rathmann and Lilien (2022a), which produces a good match to shear enhancements found in laboratory deformation tests of approximately perfect single-maximum fabrics (Shoji and Langway, 1985; Pimienta and Duval, 1987). Such strong fabrics cannot be reproduced using L = 6, but require larger L due to regularization not being spectrally sharp (i.e. its effect is not concentrated exclusively at the largest wavenumber modes, l = L). As a result, the greatest possible shear enhancement is somewhat limited, although less so for compressional/extensional enhancements (Fig. 1, gray contours) – see Appendix A for details. We note that while this issue would be partly relieved by calculating the E _ij's using using a more realistic nonlinear, transversely isotropic grain rheology, this introduces dependencies on the eighth-order moments of ϱ* (Rathmann and others, 2021), prohibitively increasing computation time as it requires L ≥ 8.

Appendix B compares the nonlinear rheology of Rathmann and Lilien (2022b) (used here) to the nonlinear extension to the general orthotropic law of flow (GOLF) proposed by Martín and others (2009) at an ice divide. While Rathmann and Lilien (2022b) showed a relatively close match between these rheologies for uncoupled simulations, we use a coupled model of an ice divide to determine whether feedbacks between fabric and flow lead to larger differences in realistic settings. We find that both the nonlinear viscosity, η ₀, and the directional enhancement factors, E _ij, can differ markedly for fabrics produced in the coupled model, and that these differences grow as the fabric develops. These differences are most pronounced when the principal directions of the fabric are misaligned with the deformation axes, suggesting that the full nonlinear rheology used here is more appropriate for realistic settings, where complex bed topography inevitably creates misalignment between the fabric and deformation, and differences compound through the fabric–flow coupling. Moreover, the additional computational cost of the full rheology used here is negligible.

Heat flow

Heat flow is governed by (e.g. Zwinger and others, 2007; Hunter and others, 2021)

(25)$$\rho_{{\rm i}} c \left({\partial T\over \partial t}\, + {\bf u}\,{\boldsymbol \cdot}\,\nabla T\right) = \nabla\,{\boldsymbol \cdot}\,( k\nabla T) + \xi,\; $$

where constants c and k are the heat capacity and conductivity of ice, respectively, and $\xi = {\dot {\boldsymbol \epsilon }}\, \colon \, {\boldsymbol \sigma }$ is the rate of strain heating. Equation (25) is solved subject to the limitation that the ice does not exceed the pressure melting point following Zwinger and others (2007).

Free surface

For the top surface, we used a kinematic free surface boundary condition that adds another equation to be solved. Free surface evolution follows the usual 2-D problem

(26)$${\partial S\over \partial t} = \dot b + u_z - u_x{\partial S\over \partial x},\; $$

where S(x) is the surface elevation and $\dot b$ is the ice-equivalent accumulation rate.

Model across Dome C

As a first application of the model, we simulate a transect across Dome C, East Antarctica (Fig. 3). The transect follows a cross section acquired with pRES (Corr and others, 2021; Ershadi and others, 2022), allowing the simulated fabric field to be compared with radar-inferred fabrics in the top ~2000 m. The transect also crosses the EDC core site, which gives direct fabric measurements through ~3200 m depth (Durand and others, 2009). Model runs spanned 250 ka, which is longer than the characteristic timescale (thickness over accumulation) with as low as ~1.5 cm a⁻¹ accumulation and up to ~3400 m ice thickness. Even this length simulation may not be sufficient for the fabric to reach steady state, which could take 10 × the characteristic timescale (Martín and others, 2009). However, such a long integration time both exceeds the length of available forcing data and is computationally prohibitive; to mitigate some of this effect we initialize the model with a simple, non-isotropic fabric profile (see initial conditions below).

Figure 3. Map of Dome C area, showing model domain and location of validation data. Black line indicates model domain. Circles show pRES acquisition sites used in the text, with letter indicating panel of Figure 7 in which the results are plotted. The green star shows EDC core location. The colors show bed elevation from BedMachine v2 (Morlighem, 2020), while gray contours show surface elevation from REMA (Howat and others, 2019). Overview map shows location in Antarctica, with shading showing surface elevation from REMA.

Because there is significant lateral convergence and divergence along flowlines, we modify the flow equations to describe a 2.5-D flowband. In a coordinate system with x along flow and z vertical, Eqns (14) and (15) become

(27)$${\partial\sigma_{xx}\over \partial x} + {\partial\sigma_{xz}\over \partial z} + {\sigma_{xx}-\sigma_{yy}\over \omega( x) }{\partial \omega( x) \over \partial x} = 0$$

and

(28)$${\partial\sigma_{xz}\over \partial x} + {\partial\sigma_{zz}\over \partial z} + {\sigma_{xz}\over \omega( x) }{\partial \omega( x) \over \partial x} = \rho_{{\rm i}} g$$

for a flowband with width ω(x) (e.g. Hvidberg, 1996). The model equations were presented above, so we are left to describe the domain, and initial and boundary conditions.

Model domain

The model domain is a flowline extending 30 km in each direction from the EDC core site (Fig. 3). The domain runs through the pRES sites presented of Corr and others (2021), passing ~1 km from the EDC site, and approximately follows the surface gradient beyond the pRES locations. The size of the domain was chosen to prevent edge effects from impinging upon our results, and to ensure that the divide never migrated outside the model domain (which was not explicitly precluded by the boundary conditions). Surface elevation is from the REference Mosaic of Antarctica (REMA; Howat and others, 2019) and bed elevation is from BedMachine v2 (Morlighem, 2020; Morlighem and others, 2020). The model uses a triangular mesh with 500 m horizontal and 125 m vertical resolution, sufficient to capture large vertical gradients while keeping reasonable aspect ratios for the mesh elements (20 ka tests were run at double that resolution, and resulting fabrics were found to be indistinguishable from the resolution used here).

Initial conditions

We initialized the model to match present conditions at the EDC site. We used the EDC borehole temperatures (Buizert and others, 2021) and depth–age scale (Parrenin and others, 2007a) initially. For the fabric, we crudely approximated the depth profile (Durand and others, 2009) as a linear transition from the shallowest measured values in EDC (${\rm a}^{( 2) }_{zz}$=0.44), to a vertical single maximum at 2000 m depth (${\rm a}^{( 2) }_{zz}$=0.8), and constant from there to the bed.

Boundary conditions

At the surface, ice is assumed to match the shallowest observations from the EDC core (${\rm a}^{( 2) }_{zz} = 0.44$). We use time-varying boundary conditions for temperature and accumulation based on data from the EDC core (EPICA community members, 2004). Temperature and accumulation are assumed to vary spatially according to the modern patterns derived from reanalysis (Mottram and others, 2021) and temporally according to the deuterium excess in the core (Jouzel and others, 2007; Parrenin and others, 2007a) (Fig. 4). The surface temperature and accumulation at EDC therefore always match the values inferred from the core at a given time, while they vary in the rest of the domain. Spatially, the accumulation is scaled by anomaly relative to EDC while the temperature is offset by the anomaly to EDC. The age of the ice at the surface is set to zero. At the bed, we assume a spatially constant 55 mW m⁻² of geothermal heat flux and 0.2 mm a⁻¹ of basal melt (i.e. the bed-perpendicular velocity is set to a constant 0.2 mm a⁻¹) (Passalacqua and others, 2017). We assume that ice moves solely through internal deformation (so bed-parallel velocity is zero). The final boundary condition needed is some constraint on ice flow at the left and right outflows. At these boundaries, we assume that the pressure is approximately glaciostatic with the form Tr(σ)/3 = ρ gd′ where d′ is the depth relative to the modern surface; tying the pressure to the modern surface allows us to avoid specifying a velocity and prevents the divide from migrating outside the domain.

Figure 4. Forcing and modeled divide stability: (a) temperature history at EDC (Jouzel and others, 2007), (b) accumulation rate at EDC Parrenin and others, 2007a, (c) modeled divide position (positive northwestward, toward A′ in Fig. 3) and (d) ice thickness at the modeled divide position.

Results

We focus on the results of the model with rate factors calibrated from the EDC core for simplicity. Figure 5 shows the modeled velocity, age, and fabric on the transect across Dome C. The modeled fabric shows a typical divide profile; fabric generally strengthens toward a vertical single maximum with depth, with varying strength depending on the horizontal position. The single maximum is strongest in the first couple of kilometers to either side of the modeled divide. Farther away, there is weakening of the vertical component of the fabric. This weakening is associated with the tilting of the fabric away from vertical near the bed in areas with significant shear (Fig. 5c). The across-flow horizontal component of the fabric, $a^{( 2) }_{yy}$, is generally stronger than the along-flow component; however, the difference is small enough that by the classical test (the Woodcock parameter) it is still a vertical single maximum rather than a girdle near the divide.

Figure 5. Model output along Dome C transect. (a) Horizontal speeds. Contours show ages of the ice. (b) Vertical component of the fabric. Vertical, red lines show locations of pRES acquisitions plotted in Figure 7; these data are a representative subset from Ershadi and others (2022). (c) Rotation of the fabric from vertical (degrees counter clockwise). In all panels, dotted black line shows the modeled, modern divide position.

The modeled horizontal ice-flow speed is very slow, reaching only 0.24 m a⁻¹ at the edges of the model domain. These speeds agree well with GPS measurements of velocity in the area, which found speeds up to 0.21 m a⁻¹ 25 km from the divide (Vittuari and others, 2004). However, those GPS measurements were not directly on the transect considered here, so more quantitative evaluation of the velocities is precluded.

The modeled ice divide was extremely stable in both position and thickness (Fig. 4). The ice thickness varied by ~10 m through the simulation, significantly <~250 m inferred from models that were specifically targeted at reconstructing surface elevations (Parrenin and others, 2007b), while the divide position moved by only ~300 m (which is at the limit of the horizontal mesh resolution). The smaller variation in divide thickness is likely a result of the boundary conditions; regional ice thickness changes were not captured by the pseudo-pressure conditions used, but actual ice-thicknesses changes during the last 250 ka likely affected the entire continent. The lack of thickness change is thus a consequence of using a model domain with a finite width smaller than that of the continent; a model with no horizontal dimension can evolve freely as it is subject to only local conditions (e.g. Parrenin and others, 2007b), and a continental-scale model can thicken or thin as a whole, but our limited domain depends on boundary conditions which do not permit such changes.

The modeled age profile shows a very minor double Raymond bump (maximum bump height 60 m, or 2% of the total thickness). Despite the extremely long characteristic timescale of the divide, the bump is allowed to form since model divide is so stable throughout the simulation. The double bump is characteristic of directional hardening due to fabric development (Martín and others, 2009), and has been observed over domes with shorter characteristic timescales (Goel and others, 2020). Though there is not a Raymond bump in radar data in the area (Cavitte and others, 2021), given the small amplitude of the modeled bump, and the overly stable modeled divide, this difference is unsurprising. Despite the relatively stable modeled divide, the bump amplitude is much smaller than the 10% of the ice thickness that would be expected in steady state (Goel and others, 2020).

Model/data comparison

Ice-core measurements

We evaluate the modeled fabric at the modeled divide position, which is ~1 km southeast of the true position in all simulations (Fig. 6). Comparison at the closest approach to EDC yields similar results, but any effect of nonlinearity of the viscosity (i.e. the Raymond bump) is expected at the modeled rather than true divide, so the former provides the more fair comparison. At both the point of closest approach and at the modeled divide, the modeled fabric shows a weaker vertical component (λ ₃) and one stronger horizontal component (λ ₂) than was measured in the EDC core (Durand and others, 2009). This difference is seen in the model both with laboratory and ice-core-calibrated rate factors, though it is significantly lessened for the ice-core-calibrated values. The model-data mismatch could be attributable to imperfect model physics describing the fabric development (due either to missing processes or poor parameterization), insufficiency of the simplified 2.5-D model to capture the stress state, temporal changes not considered by the model, or some combination of these (fabric regularization does not affect λ _i much even for L = 6; see Appendix A). Diverging flowlines lead to a more dome-like stress state, and the flow direction is not well constrained right at the divide, so the assumptions of the 2.5-D model might break down in the center of the model domain (Passalacqua and others, 2016). Temporal changes in the large-scale flow could significantly alter our results; simulations extend back 250 ka, over which ice thickness is thought to have changed by >250 m (Parrenin and others, 2007b), and accumulation patterns have likely changed even over the last 260 years, which suggests that migration of the divide is possible (Urbini and others, 2008). While we are able to incorporate local thickness and accumulation changes, such temporal changes are likely to have affected the large-scale flow of the ice sheet, and cannot be easily captured by a local model such as that used here. To further distinguish among possible causes of the model-data misfit, we compare our fabric results with pRES data along the model domain.

Figure 6. Comparison of modeled fabric (Elmer/Ice) and measured fabric from the EDC core (Durand and others, 2009). Data are shown as squares. Colors indicate eigenvalue number (blue for λ ₃, orange for λ ₂ and green for λ ₁). Model output is shown as lines, for both laboratory (solid) and ice-core (dashed) calibrated rate factors, with colors corresponding to the data.

Phase-sensitive radar

pRES returns are affected by ice-crystal fabric because the dielectric permittivity is different parallel and perpendicular to the c-axis of a single ice crystal (Fujita and others, 2006). In a typical radar geometry, returns are insensitive to the vertical component of the fabric but their strength and phase with azimuth is related to the horizontal anisotropy (e.g. Rathmann and others, 2022). From these returns, it is thus possible, with minimal assumptions, to infer the difference in the horizontal eigenvalues of the fabric distribution (generally assumed to be the ODF). Recently published pRES measurements on a transect of Dome C (Corr and others, 2021; Ershadi and others, 2022) therefore provide a way to test our modeled fabric field against observations, away from the ice-core site. Ershadi and others (2022) used these data to infer the difference between the horizontal eigenvalues and the orientation of the larger eigenvector with depth along the transect modeled here. Since our model does not capture rotation of the fabric out of the model plane, the modeled orientation of the largest eigenvalue is always transverse to flow; this agrees with all pRES measurements deeper than 500 m within error (shallower orientations are highly uncertain in the pRES due to weaker returns). We therefore compare only the difference in horizontal eigenvalues inferred from pRES and in the model output. We do so at seven of the sites (comparison at the remaining sites from Ershadi and others (2022) is shown as Supplementary Fig. S2).

Figure 7 shows the modeled and radar-inferred horizontal anisotropy at the seven sites marked in Figures 3 and 5. At the location of EDC, ice-core data are plotted as well (Durand and others, 2009). Both the ice-core- and laboratory-calibrated models match the data well near the surface, but systematically find too-large horizontal anisotropy between 500 and 2000 m throughout the domain. Simulations with ice-core- and laboratory-calibrated recrystallization produce very similar horizontal anisotropy (negligible differences except deeper than 3000 m) despite large differences in the vertical component of the fabric. Differences compared to the pRES-inferred values reach a factor of 2.5 for both simulations. The pRES-inferred values match the EDC data well overall, although they generally lie at the low end of the spread in the ice-core data (Fig. 7d). The ice-core data thus suggest that some of the difference between the pRES-inferred anisotropy and the modeled anisotropy may be due to the pRES method slightly underestimating the eigenvalue difference. The majority of the model-pRES difference is therefore attributable to the model producing fabrics with too much horizontal anisotropy rather than errors in the pRES measurements.

Figure 7. Horizontal eigenvalue difference of modeled and pRES-inferred fabric (Ershadi and others, 2022) at locations shown in Figures 3 and 5. (a–g) are E18, E12, E3, EPICA, W6, W12 and W18 from Ershadi and others (2022). pRES is shown with black circles, ice-core-calibrated model with dark gray lines, laboratory-calibrated model with dashed, light gray lines and the EDC core with red (only where the pRES and model coincide with the core). Remaining sites from Ershadi and others (2022) are shown in Supplementary Figure S2.

Discussion

Implications of model-data misfit

Using ice-core- or laboratory-calibrated rate factors, the model produces slightly too strong horizontal anisotropy throughout the model domain. Given this difference and the available data from EDC, it is reasonable to infer that the vertical component of the modeled fabric is too weak throughout the domain as well. Although we are unable to attribute the model-data misfit to a single cause, it seems likely that the differences are, in large part, due to limitations of the fabric model rather than of the ice-flow model, time-dependent forcing or boundary conditions. Since the fabric model considers the processes thought to be most relevant for fabric development, and the model can qualitatively reproduce fabrics found in ice cores and laboratory deformation experiments, we consider assumptions about the rates of these processes to be the most likely cause of the model-data misfit.

Although the rate-factor parameterizations (Γ₀ and Λ₀) of the different fabric processes are only approximations of the true physics of fabric development, they are qualitatively able to produce the diversity of natural and synthetic fabrics (Richards and others, 2021). We consider it likely that the rate factors are not accurately described by either set of coefficients (lab- or ice-core-calibrated coefficients), and that further calibration against experimental data could yield a closer match to the data, without large changes to the models of fabric processes (i.e. the functional form of Eqn (5)). Indeed, laboratory-calibrated rates were effectively tuned to three data points without uncertainty bounds (Richards and others, 2021). Additional laboratory measurements, designed to allow the rate factors to be tuned against data (at different strain rates, temperatures and stresses), might lead to better constrained coefficient functions of the recrystallization models, and might clarify whether the model-data mismatch is due to sparse experimental data.

Similarly, calibrating the coefficient functions by comparing data from EDC to a one-dimensional, time-varying model with a better-constrained strain tensor, might change the ice-core-calibrated values and elucidate whether the model is capable of exactly reproducing the data. For example, it is possible that further data will clarify whether the rate of migration recrystallization is truly linear in the strain rate or whether a more complicated relationship exists; given the difference in strain rates in the laboratory and natural settings, a rate factor that depends non-linearly on the strain rate may reconcile our estimated coefficients.

Calibration to ice cores drilled in more complex flow regimes could provide additional information about this difference. Fabric measurements on the EastGRIP core (Westhoff and others, 2021) could be particularly valuable for calibration, since it samples an active ice stream moving 50 m a⁻¹. Measurements on NEEM (Montagnat and others, 2014), drilled on a divide rather than a dome, and the South Pole ice core (Voigt, 2017), drilled on a flank moving 10 m a⁻¹ would add two more stress states for calibration. Multi-site calibration requires a more efficient approach than the brute-force method used here, and the more complex stress states may not be accurately simulated with the simple Nye and Dansgaard–Johnsen models used here. Nevertheless, such multi-site calibration may help elucidate whether site-specific characteristics of Dome C, for example the cold temperature and low strain rates, explain some of the difference between our ice-core- and laboratory-calibrated rate factors.

Large-scale modeling with a spectral fabric description

There are important differences between the spectral and tensorial representations both in terms of physics and in terms of technical implementation (Richards and others, 2021; Rathmann and others, 2021). A clear advantage of the spectral formulation is that it does not need a closure approximation, while consistency requires at least a minimal closure for a tensorial representation. Definitionally, each high-order structure tensor must have a number of non-zero components leading to the sub-traces taking values compatible with the lower-order structure tensors (e.g. Advani and Tucker, 1987). This becomes problematic when the evolution of lower-order tensors depends on the higher-order ones, since the higher-order contributions cannot ever be assumed to be zero. In the spectral representation, harmonics with l > 0 represent anomalies relative to isotropy, modes >L can safely be ignored for fabric development (an equivalent assumption can be made for the tensorial representation, but since for isotropic fabrics the higher-order tensors are non-zero it is more complicated than assuming spectral components are zero).

Perhaps the most important physical disadvantage to the spectral approach is the need to balance regularization (in orientation space) with the maximum fabric strength that is possible. For example, using canonical Laplacian regularization with L = 6, λ ₃ saturates ~0.8 in unconfined compression, while λ ₃ → 1 is expected under sufficient strain; we avoided this problem with carefully calibrated hyper diffusion that permitted λ ₃ > 0.99 while keeping the model stable. Such calibration must, however, be adjusted for different values of L, or else spurious oscillations or a maximum fabric strength will result.

In addition to physical differences, there are a number of technical differences between the spectral and tensorial representations that affect the feasibility of using the two representations in a large-scale model. We briefly mention the two that we see as most important. At the same order of truncation, the spectral and tensorial representations have the same number of degrees of freedom. With no simplifications, the spectral representation always tracks two times the number of complex, independent components (plus the zero-order component), while the tensorial representation, in native form, tracks a number of degenerate components that increases proportional to 3^m. For example, a⁽⁴⁾ has 81 components, but only nine are unique (assuming a⁽²⁾ is known); in order to avoid the computational cost of tracking all 81, however, cumbersome expressions must be derived for each individual independent component. This rapidly becomes infeasible if yet higher-order structures are desired. Even if only a⁽²⁾ and a⁽⁴⁾ were explicitly tracked, modeling the evolution of a⁽⁴⁾ would require an approximation of the non-zero elements of a⁽⁶⁾, which has only 15 independent components out of its 729 elements. The tedious work of determining these relations is avoided with the spectral representation; the simplifications of Eqn (13) remain the same at each L, and can be handled automatically to reduce the problem to only the independent components. Thus, the lack of degeneracy allows us to increase the order of the fabric description up to any desired order, so future studies with greater computational power and some need for highly resolved representations of the fabric could use arbitrarily large L (we tested up to L = 12, but expect larger L to be feasible with greater computational resources).

Another technical difference is that some models of fabric process (e.g. DDRX) naturally involve a^(m), whereas others (e.g. lattice rotation and CDRX) involve ϱ*, which might influence the choice of representation. Apart from having to construct structure tensors when comparing modeled fabrics to ice-core measurements, constructing a⁽²⁾ and a⁽⁴⁾ – even when DDRX is negligible – is necessary in the spectral approach since the eigenenhancements, E _ij, (or fluidity tensor) depend on both a⁽²⁾ and a⁽⁴⁾. Reconstructing a⁽²⁾ and a⁽⁴⁾ is quick compared to the time taken to model the evolution of ϱ*(x,  t) in 2-D, so does not add significant computational cost. Relatedly, all a^(m) lie between 0 and 1 by definition, but may fall outside this range due to numerical or truncation errors. Enforcing these bounds by renormalizing a^(m) (e.g. dividing through by the sum of unnormalized eigenvalues and setting negative eigenvalues to zero) becomes increasingly tedious for structure tensors of increasing order m > 4, as determining eigenvalues becomes increasingly computationally difficult. In the spectral approach, we found it useful to cap the angular power spectrum of $\hat \varrho _l^m$ to ensure the largest eigenvalue is ≤1 for all a^(m); to do so, we renormalize the power spectrum to that of the delta function if needed (representing a perfect single maximum). Because of the way the harmonic modes describe an anomaly from isotropy, this in turn also ensures the smallest eigenvalue is ≥0. While this does not ensure the smallest eigenvalues are ≥0 for relatively strong girdle fabrics, such fabrics do not arise in the present simulations. If renormalization is not carefully addressed, this can lead to numerical issues – particularly for strong single-maximum fabrics – caused by spectral coefficients resulting in negative directional viscosities (see Appendix A).

Outlook

The spectral representation of fabric provides a convenient way to include additional processes in modeling of fabric evolution, and thus provides a path forward for integrating future small-scale model development into large-scale ice-flow models. Simply moving to L > 2 permitted us to model the effect of DDRX, and further increases in the order of approximation could allow for use of nonlinear grain rheologies in deriving the bulk anisotropic fluidity; such a transition is greatly eased by the spectral representation of the fabric. The contrast between fabric modeled with laboratory- and ice-core-calibrated rate factors suggests a disconnect between how fabric develops in the laboratory and in natural settings; presumably this difference stems from the 5–7 orders of magnitude larger strain rates at which laboratory tests are conducted compared to typical ice-core settings. Some fundamental differences have been clear before this work; for example, laboratory ice-deformation experiments do not seem to produce strong girdles, despite widespread observations of girdles in natural ice. This modeling, and the recent work on which it relies (Richards and others, 2021; Rathmann and Lilien, 2022b), allows direct comparison of the fabric evolution at laboratory scales/rates and ice-core scales/rates, potentially providing a path forward for understanding these differences.

In this work we have avoided areas with rapid ice flow, since data with which to validate the model are unavailable in such areas. However, accurately capturing rheology in areas of fast flow is critical for accurate modeling of outlet glaciers and ice streams. For example, models that do not include anisotropy sometimes find that ad-hoc weakening of shear margins (beyond the weakening expected from temperature alone) produces a better match to velocity observations (e.g. Joughin and others, 2012). Such weakening is often attributed to anisotropy (Minchew and others, 2018; Grinsted and others, 2022). The complexity of flow in such areas necessitates a relatively complex fabric model such as that used here. Shear margins have relatively short residence times, so fabric cannot be assumed to be in steady state with its in situ stress state, precluding the use of models that make such an assumption (Graham and others, 2018), yet recrystallization is thought to be active at the relatively high temperatures and high strain rates found in shear margins (Faria and others, 2014; Hunter and others, 2021), precluding a⁽²⁾-only approaches. Moreover, the multiple deformation modes active in shear margins (transition from vertical shear to plug flow coincident with horizontal shear) suggest a need to use a tensorially anisotropic viscosity. Fabric will evolve through time, and will do so differently from other potential modes of weakening, such as crevassing, so the need to capture the full fabric field and its evolution is compounded in transient simulations. The model presented here could be applied to areas with multiple active deformation modes, such as shear margins, to more accurately capture the variations in directional viscosity and the evolution of flow through time, an application for which existing anisotropic models are not calibrated. Further calibration of this model to use in such areas could take advantage of additional geophysical constraints, such as seismics (e.g. Lutz and others, 2022), in order to validate modeled fabric in these areas.

Conclusions

We incorporated a recently developed model of fabric evolution into a large-scale ice-flow model. The model includes lattice rotation, rotation recrystallization and migration recrystallization, the three processes thought to be most important for fabric development, and couples the fabric to ice flow through a rheology with fewer approximations than those used previously. We applied this new, coupled model to simulate ice flow and fabric development on a transect across Dome C, East Antarctica, where relatively plentiful validation data exist. At the small scale, previous work showed that this model does an excellent job reproducing laboratory fabric-development experiments (Richards and others, 2021), but applying it to the large scale, as we have done here, reveals that the model is less able to reproduce the large-scale fabric found in the EDC core and inferred from pRES. Calibrating rate factors to match the EDC core, using a zero-dimensional model, reduced, but did not eliminate, this misfit. This discrepancy cannot be firmly attributed to one source, but insufficient constraints on the relative rates of different processes (migration recrystallization, rotation recrystallization and lattice rotation) may be the cause. The difference found between calibrated process rate factors using laboratory and ice-core data suggests a gap exists in our understanding of how ice-crystal fabric develops differently in laboratory compared to natural (in situ) conditions. However, this type of modeling allows quantitative comparison across spatial and temporal scales, and thus provides a potential path forward for reconciling the differences of fabric development found in laboratory and natural conditions. The model is one of the first that is well-suited to simulating the effect of fabric on the complex deformation experienced in areas such as shear margins, and may thus allow more accurate simulation of such areas.