# A physically-based formulation for texture evolution during dynamic recrystallization. A case study of ice

Thomas Chauve, Maurine Montagnat, Véronique Dansereau, Pierre Saramito,
Kévin Fourteau, Andréa Tommasi

► To cite this version:

Thomas Chauve, Maurine Montagnat, Véronique Dansereau, Pierre Saramito, Kévin Fourteau, et al. A physically-based formulation for texture evolution during dynamic recrystallization. A case study of ice. Comptes Rendus. Mécanique, 2024, 352 (G1), pp.99-134. 10.5802/crmeca.243. hal-04578817

HAL Id: hal-04578817

https://hal.science/hal-04578817

Submitted on 17 May 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Comptes Rendus
***
Mécanique

Thomas Chauve, Maurine Montagnat, Véronique Dansereau, Pierre Saramito,
Kévin Fourteau and Andréa Tommasi

A physically-based formulation for texture evolution during dynamic
recrystallization. A case study of ice

Volume 1, issue 0 (0000), p. 000-000

https://doi.org/10.5802/crmeca.243

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Comptes Rendus
***
Mécanique

Thomas Chauve, Maurine Montagnat, Véronique Dansereau, Pierre Saramito,
Kévin Fourteau and Andréa Tommasi

A physically-based formulation for texture evolution during dynamic
recrystallization. A case study of ice

Volume 1, issue 0 (0000), p. 000-000

https://doi.org/10.5802/crmeca.243

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Research article / *Article de recherche*

# A physically-based formulation for texture evolution during dynamic recrystallization. A case study of ice

*Une formulation basée sur la physique pour simuler l'évolution de la texture pendant la recristallisation dynamique. Une étude de cas pour la glace*

Thomas Chauve<sup>*, a</sup>, Maurine Montagnat<sup>a, b</sup>, Véronique Dansereau<sup>c</sup>, Pierre Saramito<sup>d</sup>, Kévin Fourteau<sup>e</sup> and Andréa Tommasi<sup>f</sup>

<sup>a</sup> Univ. Grenoble Alpes, CNRS, IRD, G-INP, IGE, Grenoble, France.

<sup>b</sup> Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Études de la Neige, Grenoble, France.

<sup>c</sup> Institut des Sciences de la Terre, Université Grenoble Alpes, CNRS (5275), Grenoble, France.

<sup>d</sup> Lab. Jean Kuntzmann, CNRS, Université Grenoble-Alpes, F-38041 Grenoble, France.

<sup>e</sup> Géosciences Montpellier - CNRS, Université de Montpellier, France.

*E-mails*: thomas.chauve@univ-grenoble-alpes.fr (T. Chauve), maurine.montagnat@univ-grenoble-alpes.fr (M. Montagnat), veronique.dansereau@univ-grenoble-alpes.fr (V. Dansereau), pierre.saramito@imag.fr (P. Saramito), kevin.fourteau@meteo.fr (K. Fourteau), andrea.tommasi@umontpellier.fr (A. Tommasi)

**Abstract.** Dynamic recrystallization can have a strong impact on texture development during the deformation of polycrystalline materials at high temperatures, particularly for materials with strong viscoplastic anisotropy such as ice. Owing to this anisotropy, recrystallization is essential for ensuring strain compatibility, and the development of textures leads to anisotropic softening. Accurate prediction of the effect of recrystallization on the texture evolution of ice is therefore crucial to adequately account for texture-induced mechanical anisotropy in large-scale models of glacial ice flow. However, this prediction remains a challenge.

We propose a new formulation for modeling texture evolution due to dynamic recrystallization on the basis of observations of the evolution of the microstructure and texture of ice deforming by dislocation creep and dynamic recrystallization. This formulation relies on an orientation attractor that maximizes the resolved shear stress on the easiest slip system in the crystal. It is implemented in the equation describing the evolution of the crystal orientation with deformation and is coupled with an anisotropic viscoplastic law that provides the mechanical response of the ice crystal. This set of equations, which is the core of the R³iCe model is solved by a finite-element method with a semi-implicit scheme coded using the Rheolef library. The resulting open-source software R³iCe is validated by comparison with laboratory creep data for ice polycrystals under uniaxial compression, simple shear, and uniaxial tension. It correctly reproduces the texture evolution and mechanical softening observed in the experiment during tertiary creep. Although the

<footnote>* Corresponding author.</footnote>
Research article / *Article de recherche*

# A physically-based formulation for texture evolution during dynamic recrystallization. A case study of ice

Une formulation basée sur la physique pour simuler l'évolution de la texture pendant la recristallisation dynamique. Une étude de cas pour la glace

Thomas Chauve<sup>*, a</sup>, Maurine Montagnat<sup>a, b</sup>, Véronique Dansereau<sup>c</sup>, Pierre Saramito<sup>d</sup>, Kévin Fourteau<sup>e</sup> and Andréa Tommasi<sup>f</sup>

<sup>a</sup> Univ. Grenoble Alpes, CNRS, IRD, G-INP, IGE, Grenoble, France.

<sup>b</sup> Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Études de la Neige, Grenoble, France.

<sup>c</sup> Institut des Sciences de la Terre, Université Grenoble Alpes, CNRS (5275), Grenoble, France.

<sup>d</sup> Lab. Jean Kuntzmann, CNRS, Université Grenoble-Alpes, F-38041 Grenoble, France.

<sup>e</sup> Géosciences Montpellier - CNRS, Université de Montpellier, France.

*E-mails*: thomas.chauve@univ-grenoble-alpes.fr (T. Chauve), maurine.montagnat@univ-grenoble-alpes.fr (M. Montagnat), veronique.dansereau@univ-grenoble-alpes.fr (V. Dansereau), pierre.saramito@imag.fr (P. Saramito), kevin.fourteau@meteo.fr (K. Fourteau), andrea.tommasi@umontpellier.fr (A. Tommasi)

**Abstract.** Dynamic recrystallization can have a strong impact on texture development during the deformation of polycrystalline materials at high temperatures, particularly for materials with strong viscoplastic anisotropy such as ice. Owing to this anisotropy, recrystallization is essential for ensuring strain compatibility, and the development of textures leads to anisotropic softening. Accurate prediction of the effect of recrystallization on the texture evolution of ice is therefore crucial to adequately account for texture-induced mechanical anisotropy in large-scale models of glacial ice flow. However, this prediction remains a challenge.

We propose a new formulation for modeling texture evolution due to dynamic recrystallization on the basis of observations of the evolution of the microstructure and texture of ice deforming by dislocation creep and dynamic recrystallization. This formulation relies on an orientation attractor that maximizes the resolved shear stress on the easiest slip system in the crystal. It is implemented in the equation describing the evolution of the crystal orientation with deformation and is coupled with an anisotropic viscoplastic law that provides the mechanical response of the ice crystal. This set of equations, which is the core of the R³iCe model is solved by a finite-element method with a semi-implicit scheme coded using the Rheolef library. The resulting open-source software R³iCe is validated by comparison with laboratory creep data for ice polycrystals under uniaxial compression, simple shear, and uniaxial tension. It correctly reproduces the texture evolution and mechanical softening observed in the experiment during tertiary creep. Although the

<footnote>* Corresponding author.</footnote>
present formulation is too time-consuming for direct implementation in large-scale ice flow models, R³iCe can be used to adjust the parameterization used to implement texture-induced anisotropy in these models. The validation was performed for ice, but the R³iCe implementation is generic and applies to any material whose behavior may be adequately described using an anisotropic flow law.

**Résumé.** La recristallisation dynamique peut avoir un impact important sur le développement des textures pendant la déformation des matériaux polycristallins à haute température, en particulier pour les matériaux présentant une forte anisotropie viscoplastique comme la glace. En raison de cette anisotropie, la recristallisation est essentielle pour assurer la compatibilité des déformations et le développement des textures conduit à un adoucissement anisotrope. Une prédiction précise de l'effet de la recristallisation sur l'évolution de la texture de la glace est donc cruciale pour tenir compte de l'anisotropie mécanique induite par la texture dans les modèles à grande échelle de l'écoulement de la glace. Pourtant, cette prédiction reste un défi.

Nous proposons une nouvelle formulation pour modéliser l'évolution de la texture due à la recristallisation dynamique, basée sur des observations de l'évolution de la microstructure et de la texture de la glace se déformant par fluage de dislocations et recristallisation dynamique. Cette formulation repose sur un attrac-teur d'orientation qui maximise la contrainte de cisaillement résolue sur le système de glissement le plus fa-cile dans le cristal. Elle est mise en œuvre dans l'équation décrivant l'évolution de l'orientation du cristal avec la déformation et couplée à une loi viscoplastique anisotrope qui fournit la réponse mécanique du cristal de glace. Cet ensemble d'équations, qui constitue le cœur du modèle R3iCe, est résolu par une méthode d'élé-ments finis avec un schéma semi-implicite codé à l'aide de la bibliothèque Rheolef. Le logiciel libre R3iCe qui en résulte est validé par comparaison avec des données de fluage en laboratoire pour des polycristaux de glace soumis à une compression uniaxiale, à un cisaillement simple et à une tension uniaxiale. Il reproduit correctement l'évolution de la texture et l'adoucissement mécanique observés lors de l'expérience pendant le fluage tertiaire. Bien que la formulation actuelle soit trop coûteuse en temps pour une mise en œuvre directe dans des modèles d'écoulement de glace à grande échelle, R3iCe peut être utilisé pour ajuster la paramétri-sation utilisée pour mettre en œuvre l'anisotropie induite par la texture dans ces modèles. La validation a été effectuée pour la glace, mais l'implémentation de R3iCe est générique et s'applique à tout matériau dont le comportement peut être décrit de manière adéquate à l'aide d'une loi d'écoulement anisotrope.

**Keywords.** dynamic recrystallization, texture, viscoplastic anisotropy, finite-element method, ice.

**Mots-clés.** recristallisation dynamique, texture, anisotropie viscoplastique, méthode des éléments finis, glace.

**Funding.** This work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation programme [grant agreement No 882450 - ERC RhEoVOLUTION]. Kevin Fourteau received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (IVORI, grant no. 949516).

Manuscript received 6 October 2023, revised 22 January 2024, accepted 19 February 2024.

**Notation**

<table><caption>**Table 1.** Mathematical notations</caption><tr><td><b>a</b></td><td>Scalar</td><td></td></tr><tr><td><b>a</b></td><td>Vector</td><td></td></tr><tr><td><b>A</b></td><td>Tensor</td><td></td></tr><tr><td><b>A</b><sup>D</sup></td><td>Deviatoric part of <b>A</b></td><td><b>A</b><sup>D</sup> = <b>A</b> - <sup>1</sup>&frasl;<sub>dim(<b>A</b>)</sub> Tr(<b>A</b>)</td></tr><tr><td>&otimes;</td><td>Tensorial (dyadic) product</td><td><b>A</b> = <b>b</b> &otimes; <b>c</b></td></tr><tr><td>Tr()</td><td>Trace of a Tensor</td><td>Tr(<b>A</b>)</td></tr><tr><td>dim()</td><td>dimension of a tensor</td><td>dim(<b>A</b>)</td></tr><tr><td>&nabla;</td><td>&nabla; operator</td><td>(<sup>&part;</sup>&frasl;<sub>&part;x</sub>, <sup>&part;</sup>&frasl;<sub>&part;y</sub>, <sup>&part;</sup>&frasl;<sub>&part;z</sub>)</td></tr><tr><td>.</td><td>Scalar product</td><td><b>a</b> &bull; <b>b</b></td></tr><tr><td>:</td><td>Double contracted product</td><td><b>A</b> : <b>B</b></td></tr><tr><td>||&middot;||</td><td>Norm</td><td>||<b>a</b>||</td></tr><tr><td><b>A</b><sup>n</sup></td><td>matrix product <i>n</i> times</td><td><b>A</b> &times; <b>A</b> &middot;&middot;&middot; &times; <b>A</b></td></tr></table>
present formulation is too time-consuming for direct implementation in large-scale ice flow models, R³iCe can be used to adjust the parameterization used to implement texture-induced anisotropy in these models. The validation was performed for ice, but the R³iCe implementation is generic and applies to any material whose behavior may be adequately described using an anisotropic flow law.

**Résumé.** La recristallisation dynamique peut avoir un impact important sur le développement des textures pendant la déformation des matériaux polycristallins à haute température, en particulier pour les matériaux présentant une forte anisotropie viscoplastique comme la glace. En raison de cette anisotropie, la recristallisation est essentielle pour assurer la compatibilité des déformations et le développement des textures conduit à un adoucissement anisotrope. Une prédiction précise de l'effet de la recristallisation sur l'évolution de la texture de la glace est donc cruciale pour tenir compte de l'anisotropie mécanique induite par la texture dans les modèles à grande échelle de l'écoulement de la glace. Pourtant, cette prédiction reste un défi.

Nous proposons une nouvelle formulation pour modéliser l'évolution de la texture due à la recristallisation dynamique, basée sur des observations de l'évolution de la microstructure et de la texture de la glace se déformant par fluage de dislocations et recristallisation dynamique. Cette formulation repose sur un attrac-teur d'orientation qui maximise la contrainte de cisaillement résolue sur le système de glissement le plus fa-cile dans le cristal. Elle est mise en œuvre dans l'équation décrivant l'évolution de l'orientation du cristal avec la déformation et couplée à une loi viscoplastique anisotrope qui fournit la réponse mécanique du cristal de glace. Cet ensemble d'équations, qui constitue le cœur du modèle R3iCe, est résolu par une méthode d'élé-ments finis avec un schéma semi-implicite codé à l'aide de la bibliothèque Rheolef. Le logiciel libre R3iCe qui en résulte est validé par comparaison avec des données de fluage en laboratoire pour des polycristaux de glace soumis à une compression uniaxiale, à un cisaillement simple et à une tension uniaxiale. Il reproduit correctement l'évolution de la texture et l'adoucissement mécanique observés lors de l'expérience pendant le fluage tertiaire. Bien que la formulation actuelle soit trop coûteuse en temps pour une mise en œuvre directe dans des modèles d'écoulement de glace à grande échelle, R3iCe peut être utilisé pour ajuster la paramétri-sation utilisée pour mettre en œuvre l'anisotropie induite par la texture dans ces modèles. La validation a été effectuée pour la glace, mais l'implémentation de R3iCe est générique et s'applique à tout matériau dont le comportement peut être décrit de manière adéquate à l'aide d'une loi d'écoulement anisotrope.

**Keywords.** dynamic recrystallization, texture, viscoplastic anisotropy, finite-element method, ice.

**Mots-clés.** recristallisation dynamique, texture, anisotropie viscoplastique, méthode des éléments finis, glace.

**Funding.** This work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation programme [grant agreement No 882450 - ERC RhEoVOLUTION]. Kevin Fourteau received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (IVORI, grant no. 949516).

Manuscript received 6 October 2023, revised 22 January 2024, accepted 19 February 2024.

**Notation**

<table><caption>**Table 1.** Mathematical notations</caption><tr><td><b>a</b></td><td>Scalar</td><td></td><td></td></tr><tr><td><b>a</b></td><td>Vector</td><td></td><td></td></tr><tr><td><b>A</b></td><td>Tensor</td><td></td><td></td></tr><tr><td><b>A</b><sup>D</sup></td><td>Deviatoric part of <b>A</b></td><td><b>A</b><sup>D</sup> = <b>A</b> - <sup>1</sup>&frasl;<sub>dim(<b>A</b>)</sub> Tr(<b>A</b>)</td><td><b>A</b><sub>ij</sub><sup>D</sup> = <b>A</b><sub>ij</sub> - <sup>&delta;<sub>ij</sub></sup>&frasl;<sub>dim(<b>A</b>)</sub> Tr(<b>A</b>)</td></tr><tr><td>&otimes;</td><td>Tensorial (dyadic) product</td><td><b>A</b> = <b>b</b> &otimes; <b>c</b></td><td><b>A</b><sub>ij</sub> = <b>b</b><sub>i</sub> <b>c</b><sub>j</sub></td></tr><tr><td>Tr()</td><td>Trace of a Tensor</td><td>Tr(<b>A</b>)</td><td><b>A</b><sub>ii</sub></td></tr><tr><td>dim()</td><td>dimension of a tensor</td><td>dim(<b>A</b>)</td><td></td></tr><tr><td>&nabla;</td><td>&nabla; operator</td><td>(<sup>&part;</sup>&frasl;<sub>&part;x</sub>, <sup>&part;</sup>&frasl;<sub>&part;y</sub>, <sup>&part;</sup>&frasl;<sub>&part;z</sub>)</td><td><sup>&part;</sup>&frasl;<sub>&part;x<sub>i</sub></sub></td></tr><tr><td>.</td><td>Scalar product</td><td><b>a</b> &bull; <b>b</b></td><td><b>a</b><sub>i</sub> <b>b</b><sub>i</sub></td></tr><tr><td>:</td><td>Double contracted product</td><td><b>A</b> : <b>B</b></td><td><b>A</b><sub>ij</sub> <b>B</b><sub>ij</sub></td></tr><tr><td>||.||</td><td>Norm</td><td>||<b>a</b>||</td><td><b>a</b><sub>i</sub><sup>2</sup></td></tr><tr><td><b>A</b><sup>n</sup></td><td>matrix product <i>n</i> times</td><td><b>A</b> &times; <b>A</b> &middot;&middot;&middot; &times; <b>A</b></td><td></td></tr></table>
# 1. Introduction

Dynamic recrystallization (DRX) is a mechanism by which local stress incompatibilities and local strain heterogeneities are relaxed [1–3]. By *local*, we mean the scale of the grains that constitute crystalline materials such as ice, metals, and most rocks. Dynamic recrystallization occurs in plastically deformed materials at high temperatures. It is characterized by the nucleation of new grains and grain boundary migration and can lead to the complete re-organisation of the microstructure (grains shape and size) and drastic modification of the texture (crystallographic orientations of grains) [2, 4–7]. These changes in the microstructure and texture may lead to modification of the rheological behavior, i.e., the response of the material to imposed stresses or strains.

Ice is a hexagonal crystalline material (Figure 1) that deforms on Earth at temperatures very close to its melting point. In the ductile regime, its deformation is characterized by a strong anisotropy, with dislocations gliding almost solely on the basal plane [8, 9]. Therefore, the orientation of an ice crystal relative to the direction of solicitation strongly controls its viscosity. Consequently, the distribution of crystallographic orientations in polycrystalline ice impacts the mechanical response in the form of texture-induced viscoplastic anisotropy [10].

Microstructural evidence for recrystallization is systematically observed along deep ice cores, where it is supposed to contribute variably to the measured texture, depending on the deformation conditions (temperature, impurities). To interpret texture evolution along ice cores, a model has been proposed by De la Chapelle et al. 1998 [11]. This model states that DRX may either slow down the strain-induced texture evolution (the so-called “continuous” dynamic recrystallization regime) or drastically modify the texture, a final texture controlled by the stress state (the “dis-continuous” dynamic recrystallization regime). Whether continuous or discontinuous DRX predominates depends on the temperature and strain conditions along the ice core [12, 13]. In deep sections of ice sheets, along ice streams, and in glacier areas where temperature and strain are high, DRX leads to highly anisotropic textures, particularly when simple shear is dominant (see e.g. [14–17]). These textures result in strong viscoplastic anisotropy, which can either enhance or slow down strain rates. For instance, the vertical cluster of *c*-axes usually observed along ice cores will accelerate horizontal shearing but slow down vertical compaction [18–23]. To model the deformation of an ice sheet, it is therefore essential to consider the effect of DRX on texture evolution.

Dynamic recrystallization in ice has been extensively studied in the laboratory. These studies allowed us to characterize the basic mechanisms at play, nucleation and grain-boundary migration, and the impact of strain heterogeneities on their kinetics (see e.g. [3, 24–31]). Kamb (1972) [4] was a pioneer in performing several creep experiments in simple shear and compression to characterize the effect of DRX on both texture evolution and macroscopic mechanical behavior. He proposed that the strain rate acceleration, which is characteristic of the tertiary creep stage, could result from the development of DRX textures favorable to the activation of basal slip. This study is likely at the origin of the widely established belief that DRX favors “well-oriented” grains relative to the imposed stress. Since then, work by Grennerat et al. 2012 [32] has revealed that the strain field during deformation in polycrystalline ice is highly heterogeneous, with no simple relationship between the intensity of the local strain and the Schmid factor, a parameter that quantifies the orientation relationship between a crystal and the imposed macroscopic stress. More recently, direct measures of the evolution of the strain field with DRX by Chauve et al. 2017 [33] have illustrated the heterogeneity of the mechanical behavior at the grain scale and its impact on DRX. They showed that DRX tends to homogenize the strain field in relation to local microstructure changes.
# 1. Introduction

Dynamic recrystallization (DRX) is a mechanism by which local stress incompatibilities and local strain heterogeneities are relaxed [1–3]. By *local*, we mean the scale of the grains that constitute crystalline materials such as ice, metals, and most rocks. Dynamic recrystallization occurs in plastically deformed materials at high temperatures. It is characterized by the nucleation of new grains and grain boundary migration and can lead to the complete re-organisation of the microstructure (grains shape and size) and drastic modification of the texture (crystallographic orientations of grains) [2, 4–7]. These changes in the microstructure and texture may lead to modification of the rheological behavior, i.e., the response of the material to imposed stresses or strains.

Ice is a hexagonal crystalline material (Figure 1) that deforms on Earth at temperatures very close to its melting point. In the ductile regime, its deformation is characterized by a strong anisotropy, with dislocations gliding almost solely on the basal plane [8, 9]. Therefore, the orientation of an ice crystal relative to the direction of solicitation strongly controls its viscosity. Consequently, the distribution of crystallographic orientations in polycrystalline ice impacts the mechanical response in the form of texture-induced viscoplastic anisotropy [10].

Microstructural evidence for recrystallization is systematically observed along deep ice cores, where it is supposed to contribute variably to the measured texture, depending on the deformation conditions (temperature, impurities). To interpret texture evolution along ice cores, a model has been proposed by De la Chapelle et al. 1998 [11]. This model states that DRX may either slow down the strain-induced texture evolution (the so-called “continuous” dynamic recrystallization regime) or drastically modify the texture, a final texture controlled by the stress state (the “dis-continuous” dynamic recrystallization regime). Whether continuous or discontinuous DRX predominates depends on the temperature and strain conditions along the ice core [12, 13]. In deep sections of ice sheets, along ice streams, and in glacier areas where temperature and strain are high, DRX leads to highly anisotropic textures, particularly when simple shear is dominant (see e.g. [14–17]). These textures result in strong viscoplastic anisotropy, which can either enhance or slow down strain rates. For instance, the vertical cluster of *c*-axes usually observed along ice cores will accelerate horizontal shearing but slow down vertical compaction [18–23]. To model the deformation of an ice sheet, it is therefore essential to consider the effect of DRX on texture evolution.

Dynamic recrystallization in ice has been extensively studied in the laboratory. These studies allowed us to characterize the basic mechanisms at play, nucleation and grain-boundary migration, and the impact of strain heterogeneities on their kinetics (see e.g. [3, 24–31]). Kamb (1972) [4] was a pioneer in performing several creep experiments in simple shear and compression to characterize the effect of DRX on both texture evolution and macroscopic mechanical behavior. He proposed that the strain rate acceleration, which is characteristic of the tertiary creep stage, could result from the development of DRX textures favorable to the activation of basal slip. This study is likely at the origin of the widely established belief that DRX favors “well-oriented” grains relative to the imposed stress. Since then, work by Grennerat et al. 2012 [32] has revealed that the strain field during deformation in polycrystalline ice is highly heterogeneous, with no simple relationship between the intensity of the local strain and the Schmid factor, a parameter that quantifies the orientation relationship between a crystal and the imposed macroscopic stress. More recently, direct measures of the evolution of the strain field with DRX by Chauve et al. 2017 [33] have illustrated the heterogeneity of the mechanical behavior at the grain scale and its impact on DRX. They showed that DRX tends to homogenize the strain field in relation to local microstructure changes.
Although this recent data show that DRX textures result from a complex interplay between the imposed stress or strain states, the local microstructure, and the resulting redistribution of stresses, it also corroborates the notion that the development of DRX textures is systematically associated with marked strain softening that signs the onset of tertiary creep [22, 25, 27]. Further texture evolution during tertiary creep results in a quasi-constant strain rate for bulk strains higher than 10%; this “stationary” strain rate is independent of the initial texture and controlled by the imposed stress [34–37]. This contrasts with the minimum strain rate observed at the onset of secondary creep, which depends on the initial texture and can be described by the classical Glen’s law (a Norton-Hoff type law), which nonlinearly relates the strain rate to the applied stress (through an exponent of 3 and an Arhenius type of dependence to the temperature) [9, 38]. This law does not hold in tertiary creep, but it has been adapted by modifying the stress exponent to a value slightly higher than 3 and adding a strain-rate enhancement factor (the ratio between the strain rate measured at tertiary creep and the minimum strain rate at secondary creep). This enhancement factor, between 4 in compression and 8 in simple shear [22, 37] accounts for the impact of the DRX texture on strain softening during tertiary creep.

Textures that develop in the ice during DRX are well characterized. Under uniaxial compression, the textures measured at tertiary creep are characterized by *c*-axes oriented within a girdle at an angle between 30° and 50° to the compression axis depending on the experimental temperature *T* and stress. These textures are stable, i.e., they evolve very slowly for strains higher than 10% [4, 24, 27, 29, 30, 37, 39]. During simple shear, tertiary creep textures seem to evolve with increasing strain from a two-maxima pattern, which is symmetric relative to the principal finite stretching direction, toward a strongly clustered texture, characterized by *c*-axes aligned in a direction perpendicular to the shear plane [25, 26, 31, 40, 41].

Texture development in ice deforming by dislocation creep has been successfully simulated using approaches of variable complexity [18, 19, 42, 43]. However, these approaches fail to reproduce the tertiary creep experimental textures, and the rare models simulating the effect of dynamic recrystallization are rather unsatisfactory in terms of either texture patterns or kinetics of texture evolution [44–48]. The discrepancies between models and observations result from the use of a simplified, mean-field approach that does not enable DRX mechanisms to be properly represented [44, 45], or from limitations in the deformation modeling frame despite the high complexity of microstructure evolution [46]. In [47, 48], the simulated textures are very close to the observed ones, but the kinetics of evolution are not. This discrepancy may result from the use of a Sachs approximation for modeling the mechanical response of the polycrystal, which results in a very soft mechanical response and incoherence in the grain rotation calculation [49].

In this work, we propose a new formulation to predict the impact of DRX on the texture evolution of polycrystalline materials such as ice. This formulation couples a continuous transverse isotropic (CTI) law, which models the anisotropic response of the ice crystal, to a crystal orientation evolution equation that accounts for DRX. The contribution of DRX to texture evolution is modeled by a physically based orientation attractor that maximizes the resolved shear stress on the easiest slip system in the crystal (basal slip for ice). This set of equations, named CTI-RX, is solved in a full-field model using the finite-element method (FEM).

In Section 2, we present the CTI-RX model: its equations for the CTI law (Section 2.1) and the proposed formulation for the evolution of the *c*-axis (Section 2.2). The numerical scheme designed to solve these equations is described in Section 3. Section 4 presents the calibration and validation for cases without dynamic recrystallization. Section 5 evaluates the capability of the model to reproduce the evolution of the texture and the mechanical response by comparing the predictions with experimental data for tertiary creep of ice polycrystals subjected to uniaxial compression, uniaxial tension, and simple shear.
Although this recent data show that DRX textures result from a complex interplay between
the imposed stress or strain states, the local microstructure, and the resulting redistribution of
stresses, it also corroborates the notion that the development of DRX textures is systematically
associated with marked strain softening that signs the onset of tertiary creep [22, 25, 27]. Further
texture evolution during tertiary creep results in a quasi-constant strain rate for bulk strains
higher than 10%; this "stationary" strain rate is independent of the initial texture and controlled
by the imposed stress [34–37]. This contrasts with the minimum strain rate observed at the onset
of secondary creep, which depends on the initial texture and can be described by the classical
Glen's law (a Norton-Hoff type law), which nonlinearly relates the strain rate to the applied stress
(through an exponent of 3 and an Arhenius type of dependence to the temperature) [9, 38]. This
law does not hold in tertiary creep, but it has been adapted by modifying the stress exponent
to a value slightly higher than 3 and adding a strain-rate enhancement factor (the ratio between
the strain rate measured at tertiary creep and the minimum strain rate at secondary creep). This
enhancement factor, between 4 in compression and 8 in simple shear [22, 37] accounts for the
impact of the DRX texture on strain softening during tertiary creep.

Textures that develop in the ice during DRX are well characterized. Under uniaxial compres-
sion, the textures measured at tertiary creep are characterized by **c**-axes oriented within a girdle
at an angle between 30° and 50° to the compression axis depending on the experimental tem-
perature *T* and stress. These textures are stable, i.e., they evolve very slowly for strains higher
than 10% [4, 24, 27, 29, 30, 37, 39]. During simple shear, tertiary creep textures seem to evolve
with increasing strain from a two-maxima pattern, which is symmetric relative to the principal
finite stretching direction, toward a strongly clustered texture, characterized by **c**-axes aligned in
a direction perpendicular to the shear plane [25, 26, 31, 40, 41].

Texture development in ice deforming by dislocation creep has been successfully simulated using approaches of variable complexity [18, 19, 42, 43]. However, these approaches fail to reproduce the tertiary creep experimental textures, and the rare models simulating the effect of dynamic recrystallization are rather unsatisfactory in terms of either texture patterns or kinetics of texture evolution [44–48]. The discrepancies between models and observations result from the use of a simplified, mean-field approach that does not enable DRX mechanisms to be properly represented [44, 45], or from limitations in the deformation modeling frame despite the high complexity of microstructure evolution [46]. In [47, 48], the simulated textures are very close to the observed ones, but the kinetics of evolution are not. This discrepancy may result from the use of a Sachs approximation for modeling the mechanical response of the polycrystal, which results in a very soft mechanical response and incoherence in the grain rotation calculation [49].

In this work, we propose a new formulation to predict the impact of DRX on the texture evo-
lution of polycrystalline materials such as ice. This formulation couples a continuous transverse
isotropic (CTI) law, which models the anisotropic response of the ice crystal, to a crystal orien-
tation evolution equation that accounts for DRX. The contribution of DRX to texture evolution is
modeled by a physically based orientation attractor that maximizes the resolved shear stress on
the easiest slip system in the crystal (basal slip for ice). This set of equations, named CTI-RX, is
solved in a full-field model using the finite-element method (FEM).

In Section 2, we present the CTI-RX model: its equations for the CTI law (Section 2.1) and
the proposed formulation for the evolution of the **c**-axis (Section 2.2). The numerical scheme
designed to solve these equations is described in Section 3. Section 4 presents the calibration
and validation for cases without dynamic recrystallization. Section 5 evaluates the capability of
the model to reproduce the evolution of the texture and the mechanical response by comparing
the predictions with experimental data for tertiary creep of ice polycrystals subjected to uniaxial
compression, uniaxial tension, and simple shear.
## 2. CTI-RX- a model for the ice single crystal behavior

### 2.1. Continuous transverse isotropic behavior

Following [50–54], we account for the viscoplastic anisotropy of the ice single crystal by modeling its mechanical behavior using the continuous transverse isotropic (CTI) law. In this framework, the viscoplastic response of the ice single crystal depends only on the orientation of the c-axis, being isotropic in the basal plane (see Figure 1). This formulation is based on theoretical work by Boehler (1987) [55], which has shown, through geometrical arguments, that the relation between the deviatoric stress S and the deviatoric strain rate D in a transverse isotropic symmetry may be described by a potential $\phi_S$:

$$D = \left( \frac{\partial \phi_S}{\partial S} \right)^D \qquad (1)$$

where $\phi_S$ is expressed on the $\mathcal{B}_{CTI}$ base, which is built from a combination of the following four invariants:

$$\mathcal{B}_{CTI} = [\mathrm{Tr}(\mathbf{S}^2), \mathrm{Tr}(\mathbf{S}^3), \mathrm{Tr}(\mathbf{MS}), \mathrm{Tr}(\mathbf{MS}^2)], \qquad (2)$$

With $\mathbf{M} = \mathbf{e}_3 \otimes \mathbf{e}_3$ where $\mathbf{e}_3$ is the direction of the plane of isotropy (Figure 1). The potential $\phi_S$ is described below in the linear ($\phi_S^{(1)}$) and nonlinear ($\phi_S^{(n)}$) cases.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** The transverse isotropic geometry is fully described by one direction, $e_3$, which defines the plane of isotropy. For an ice single crystal, this direction corresponds to that of the c-axis (i.e. the normal to the basal (0001) plane).</figcaption></figure>

#### 2.1.1. The linear CTI formulation

To obtain a linear CTI formulation, the potential, $\phi_S^{(1)}$, must be quadratic. It is therefore constructed from the quadratic invariants taken from the $\mathcal{B}_{CTI}$ base as:

$$\phi_S^{(1)} = \delta_1 \mathrm{Tr}(\mathbf{S}^2) + \delta_2 \mathrm{Tr}(\mathbf{MS})^2 + \delta_3 \mathrm{Tr}(\mathbf{MS}^2). \qquad (3)$$

Considering the following derivation rules,

$$\frac{\partial \mathrm{Tr}(\mathbf{MS})}{\partial \mathbf{S}} = \mathbf{M}, \quad \frac{\partial \mathrm{Tr}(\mathbf{S}^2)}{\partial \mathbf{S}} = 2\mathbf{S}, \quad \frac{\partial \mathrm{Tr}(\mathbf{MS}^2)}{\partial \mathbf{S}} = \mathbf{MS} + \mathbf{SM}, \qquad (4)$$

the constitutive relation is given by

$$D = \left( \frac{\partial \phi_S^{(1)}}{\partial S} \right)^D = 2\delta_1 S + 2\delta_2 M^D \mathrm{Tr}(\mathbf{MS}) + \delta_3 (\mathbf{MS} + \mathbf{SM})^D. \qquad (5)$$

The three independent parameters $\delta_i$ in Equation (5) can be fitted from the experimental data as detailed below.
## 2. CTI-RX- a model for the ice single crystal behavior

### 2.1. Continuous transverse isotropic behavior

Following [50–54], we account for the viscoplastic anisotropy of the ice single crystal by modeling its mechanical behavior using the continuous transverse isotropic (CTI) law. In this framework, the viscoplastic response of the ice single crystal depends only on the orientation of the c-axis, being isotropic in the basal plane (see Figure 1). This formulation is based on theoretical work by Boehler (1987) [55], which has shown, through geometrical arguments, that the relation between the deviatoric stress S and the deviatoric strain rate D in a transverse isotropic symmetry may be described by a potential $\phi_S$:

$$D = \left( \frac{\partial \phi_S}{\partial S} \right)^D \qquad (1)$$

where $\phi_S$ is expressed on the $\mathcal{B}_{CTI}$ base, which is built from a combination of the following four invariants:

$$\mathcal{B}_{CTI} = [\mathrm{Tr}(\mathbf{S}^2), \mathrm{Tr}(\mathbf{S}^3), \mathrm{Tr}(\mathbf{MS}), \mathrm{Tr}(\mathbf{MS}^2)], \qquad (2)$$

With $\mathbf{M} = \mathbf{e}_3 \otimes \mathbf{e}_3$ where $\mathbf{e}_3$ is the direction of the plane of isotropy (Figure 1). The potential $\phi_S$ is described below in the linear ($\phi_S^{(1)}$) and nonlinear ($\phi_S^{(n)}$) cases.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The transverse isotropic geometry is fully described by one direction, $\mathbf{e}_3$, which defines the plane of isotropy. For an ice single crystal, this direction corresponds to that of the c-axis (i.e. the normal to the basal (0001) plane).</figcaption></figure>

#### 2.1.1. The linear CTI formulation

To obtain a linear CTI formulation, the potential, $\phi_S^{(1)}$, must be quadratic. It is therefore constructed from the quadratic invariants taken from the $\mathcal{B}_{CTI}$ base as:

$$\phi_S^{(1)} = \delta_1 \mathrm{Tr}(\mathbf{S}^2) + \delta_2 \mathrm{Tr}(\mathbf{MS})^2 + \delta_3 \mathrm{Tr}(\mathbf{MS}^2). \qquad (3)$$

Considering the following derivation rules,

$$\frac{\partial \mathrm{Tr}(\mathbf{MS})}{\partial \mathbf{S}} = \mathbf{M}, \quad \frac{\partial \mathrm{Tr}(\mathbf{S}^2)}{\partial \mathbf{S}} = 2\mathbf{S}, \quad \frac{\partial \mathrm{Tr}(\mathbf{MS}^2)}{\partial \mathbf{S}} = \mathbf{MS} + \mathbf{SM}, \qquad (4)$$

the constitutive relation is given by

$$D = \left( \frac{\partial \phi_S^{(1)}}{\partial S} \right)^D = 2\delta_1 S + 2\delta_2 M^D \mathrm{Tr}(\mathbf{MS}) + \delta_3 (\mathbf{MS} + \mathbf{SM})^D. \qquad (5)$$

The three independent parameters $\delta_i$ in Equation (5) can be fitted from the experimental data as detailed below.
In the present implementation of the model, the anisotropy of the ice single crystal (or grain) behavior is described using the following parameters, where $^gX$ indicates a tensor expressed in the grain reference frame:

* $\psi_1$ is the fluidity of the ice crystal subjected to shear on the basal plane:

$$^gD_{13} = \frac{\psi_1}{2} g S_{13}, \quad ^gD_{23} = \frac{\psi_1}{2} g S_{23}, \qquad (6)$$

* $\beta$ is the ratio between the fluidity for shear parallel to the basal plane and that for shear normal to the basal plane:

$$^gD_{12} = \beta \frac{\psi_1}{2} g S_{12}, \qquad (7)$$

* $\gamma$ is the ratio of the viscosity for compression (or traction) Along the c-axis $\eta_c$ and for compression (or traction) in any direction $^g e_r$ within the basal plane:

$$^gS_{33} = 2\eta_c^g D_{33}, \quad \gamma^g S_{rr} = 2\eta_c^g D_{rr} \qquad (8)$$

For solicitations of similar amplitude in different directions, for instance along $e_3$ or perpendicular to $e_3$ (i.e. along $e_r$) $S = {}^gS_{33} = {}^gS_{rr}$, in the linear case, the strain rates are related by:

$$^gD_{rr} = \gamma^g D_{33} \qquad (9)$$

Using the above constraints on single crystal behavior, identification gives (see Appendix A for details):

$$\delta_1 = \frac{\psi_1 \beta}{4}, \quad \delta_2 = \frac{\psi_1}{2} \left( \beta \frac{\gamma + 2}{4\gamma - 1} - 1 \right), \quad \delta_3 = \frac{\psi_1}{2} (1 - \beta) \qquad (10)$$

The constitutive equation (5) can also be expressed in an inverse manner as follows:

$$S = 2\alpha_1 D + 2\alpha_2 M^D \mathrm{Tr}(MD) + \alpha_3 (MD + DM)^D \qquad (11)$$

and in terms of the potential

$$\phi_D^{(1)} = \alpha_1 \mathrm{Tr}(D^2) + \alpha_2 \mathrm{Tr}(MD)^2 + \alpha_3 \mathrm{Tr}(MD^2) \qquad (12)$$

with: $\alpha_1 = \frac{\eta_1}{\beta}$, $\alpha_2 = 2\eta_1(\frac{\gamma}{\beta} - 1)$, $\alpha_3 = 2\eta_1(1 - \frac{1}{\beta})$ and $\eta_1 = \frac{1}{\psi_1}$ (see [19] and appendix A for details).

### 2.1.2. The non-linear CTI formulation

A non-linear potential $\phi_D^{(n)}$ linking $S$ to $D^n$ is obtained by taking the linear potential $\tilde{\phi}_D^{(1)} = \frac{1}{2\eta_1}\phi_D^{(1)}$ to a higher order, $k$:

$$\phi_D^{(n)} = \frac{1}{k} (\tilde{\phi}_D^{(1)})^k . \qquad (13)$$

This hypothesis limits the number of independent parameters to 3 instead of 7 in the general case. A similar hypothesis was made in previous studies [19, 48, 56]. Using the chain rule for deriving $\phi_D^n$ gives:

$$S = \frac{\partial \phi_D^{(n)}}{\partial D} = (\tilde{\phi}_D^{(1)})^{k-1} \frac{\partial \tilde{\phi}_D^{(1)}}{\partial D} \qquad (14)$$

$k$ is further obtained by identification of the power exponent on Equation (14) in order to have $S \sim D^n$:

$$2(k-1)+2 = \frac{1}{n} \Rightarrow k = \frac{n+1}{2n} \qquad (15)$$

Therefore, a non-linear CTI law can be written as:

$$\phi_D^{(n)} = \frac{4n}{n+1} (\alpha_1 \mathrm{tr}(D^2) + \alpha_2 \mathrm{tr}(MD)^2 + \alpha_3 \mathrm{tr}(MD^2))^{\frac{n+1}{2n}}, \qquad (16)$$

Leading to the constitutive equation:

$$S = \eta_n^* (2\alpha_1 D + 2\alpha_2 M^D \mathrm{Tr}(MD) + \alpha_3 (MD + DM)^D), \qquad (17)$$
In the present implementation of the model, the anisotropy of the ice single crystal (or grain) behavior is described using the following parameters, where $^gX$ indicates a tensor expressed in the grain reference frame:

* $\psi_1$ is the fluidity of the ice crystal subjected to shear on the basal plane:

$$^gD_{13} = \frac{\psi_1}{2} g S_{13}, \quad ^gD_{23} = \frac{\psi_1}{2} g S_{23}, \qquad (6)$$

* $\beta$ is the ratio between the fluidity for shear parallel to the basal plane and that for shear normal to the basal plane:

$$^gD_{12} = \beta \frac{\psi_1}{2} g S_{12}, \qquad (7)$$

* $\gamma$ is the ratio of the viscosity for compression (or traction) Along the c-axis $\eta_c$ and for compression (or traction) in any direction $^g e_r$ within the basal plane:

$$^gS_{33} = 2\eta_c^g D_{33}, \quad \gamma^g S_{rr} = 2\eta_c^g D_{rr} \qquad (8)$$

For solicitations of similar amplitude in different directions, for instance along $e_3$ or perpendicular to $e_3$ (i.e. along $e_r$) $S = {}^gS_{33} = {}^gS_{rr}$, in the linear case, the strain rates are related by:

$$^gD_{rr} = \gamma^g D_{33} \qquad (9)$$

Using the above constraints on single crystal behavior, identification gives (see Appendix A for details):

$$\delta_1 = \frac{\psi_1 \beta}{4}, \quad \delta_2 = \frac{\psi_1}{2} \left( \beta \frac{\gamma + 2}{4\gamma - 1} - 1 \right), \quad \delta_3 = \frac{\psi_1}{2} (1 - \beta) \qquad (10)$$

The constitutive equation (5) can also be expressed in an inverse manner as follows:

$$S = 2\alpha_1 \mathbf{D} + 2\alpha_2 \mathbf{M}^D \mathrm{Tr}(\mathbf{MD}) + \alpha_3 (\mathbf{MD} + \mathbf{DM})^D \qquad (11)$$

and in terms of the potential

$$\phi_D^{(1)} = \alpha_1 \mathrm{Tr}(\mathbf{D}^2) + \alpha_2 \mathrm{Tr}(\mathbf{MD})^2 + \alpha_3 \mathrm{Tr}(\mathbf{MD}^2) \qquad (12)$$

with: $\alpha_1 = \frac{\eta_1}{\beta}$, $\alpha_2 = 2\eta_1(\frac{\gamma}{\beta} - 1)$, $\alpha_3 = 2\eta_1(1 - \frac{1}{\beta})$ and $\eta_1 = \frac{1}{\psi_1}$ (see [19] and appendix A for details).

### 2.1.2. The non-linear CTI formulation

A non-linear potential $\phi_D^{(n)}$ linking $S$ to $D^n$ is obtained by taking the linear potential $\tilde{\phi}_D^{(1)} = \frac{1}{2\eta_1} \phi_D^{(1)}$ to a higher order, $k$:

$$\phi_D^{(n)} = \frac{1}{k} (\tilde{\phi}_D^{(1)})^k . \qquad (13)$$

This hypothesis limits the number of independent parameters to 3 instead of 7 in the general case. A similar hypothesis was made in previous studies [19, 48, 56]. Using the chain rule for deriving $\phi_D^n$ gives:

$$S = \frac{\partial \phi_D^{(n)}}{\partial \mathbf{D}} = (\tilde{\phi}_D^{(1)})^{k-1} \frac{\partial \tilde{\phi}_D^{(1)}}{\partial \mathbf{D}} \qquad (14)$$

$k$ is further obtained by identification of the power exponent on Equation (14) in order to have $S \sim D^n$:

$$2(k-1)+2 = \frac{1}{n} \Rightarrow k = \frac{n+1}{2n} \qquad (15)$$

Therefore, a non-linear CTI law can be written as:

$$\phi_D^{(n)} = \frac{4n}{n+1} \left( \alpha_1 \mathrm{tr}(\mathbf{D}^2) + \alpha_2 \mathrm{tr}(\mathbf{MD})^2 + \alpha_3 \mathrm{tr}(\mathbf{MD}^2) \right)^{\frac{n+1}{2n}}, \qquad (16)$$

Leading to the constitutive equation:

$$S = \eta_n^* (2\alpha_1 \mathbf{D} + 2\alpha_2 \mathbf{M}^D \mathrm{Tr}(\mathbf{MD}) + \alpha_3 (\mathbf{MD} + \mathbf{DM})^D), \qquad (17)$$
with effective viscosity

$$ \eta_n^* = 2\eta_n(\alpha_1 \operatorname{tr}(\mathbf{D}^2) + \alpha_2 \operatorname{tr}(\mathbf{MD})^2 + \alpha_3 \operatorname{tr}(\mathbf{MD}^2))^{\frac{1-n}{2n}} \quad (18) $$

and $\alpha_i$ parameters:

$$ \alpha_1 = \frac{1}{2\beta}, \quad \alpha_2 = \left(\frac{\gamma}{\beta} - 1\right), \quad \alpha_3 = \left(1 - \frac{1}{\beta}\right). \quad (19) $$

It is worth noting that $\beta$ links $gD_{12}$ to $gS_{12}$ through

$$ gD_{12} = \frac{\psi_n}{2} \beta^{\frac{n+1}{2}} gS_{12}^n. \quad (20) $$

Therefore, $\beta$ depends on non-linearity, i.e., on the value of $n$. A similar enhancement is obtained for a linear behavior $n=1$ and a non-linear behavior with $n=3$ by setting $\beta_{n=3} = \sqrt{\beta_{n=1}}$.

## 2.2. Evolution of crystal rotation in response to deformation and recrystallization

Deformation resulting from dislocation glide in the different slip plane results in a rotation of the crystal whose rate is described as in [57, 58] by

$$ \frac{\partial \mathbf{c}}{\partial t} = \mathbf{W}\mathbf{c} - \lambda [\mathbf{D}\mathbf{c} - (\mathbf{c}^T \mathbf{D}\mathbf{c})\mathbf{c}], \quad (21) $$

where $\mathbf{W}(\mathbf{u}) = \frac{\nabla(\mathbf{u}) - \nabla^T(\mathbf{u})}{2}$ and $\mathbf{D}(\mathbf{u}) = \frac{\nabla(\mathbf{u}) + \nabla^T(\mathbf{u})}{2}$ are the spin- and strain-rate tensors. According to Equation (21), the rotation of the c-axis is due to the bulk spin, $\mathbf{Wc}$, and the viscoplastic spin expressed by $\lambda[\mathbf{Dc} - (\mathbf{c}^T \mathbf{Dc})\mathbf{c}]$. $\lambda$ is a parameter that depends on the plasticity model chosen for the ice crystal. In this frame, we assume that deformation is only due to the contribution of dislocations gliding on the basal plane, $\lambda$ is set to 1. A value of $\lambda < 1$ must match a plasticity model that include prismatic and pyramidal slip [57, 58]. This assumption is well funded for ice and has often been made in previous studies [45, 59–62].

To account for the impact of dynamic recrystallization on the evolution of the c-axis orientation, we add a new term to this equation based on the assumption that dynamic recrystallization induces a local re-orientation of the c-axis, either through nucleation of a grain with a new orientation or by migration of a grain boundary, and that this re-orientation leads to local softening. This continuous description of the impact of DRX on the evolution of the c-axis orientation is supported by experimental observations of DRX nucleation in ice. These observations show that nucleation by bulging and polygonization results in progressive misorientation of the recrystallized volumes relative to the parent grains, with misorientations occurring most often < 20° [33].

With the addition of the dynamic recrystallization term, the c-axis rotation is described as follows:

$$ \frac{\partial \mathbf{c}}{\partial t} = \mathbf{Wc} - \lambda [\mathbf{Dc} - (\mathbf{c}^T \mathbf{Dc})\mathbf{c}] + \frac{1}{\Gamma_{RX}} (\mathbf{c}_0 - \mathbf{c}), \quad (22) $$

where $\Gamma_{RX}$ is a parameter that controls the rate of rotation of the c-axis toward an attractor, $\mathbf{c}_0$, which represents the “ideal” orientation produced by recrystallization and whose formulation is described in Section 2.3. Therefore, the evolution of the c-axis orientation is controlled by the balance between (i) the deformation-induced rotation and (ii) the rotation toward a recrystallization-driven attractor, $\mathbf{c}_0$.

The meaning of $\Gamma_{RX}$ can be illustrated by considering the artificial case in which only dynamic recrystallization contributes to the c-axis rotation. In this case, Equation (22) becomes

$$ \frac{\partial \mathbf{c}}{\partial t} = \frac{1}{\Gamma_{RX}} (\mathbf{c}_0 - \mathbf{c}), \quad (23) $$
with effective viscosity

$$ \eta_n^* = 2\eta_n (\alpha_1 \operatorname{tr}(\mathbf{D}^2) + \alpha_2 \operatorname{tr}(\mathbf{MD})^2 + \alpha_3 \operatorname{tr}(\mathbf{MD}^2))^{\frac{1-n}{2n}} \quad (18) $$

and $\alpha_i$ parameters:

$$ \alpha_1 = \frac{1}{2\beta}, \quad \alpha_2 = \left(\frac{\gamma}{\beta} - 1\right), \quad \alpha_3 = \left(1 - \frac{1}{\beta}\right). \quad (19) $$

It is worth noting that $\beta$ links $gD_{12}$ to $gS_{12}$ through

$$ gD_{12} = \frac{\psi_n}{2} \beta^{\frac{n+1}{2}} gS_{12}^n. \quad (20) $$

Therefore, $\beta$ depends on non-linearity, i.e., on the value of $n$. A similar enhancement is obtained for a linear behavior $n=1$ and a non-linear behavior with $n=3$ by setting $\beta_{n=3} = \sqrt{\beta_{n=1}}$.

## 2.2. Evolution of crystal rotation in response to deformation and recrystallization

Deformation resulting from dislocation glide in the different slip plane results in a rotation of the crystal whose rate is described as in [57, 58] by

$$ \frac{\partial \mathbf{c}}{\partial t} = \mathbf{W}\mathbf{c} - \lambda [\mathbf{D}\mathbf{c} - (\mathbf{c}^T \mathbf{D}\mathbf{c})\mathbf{c}], \quad (21) $$

where $\mathbf{W}(\mathbf{u}) = \frac{\nabla(\mathbf{u}) - \nabla^T(\mathbf{u})}{2}$ and $\mathbf{D}(\mathbf{u}) = \frac{\nabla(\mathbf{u}) + \nabla^T(\mathbf{u})}{2}$ are the spin- and strain-rate tensors. According to Equation (21), the rotation of the c-axis is due to the bulk spin, $\mathbf{Wc}$, and the viscoplastic spin expressed by $\lambda[\mathbf{Dc} - (\mathbf{c}^T \mathbf{Dc})\mathbf{c}]$. $\lambda$ is a parameter that depends on the plasticity model chosen for the ice crystal. In this frame, we assume that deformation is only due to the contribution of dislocations gliding on the basal plane, $\lambda$ is set to 1. A value of $\lambda < 1$ must match a plasticity model that include prismatic and pyramidal slip [57, 58]. This assumption is well funded for ice and has often been made in previous studies [45, 59–62].

To account for the impact of dynamic recrystallization on the evolution of the c-axis orientation, we add a new term to this equation based on the assumption that dynamic recrystallization induces a local re-orientation of the c-axis, either through nucleation of a grain with a new orientation or by migration of a grain boundary, and that this re-orientation leads to local softening. This continuous description of the impact of DRX on the evolution of the c-axis orientation is supported by experimental observations of DRX nucleation in ice. These observations show that nucleation by bulging and polygonization results in progressive misorientation of the recrystallized volumes relative to the parent grains, with misorientations occurring most often < 20° [33].

With the addition of the dynamic recrystallization term, the c-axis rotation is described as follows:

$$ \frac{\partial \mathbf{c}}{\partial t} = \mathbf{Wc} - \lambda [\mathbf{Dc} - (\mathbf{c}^T \mathbf{Dc})\mathbf{c}] + \frac{1}{\Gamma_{RX}} (\mathbf{c}_0 - \mathbf{c}), \quad (22) $$

where $\Gamma_{RX}$ is a parameter that controls the rate of rotation of the c-axis toward an attractor, $\mathbf{c}_0$, which represents the “ideal” orientation produced by recrystallization and whose formulation is described in Section 2.3. Therefore, the evolution of the c-axis orientation is controlled by the balance between (i) the deformation-induced rotation and (ii) the rotation toward a recrystallization-driven attractor, $\mathbf{c}_0$.

The meaning of $\Gamma_{RX}$ can be illustrated by considering the artificial case in which only dynamic recrystallization contributes to the c-axis rotation. In this case, Equation (22) becomes

$$ \frac{\partial \mathbf{c}}{\partial t} = \frac{1}{\Gamma_{RX}} (\mathbf{c}_0 - \mathbf{c}), \quad (23) $$
The norm of $\frac{\partial c}{\partial t}$ can be interpreted as the rotation rate due only to dynamic recrystallization. This rate increases as $||c_0 - c||$. Considering a constant $c_0$, the solution for the c-axis rotation is:

$$c(t) = (c(0) - c_0) e^{-\frac{t}{\Gamma_{RX}}} + c(0). \qquad (24)$$

In this case, c rotates towards $c_0$ at a rate following an exponential decay with a characteristic time $\Gamma_{RX}$, such that 95% (99%) of the rotation is obtained after $4\Gamma_{RX}$ ($5\Gamma_{RX}$). In the CTI-RX model, $c_0$ is not constant but depends on the local diatoric stress tensor, S.

## 2.3. Formulation of the recrystallization attractor $c_0$

The formulation of $c_0$ stems from the assumption that a grain (or part of a grain) rotates toward, or/and nucleates in, an orientation that facilitates local viscoplastic deformation (dislocation glide). This assumption is supported by experimental observations that show that nucleation and grain boundary migration occur in areas of high strain incompatibility and produce orientations that tend to reduce this incompatibility [3, 33, 63]. Considering that the ice single crystal deforms mainly by dislocation glide on the basal plane, the recrystallization-induced rotation should maximize the resolved shear stress (RSS) on the basal plane. RSS for a dislocation with a Burgers vector $a_i$ (Figure 1) gliding on the basal plane c in response to a diatoric stress S is given by

$$RSS(a_i, c, S) = S : \mu, \qquad (25)$$

with $\mu = \frac{1}{2} \frac{c \otimes a_i + a_i \otimes c}{||c|| \cdot ||a_i||}$.

In the CTI formulation, the orientation for the maximum RSS is the orientation where the shear stress is maximum in the basal plane and the behavior is isotropic in this plane. Considering $s_1 > s_2 > s_3$ as the eigenvalues and $v_1$ ($i = 1-3$) as the corresponding eigenvectors of the diatoric stress tensor S, the maximum resolved shear stress is $||s_1 - s_3||$ and the associated orientations are [64]:

$$c_0 = \frac{1}{2} (v_1 \pm v_3) \qquad (26)$$

Degenerated cases exist if at least two eigenvalues are equal:

* $s_1 = s_2$, the possible attractors $c_0$ are the orientations at 45° from $v_3$.

* $s_2 = s_3$, the possible attractors $c_0$ are the orientations at 45° from $v_1$.

* $s_1 = s_2 = s_3 = 0$, the possible attractors $c_0$ are all orientations.

In all other cases, there are two possible orientations for the attractor $c_0$. To solve Equation (22) we define $c_0$ as which of the two solutions is closer to the actual c axis, in order to minimize the path between c and $c_0$. In practice, $c_0$ is chosen so that the scalar product between $c_0$ and c is maximal.

## 2.4. Interpretation of the orientation evolution equation

Equation (22) can be interpreted in terms of the relative impact of deformation versus recrystallization on the orientation of the c-axis. In the following section, we describe three end-member cases.

For uniaxial compression along the y-axis (Figure 2 (a)), the stable solution for the deformation term $(Wc - \lambda|Dc - (c^T Dc)c|)$ corresponds to c aligned with y (red dot, Figure 2 (a)) [19]. For recrystallization $(\frac{1}{\Gamma_{RX}}(c_0 - c))$, the stable solutions are the orientations at 45° to the y-axis (blue line, Figure 2 (a)). When both deformation and recrystallization are active, the stable orientations of the c-axis are located on a girdle between these two end-members (black line, Figure 2 (a)). The position of this girdle depends on the value of $\Gamma_{RX}$ (black arrows, Figure 2 (a)). The higher is $\Gamma_{RX}$ the closer is the girdle to the compression axis.
The norm of $\frac{\partial c}{\partial t}$ can be interpreted as the rotation rate due only to dynamic recrystallization. This rate increases as $\|c_0 - c\|$. Considering a constant $c_0$, the solution for the c-axis rotation is:

$$c(t) = (c(0) - c_0) e^{-\frac{t}{\Gamma_{RX}}} + c(0). \quad (24)$$

In this case, c rotates towards $c_0$ at a rate following an exponential decay with a characteristic time $\Gamma_{RX}$, such that 95% (99%) of the rotation is obtained after $4\Gamma_{RX}$ ($5\Gamma_{RX}$). In the CTI-RX model, $c_0$ is not constant but depends on the local diatoric stress tensor, S.

## 2.3. Formulation of the recrystallization attractor $c_0$

The formulation of $c_0$ stems from the assumption that a grain (or part of a grain) rotates toward, or/and nucleates in, an orientation that facilitates local viscoplastic deformation (dislocation glide). This assumption is supported by experimental observations that show that nucleation and grain boundary migration occur in areas of high strain incompatibility and produce orientations that tend to reduce this incompatibility [3, 33, 63]. Considering that the ice single crystal deforms mainly by dislocation glide on the basal plane, the recrystallization-induced rotation should maximize the resolved shear stress (RSS) on the basal plane. RSS for a dislocation with a Burgers vector $a_i$ (Figure 1) gliding on the basal plane c in response to a diatoric stress S is given by

$$RSS(a_i, c, S) = S : \mu, \quad (25)$$

with $\mu = \frac{1}{2} \frac{c \otimes a_i + a_i \otimes c}{\|c\| \|a_i\|}$.

In the CTI formulation, the orientation for the maximum RSS is the orientation where the shear stress is maximum in the basal plane and the behavior is isotropic in this plane. Considering $s_1 > s_2 > s_3$ as the eigenvalues and $v_1$ ($i = 1-3$) as the corresponding eigenvectors of the diatoric stress tensor S, the maximum resolved shear stress is $\|s_1 - s_3\|$ and the associated orientations are [64]:

$$c_0 = \frac{1}{2} (v_1 \pm v_3) \quad (26)$$

Degenerated cases exist if at least two eigenvalues are equal:

* $s_1 = s_2$, the possible attractors $c_0$ are the orientations at 45° from $v_3$.

* $s_2 = s_3$, the possible attractors $c_0$ are the orientations at 45° from $v_1$.

* $s_1 = s_2 = s_3 = 0$, the possible attractors $c_0$ are all orientations.

In all other cases, there are two possible orientations for the attractor $c_0$. To solve Equation (22) we define $c_0$ as which of the two solutions is closer to the actual c axis, in order to minimize the path between c and $c_0$. In practice, $c_0$ is chosen so that the scalar product between $c_0$ and c is maximal.

## 2.4. Interpretation of the orientation evolution equation

Equation (22) can be interpreted in terms of the relative impact of deformation versus recrystallization on the orientation of the c-axis. In the following section, we describe three end-member cases.

For uniaxial compression along the y-axis (Figure 2 (a)), the stable solution for the deformation term $(Wc - \lambda|Dc - (c^T Dc)c|)$ corresponds to c aligned with y (red dot, Figure 2 (a)) [19]. For recrystallization $(\frac{1}{\Gamma_{RX}}(c_0 - c))$, the stable solutions are the orientations at 45° to the y-axis (blue line, Figure 2 (a)). When both deformation and recrystallization are active, the stable orientations of the c-axis are located on a girdle between these two end-members (black line, Figure 2 (a)). The position of this girdle depends on the value of $\Gamma_{RX}$ (black arrows, Figure 2 (a)). The higher is $\Gamma_{RX}$ the closer is the girdle to the compression axis.
For uniaxial tension along y, the stable solutions for the deformation term are all the orientations located in the xOz plane (red line, Figure 2 (b)) [19]. When both deformation and recrystallization are active, the stable orientations (black circle, Figure 2 (b)) form a girdle between 45° and 90° to the y-axis depending on the value of $\Gamma_{RX}$.

For simple shear parallel to the xOz plane (with $\frac{\partial u_x}{\partial y} < 0$, Figure 2 (c)), there are four stable solutions of the recrystallization term that are $\pm y$ and $\pm x$ (Figure 2 (c), blue dots), and two for the deformation term, $\pm y$ (Figure 2 (c), red dots). The resulting stable orientations of the c-axis when both deformation and recrystallization are active are located at $\pm y$ and between x and $x+y$ (and the symmetrical equivalent) (black dashed line dots, Figure 2 (c)). This latter stable orientation, whose exact position depends on the value of $\Gamma_{RX}$, only exists for low $\Gamma_{RX}$ and vanishes for large $\Gamma_{RX}$ ($\Gamma_{RX} \to +\infty$).

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Pole figure schematic representation of the stable orientations for the evolution equation (eq. (22)) for various stress configurations S: (a) uniaxial compression along y, (b) uniaxial tension along y and (c) simple shear in the xOz plane, $\frac{\partial u_x}{\partial y} < 0$ (upper hemisphere).</figcaption></figure>

## 2.5. CTI-RX flow model description

The system of equations constituting the CTI-RX model is based on the approximation of an incompressible fluid behavior, to which we added equations accounting for the single-crystal anisotropic behavior and the effect of recrystallization on the evolution of the crystal orientation. It is composed of

$$ (1) \text{ The momentum equation under the Stokes assumption: } \qquad -\nabla \cdot \mathbf{S} + \nabla \cdot \mathbf{p} = 0 \qquad (27) $$

$$ (2) \text{ The non-linear CTI law: } \qquad \mathbf{S} = \eta_n^* (2\alpha_1 \mathbf{D} + 2\alpha_2 \mathbf{M}^D \mathrm{Tr}(\mathbf{M}\mathbf{D}) + \alpha_3 (\mathbf{M}\mathbf{D} + \mathbf{D}\mathbf{M})^D) \qquad (28) $$

$$ \eta_n^* = 2\eta_n (\alpha_1 \mathrm{tr}(\mathbf{D}^2) + \alpha_2 \mathrm{tr}(\mathbf{M}\mathbf{D})^2 + \alpha_3 \mathrm{tr}(\mathbf{M}\mathbf{D}^2))^{\frac{1-n}{2n}} $$
For uniaxial tension along y, the stable solutions for the deformation term are all the orientations located in the xOz plane (red line, Figure 2 (b)) [19]. When both deformation and recrystallization are active, the stable orientations (black circle, Figure 2 (b)) form a girdle between 45° and 90° to the y-axis depending on the value of $\Gamma_{RX}$.

For simple shear parallel to the xOz plane (with $\frac{\partial u_x}{\partial y} < 0$, Figure 2 (c)), there are four stable solutions of the recrystallization term that are $\pm y$ and $\pm x$ (Figure 2 (c), blue dots), and two for the deformation term, $\pm y$ (Figure 2 (c), red dots). The resulting stable orientations of the c-axis when both deformation and recrystallization are active are located at $\pm y$ and between x and $x+y$ (and the symmetrical equivalent) (black dashed line dots, Figure 2 (c)). This latter stable orientation, whose exact position depends on the value of $\Gamma_{RX}$, only exists for low $\Gamma_{RX}$ and vanishes for large $\Gamma_{RX}$ ($\Gamma_{RX} \to +\infty$).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Pole figure schematic representation of the stable orientations for the evolution equation (eq. (22)) for various stress configurations S: (a) uniaxial compression along y, (b) uniaxial tension along y and (c) simple shear in the xOz plane, $\frac{\partial u_x}{\partial y} < 0$ (upper hemisphere).</figcaption></figure>

## 2.5. CTI-RX flow model description

The system of equations constituting the CTI-RX model is based on the approximation of an incompressible fluid behavior, to which we added equations accounting for the single-crystal anisotropic behavior and the effect of recrystallization on the evolution of the crystal orientation. It is composed of

$$ (1) \text{ The momentum equation under the Stokes assumption: } \qquad -\nabla \cdot \mathbf{S} + \nabla p = 0 \qquad (27) $$

$$ (2) \text{ The non-linear CTI law: } \qquad \mathbf{S} = \eta_n^* (2\alpha_1 \mathbf{D} + 2\alpha_2 \mathbf{M}^D \mathrm{Tr}(\mathbf{M}\mathbf{D}) + \alpha_3 (\mathbf{M}\mathbf{D} + \mathbf{D}\mathbf{M})^D) \qquad (28) $$

$$ \eta_n^* = 2\eta_n (\alpha_1 \mathrm{tr}(\mathbf{D}^2) + \alpha_2 \mathrm{tr}(\mathbf{M}\mathbf{D})^2 + \alpha_3 \mathrm{tr}(\mathbf{M}\mathbf{D}^2))^{\frac{1-n}{2n}} $$
(3) The equation for the evolution of the c-axis:

$$ \frac{\partial \mathbf{c}}{\partial t} = \mathbf{W}\mathbf{c} - \lambda [\mathbf{D}\mathbf{c} - (\mathbf{c}^T\mathbf{D}\mathbf{c})\mathbf{c}] + \frac{1}{\Gamma_{RX}}(\mathbf{c}_0 - \mathbf{c}), \quad (29) $$

With $\mathbf{c}_0 = \frac{v_1 \pm v_3}{2}$ where $v_i$ are the eigenvectors of $\mathbf{S}$ and $\|\mathbf{c} \cdot \mathbf{c}_0\|$ is maximum.

### 3. Numerical scheme and implementation

The CTI-RX equations (Section 2.5) are solved on a cuboid of polycrystalline ice deformed by creep (stress boundary conditions) under either uniaxial compression, uniaxial tension, or simple shear (Figure 3). A full-field approach using FEM was chosen to describe the field of orientations $\mathbf{c}$ at the element size. This implementation results in the $R^3iCe$ model (Rheology, Recrystallization, Rheolef, in Continuous Transverse Isotropic material). This design allows direct comparisons of texture evolution and macroscopic mechanical behavior with data from laboratory creep experiments [29–31].

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Schematic representation of the rectangular cubic simulation domain, with surface and point boundaries labeled.</figcaption></figure>

#### 3.1. Numerical experiment setup

##### 3.1.1. Uniaxial creep tests

In uniaxial creep, a zero vertical velocity ($u_y(\text{bottom}) = 0$) condition is applied at the bottom surface to simulate a fixed plateau, as in laboratory creep tests. To avoid the translation of the sample, the position of one bottom corner is fully fixed ($u(\text{left\_bottom\_back}) = (0,0,0)$). To avoid rotation around the y-axis, a zero-velocity condition along the x-axis is also applied on a second summit ($u_x(\text{left\_bottom\_front}) = 0$). The sample is deformed by uniaxial compression (or tension) by applying a constant normal stress on the top surface ($F_y(\text{top}) = F$). By convention, $F < 0$ ($F > 0$) corresponds to compression (tension). This Neumann condition leads to a non-uniform displacement over the surface as weaker elements deform faster. This may hinder direct comparisons with laboratory experiments or other numerical models, in which homogeneous displacement conditions are imposed on the surface to which compression is applied. This issue has been addressed in Section 4.1.
(3) The equation for the evolution of the c-axis:

$$ \frac{\partial \mathbf{c}}{\partial t} = \mathbf{W}\mathbf{c} - \lambda [\mathbf{D}\mathbf{c} - (\mathbf{c}^T\mathbf{D}\mathbf{c})\mathbf{c}] + \frac{1}{\Gamma_{RX}}(\mathbf{c}_0 - \mathbf{c}), \quad (29) $$

With $\mathbf{c}_0 = \frac{v_1 \pm v_3}{2}$ where $v_i$ are the eigenvectors of $\mathbf{S}$ and $\|\mathbf{c} \cdot \mathbf{c}_0\|$ is maximum.

### 3. Numerical scheme and implementation

The CTI-RX equations (Section 2.5) are solved on a cuboid of polycrystalline ice deformed by creep (stress boundary conditions) under either uniaxial compression, uniaxial tension, or simple shear (Figure 3). A full-field approach using FEM was chosen to describe the field of orientations $\mathbf{c}$ at the element size. This implementation results in the $R^3iCe$ model (Rheology, Recrystallization, Rheolef, in Continuous Transverse Isotropic material). This design allows direct comparisons of texture evolution and macroscopic mechanical behavior with data from laboratory creep experiments [29–31].

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Schematic representation of the rectangular cubic simulation domain, with surface and point boundaries labeled.</figcaption></figure>

#### 3.1. Numerical experiment setup

##### 3.1.1. Uniaxial creep tests

In uniaxial creep, a zero vertical velocity ($u_y(\text{bottom}) = 0$) condition is applied at the bottom surface to simulate a fixed plateau, as in laboratory creep tests. To avoid the translation of the sample, the position of one bottom corner is fully fixed ($u(\text{left\_bottom\_back}) = (0,0,0)$). To avoid rotation around the y-axis, a zero-velocity condition along the x-axis is also applied on a second summit ($u_x(\text{left\_bottom\_front}) = 0$). The sample is deformed by uniaxial compression (or tension) by applying a constant normal stress on the top surface ($F_y(\text{top}) = F$). By convention, $F < 0$ ($F > 0$) corresponds to compression (tension). This Neumann condition leads to a non-uniform displacement over the surface as weaker elements deform faster. This may hinder direct comparisons with laboratory experiments or other numerical models, in which homogeneous displacement conditions are imposed on the surface to which compression is applied. This issue has been addressed in Section 4.1.
### 3.1.2. Simple shear creep tests

In the simple shear tests, forcing is applied as a constant tangential stress on the top boundary ($F_x$($top$) = $F$) of the sample. Other boundary conditions are:

* $\mathbf{u}$($bottom$) = 0,

* $u_y = u_z = 0$ for the lateral surfaces (left, right, front, back) to avoid solid rotation,

* $u_y$($top$) = $u_z$($top$) = 0.

### 3.2. Characteristic Numbers and time

When adimensionalized with respect to the material viscosity $\eta_n$ in Pa s$\frac{1}{n}$, the vertical extent of the cuboid $L$ in meters (Figure 3) and the applied macroscopic stress $\Sigma$ (in Pa), the system of equations becomes:

$$ -\tilde{\nabla}\tilde{\mathbf{S}} + \tilde{\nabla}\tilde{\mathbf{p}} = 0 \qquad (30) $$

$$ \begin{gathered} \tilde{\mathbf{S}} = \tilde{\eta}_n^* \left( 2\alpha_1 \tilde{\mathbf{D}} + 2\alpha_2 \mathbf{M}^D \mathrm{Tr}(\mathbf{M}\tilde{\mathbf{D}}) + \alpha_3 (\mathbf{M}\tilde{\mathbf{D}} + \tilde{\mathbf{D}}\mathbf{M})^D \right) \\ \tilde{\eta}_n^* = 2 \left( \alpha_1 \mathrm{Tr}(\tilde{\mathbf{D}}^2) + \alpha_2 \mathrm{Tr}(\mathbf{M}\tilde{\mathbf{D}})^2 + \alpha_3 \mathrm{Tr}(\mathbf{M}\tilde{\mathbf{D}}^2) \right)^{\frac{1-n}{2n}} \end{gathered} \qquad (31) $$

$$ \frac{\partial c}{\partial \tilde{t}} = \tilde{\mathbf{W}}(\tilde{\mathbf{u}}) \cdot \mathbf{c} - \lambda [\tilde{\mathbf{D}}(\tilde{\mathbf{u}}) \cdot \mathbf{c} - (\mathbf{c}^T \cdot \tilde{\mathbf{D}}(\tilde{\mathbf{u}}) \cdot \mathbf{c}) \cdot \mathbf{c}] + \mathcal{M}_o (\mathbf{c}_0 - \mathbf{c}) \qquad (32) $$

where the superscript (.) is used for all non-dimensional variables and operators.

The dimensionless number that arises when scaling the c-axis evolution equation is $\mathcal{M}_o$:

$$ \mathcal{M}_o = \left( \frac{\eta_n}{\Sigma} \right)^n \frac{1}{\Gamma_{RX}} \qquad (33) $$

$\mathcal{M}_o$ is the ratio between the deformation characteristic time $\tau = (\eta_n/\Sigma)^n$ and the recrystallization characteristic time $\Gamma_{RX}$, which controls the rate of rotation of c toward the attractor $\mathbf{c}_0$. This dimensionless number quantifies the relative weight of the contributions of recrystallization versus deformation to the rotation of the c-axis.

### 3.3. Numerical algorithm

Because the evolution of the c-axis orientation depends on the strain rate ($\tilde{\mathbf{D}}(\tilde{\mathbf{u}})$) and the effective viscosity $\eta^*$ depends on both the strain rate and the orientation of the c-axis, the system of coupled Equations (31) and (32) is non-linear. To solve this problem, we linearize the equation describing the evolution of the c-axis using a second-order backward differentiation formula (BDF2). The time discretization of the CTI-RX model results in the following system of equations:

$$ -\tilde{\nabla}\tilde{\mathbf{S}}^{t+1} + \tilde{\nabla}\tilde{\mathbf{p}}^{t+1} = 0 \qquad (34) $$

$$ \begin{gathered} \tilde{\mathbf{S}}^{t+1} = \tilde{\eta}_n^* \left( 2\alpha_1 \tilde{\mathbf{D}}^{t+1} + 2\alpha_2 \mathbf{M}^{t+1} \mathrm{Tr}(\mathbf{M}^{t+1} \tilde{\mathbf{D}}^{t+1}) + \alpha_3 (\mathbf{M}^{t+1} \tilde{\mathbf{D}}^{t+1} + \tilde{\mathbf{D}}^{t+1} \mathbf{M}^{t+1})^D \right) \\ \tilde{\eta}_n^* \left( \alpha_1 \mathrm{Tr}(\tilde{\mathbf{D}}^{t+1})^2 + \alpha_2 \mathrm{Tr}(\mathbf{M}^{t+1} \tilde{\mathbf{D}}^{t+1})^2 + \alpha_3 \mathrm{Tr}(\mathbf{M}^{t+1} (\tilde{\mathbf{D}}^{t+1})^2) \right)^{\frac{1-n}{2n}} \end{gathered} \qquad (35) $$

$$ \begin{gathered} \frac{3}{2\Delta\tilde{t}} \left( \mathbf{c}^{t+1} - \frac{4}{3}\mathbf{c}^t + \frac{1}{3}\mathbf{c}^{t-1} \right) \\ = \tilde{\mathbf{W}}(\tilde{\mathbf{u}}^{t+1}) \cdot \mathbf{c}^{t+1} - \lambda \left[ \tilde{\mathbf{D}}(\tilde{\mathbf{u}}^{t+1}) \cdot \mathbf{c}^{t+1} - (\mathbf{c}^{t+1})^T \cdot \tilde{\mathbf{D}}(\tilde{\mathbf{u}}^{t+1}) \cdot \mathbf{c}^{t+1} \right] \cdot \mathbf{c}^{t+1} + \mathcal{M}_o (\mathbf{c}_0^{t+1} - \mathbf{c}^{t+1}) \end{gathered} \qquad (36) $$

Finite elements and variational methods were used to solve the time-discretized problem on a Lagrangian grid within the C++ environment Rheolef [65]. The domain is spatially discretize using hexahedral finite elements (unless otherwise mentioned), and an initial c-axis orientation
### 3.1.2. Simple shear creep tests

In the simple shear tests, forcing is applied as a constant tangential stress on the top boundary ($F_x$($top$) = $F$) of the sample. Other boundary conditions are:

* $\mathbf{u}$($bottom$) = 0,

* $u_y = u_z = 0$ for the lateral surfaces (left, right, front, back) to avoid solid rotation,

* $u_y$($top$) = $u_z$($top$) = 0.

### 3.2. Characteristic Numbers and time

When adimensionalized with respect to the material viscosity $\eta_n$ in Pa s$\frac{1}{n}$, the vertical extent of the cuboid $L$ in meters (Figure 3) and the applied macroscopic stress $\Sigma$ (in Pa), the system of equations becomes:

$$ -\tilde{\nabla}\tilde{\mathbf{S}} + \tilde{\nabla}\tilde{\mathbf{p}} = 0 \qquad (30) $$

$$ \tilde{\mathbf{S}} = \tilde{\eta}_n^* \left( 2\alpha_1 \tilde{\mathbf{D}} + 2\alpha_2 \mathbf{M}^D \mathrm{Tr}(\mathbf{M}\tilde{\mathbf{D}}) + \alpha_3 (\mathbf{M}\tilde{\mathbf{D}} + \tilde{\mathbf{D}}\mathbf{M})^D \right) \\ \tilde{\eta}_n^* = 2 \left( \alpha_1 \mathrm{Tr}(\tilde{\mathbf{D}}^2) + \alpha_2 \mathrm{Tr}(\mathbf{M}\tilde{\mathbf{D}})^2 + \alpha_3 \mathrm{Tr}(\mathbf{M}\tilde{\mathbf{D}}^2) \right)^{\frac{1-n}{2n}} \qquad (31) $$

$$ \frac{\partial c}{\partial \tilde{t}} = \tilde{\mathbf{W}}(\tilde{\mathbf{u}}) \cdot \mathbf{c} - \lambda [\tilde{\mathbf{D}}(\tilde{\mathbf{u}}) \cdot \mathbf{c} - (\mathbf{c}^T \cdot \tilde{\mathbf{D}}(\tilde{\mathbf{u}}) \cdot \mathbf{c}) \cdot \mathbf{c}] + \mathcal{M}_o (\mathbf{c}_0 - \mathbf{c}) \qquad (32) $$

where the superscript (.) is used for all non-dimensional variables and operators.

The dimensionless number that arises when scaling the c-axis evolution equation is $\mathcal{M}_o$:

$$ \mathcal{M}_o = \left( \frac{\eta_n}{\Sigma} \right)^n \frac{1}{\Gamma_{RX}} \qquad (33) $$

$\mathcal{M}_o$ is the ratio between the deformation characteristic time $\tau = (\eta_n/\Sigma)^n$ and the recrystallization characteristic time $\Gamma_{RX}$, which controls the rate of rotation of c toward the attractor $\mathbf{c}_0$. This dimensionless number quantifies the relative weight of the contributions of recrystallization versus deformation to the rotation of the c-axis.

### 3.3. Numerical algorithm

Because the evolution of the c-axis orientation depends on the strain rate ($\tilde{\mathbf{D}}(\tilde{\mathbf{u}})$) and the effective viscosity $\eta^*$ depends on both the strain rate and the orientation of the c-axis, the system of coupled Equations (31) and (32) is non-linear. To solve this problem, we linearize the equation describing the evolution of the c-axis using a second-order backward differentiation formula (BDF2). The time discretization of the CTI-RX model results in the following system of equations:

$$ -\tilde{\nabla}\tilde{\mathbf{S}}^{t+1} + \tilde{\nabla}\tilde{\mathbf{p}}^{t+1} = 0 \qquad (34) $$

$$ \tilde{\mathbf{S}}^{t+1} = \tilde{\eta}_n^* \left( 2\alpha_1 \tilde{\mathbf{D}}^{t+1} + 2\alpha_2 \mathbf{M}^{t+1} \mathrm{Tr}(\mathbf{M}^{t+1} \tilde{\mathbf{D}}^{t+1}) + \alpha_3 (\mathbf{M}^{t+1} \tilde{\mathbf{D}}^{t+1} + \tilde{\mathbf{D}}^{t+1} \mathbf{M}^{t+1})^D \right) \qquad (35) $$

$$ \tilde{\eta}_n^* \quad t+1 = 2 \left( \alpha_1 \mathrm{Tr} \left( (\tilde{\mathbf{D}}^{t+1})^2 \right) + \alpha_2 \mathrm{Tr} \left( \mathbf{M}^{t+1} \tilde{\mathbf{D}}^{t+1} \right)^2 + \alpha_3 \mathrm{Tr} \left( \mathbf{M}^{t+1} (\tilde{\mathbf{D}}^{t+1})^2 \right) \right)^{\frac{1-n}{2n}} $$

$$ \frac{3}{2\Delta\tilde{t}} \left( \mathbf{c}^{t+1} - \frac{4}{3}\mathbf{c}^t + \frac{1}{3}\mathbf{c}^{t-1} \right) = \tilde{\mathbf{W}}(\tilde{\mathbf{u}}^{t+1}) \cdot \mathbf{c}^{t+1} - \lambda \left[ \tilde{\mathbf{D}}(\tilde{\mathbf{u}}^{t+1}) \cdot \mathbf{c}^{t+1} - (\mathbf{c}^{t+1})^T \cdot \tilde{\mathbf{D}}(\tilde{\mathbf{u}}^{t+1}) \cdot \mathbf{c}^{t+1} \right] \cdot \mathbf{c}^{t+1} + \mathcal{M}_o (\mathbf{c}_0^{t+1} - \mathbf{c}^{t+1}) \qquad (36) $$

Finite elements and variational methods were used to solve the time-discretized problem on a Lagrangian grid within the C++ environment Rheolef [65]. The domain is spatially discretize using hexahedral finite elements (unless otherwise mentioned), and an initial c-axis orientation
is prescribed for each element. The FE approximation for the tensorial mechanical parameter **M** and the scalar mechanical parameter **c** is computed with polynomials of the order 0 ($P_0$). The FE approximation for the tensorial strain rate $\tilde{\mathbf{D}}$, spin rate $\tilde{\mathbf{W}}$ and stress $\tilde{\mathbf{S}}$ fields is computed with polynomials of order 1 ($P_1$). The FE approximation for the displacement field is computed with polynomials of the order 2 ($P_2$).

The Rheolef *problem_mixed* solver is used to solve the displacement field $\tilde{\mathbf{u}}$ and the pressure field $\tilde{\mathbf{p}}$ under the incompressibility hypothesis for the linearized version of Equation (35), which is obtained by fixing the apparent viscosity $\eta^*$. A first fix loop is used to solve the non-linearity associated with $\eta^*$. This loop is nested in a second fix loop that solves the temporal evolution of the **c**-axis orientation using Equation (36). In practice, it is recommended for numerical efficiency to flatten the two nested fix loops.

?? 1 presents the numerical algorithm of two nested fix points implemented in R³iCe.

This algorithm is a semi-implicit solver for the coupled Equations (31) and (32). The simulated orientation field **c** is thus consistent with the displacement field $\tilde{\mathbf{u}}$ and therefore with the diatoriic strain rate $\tilde{\mathbf{D}}$ and the deviatoric stress $\tilde{\mathbf{S}}$. This implementation of CTI-RX model is distributed as the R³iCe code.

**Algorithm 1** Numerical algorithm, based on two nested fix loops: one for solving the viscosity, $\eta_n^*$, nested in a second one solving the temporal evolution of the orientation of **c**-axis.

1: **Inputs:**
   $\mathbf{c}^{\text{init}}$
2: **Initialize:**
   $\mathbf{c}^{-1} \leftarrow \mathbf{c}^{\text{init}}$
   $\mathbf{c}^{0} \leftarrow \mathbf{c}^{\text{init}}$
   $\mathbf{u}^{0} \leftarrow 0$
3: **for** $t$ from 0 to $N$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;$k_c \leftarrow 0$
5: &nbsp;&nbsp;&nbsp;&nbsp;$\text{res}_c \leftarrow +\infty$
6: &nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{u}^{t+1}, k_c=0, k_{nl}=0 \leftarrow \mathbf{u}^t$
7: &nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{c}^{t+1}, k_c=0 \leftarrow \mathbf{c}^t$
8: &nbsp;&nbsp;&nbsp;&nbsp;**while** $\text{res}_c > \text{tol}_c$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k_{nl} \leftarrow 0$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{res}_{nl} \leftarrow +\infty$
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{u}^{t+1}, k_c, k_{nl}=0 \leftarrow \mathbf{u}^{t+1}, k_c$
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**while** $\text{res}_{nl} > \text{tol}_{nl}$ **do**
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$H_n^* \leftarrow \eta_n^*(\mathbf{u}^{t+1}, k_c, k_{nl}, \mathbf{c}^{t+1}, k_c)$
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;find $\mathbf{u}^{t+1}, k_c, k_{nl}+1 \leftarrow \mathbf{u}^{t+1}, k_c, k_{nl}+1 (H_n^*, \mathbf{c}^{t+1}, k_c)$
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{res}_{nl} \leftarrow \|\mathbf{S}(\eta_n^*(\mathbf{u}^{t+1}, k_c, k_{nl}+1, \mathbf{c}^{t+1}, k_c), \mathbf{u}^{t+1}, k_c, k_{nl}+1, \mathbf{c}^{t+1}, k_c) - \mathbf{S}(H_n^*, \mathbf{u}^{t+1}, k_c, k_{nl}+1, \mathbf{c}^{t+1}, k_c)\|$}
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k_{nl} \leftarrow k_{nl} + 1$
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end while**
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$H_n^* \leftarrow \eta_n^*(\mathbf{u}^{t+1}, k_c, k_{nl}, \mathbf{c}^{t+1}, k_c)$
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;compute $\mathbf{S} \leftarrow \mathbf{S}(H_n^*, \mathbf{u}^{t+1}, k_c, k_{nl})$
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;find $\mathbf{c}_0 \leftarrow \mathbf{c}_0(\mathbf{S})$
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;compute $\mathbf{c}^{t+1}, k_c+1 \leftarrow \mathbf{c}^{t+1}, k_c+1 (\mathbf{c}_0, \mathbf{u}^{t+1}, k_c, k_{nl})$
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{res}_c \leftarrow \|\mathbf{S}(\eta_n^*(\mathbf{u}^{t+1}, k_c, k_{nl}, \mathbf{c}^{t+1}, k_c+1), \mathbf{u}^{t+1}, k_c, k_{nl}+1, \mathbf{c}^{t+1}, k_c+1) - \mathbf{S}(H_n^*, \mathbf{u}^{t+1}, k_c, k_{nl}, \mathbf{c}^{t+1}, k_c)\|$}
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k_c \leftarrow k_c + 1$
24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end while**
25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{u}^{t+1} \leftarrow \mathbf{u}^{t+1}, k_c, k_{nl}$
26: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{c}^{t+1} \leftarrow \mathbf{c}^{t+1}, k_c$
is prescribed for each element. The FE approximation for the tensorial mechanical parameter **M** and the scalar mechanical parameter **c** is computed with polynomials of the order 0 ($P_0$). The FE approximation for the tensorial strain rate $\tilde{\mathbf{D}}$, spin rate $\tilde{\mathbf{W}}$ and stress $\tilde{\mathbf{S}}$ fields is computed with polynomials of order 1 ($P_1$). The FE approximation for the displacement field is computed with polynomials of the order 2 ($P_2$).

The Rheolef *problem_mixed* solver is used to solve the displacement field $\tilde{\mathbf{u}}$ and the pressure field $\tilde{\mathbf{p}}$ under the incompressibility hypothesis for the linearized version of Equation (35), which is obtained by fixing the apparent viscosity $\eta^*$. A first fix loop is used to solve the non-linearity associated with $\eta^*$. This loop is nested in a second fix loop that solves the temporal evolution of the **c**-axis orientation using Equation (36). In practice, it is recommended for numerical efficiency to flatten the two nested fix loops.

?? 1 presents the numerical algorithm of two nested fix points implemented in R³iCe.

This algorithm is a semi-implicit solver for the coupled Equations (31) and (32). The simulated orientation field **c** is thus consistent with the displacement field $\tilde{\mathbf{u}}$ and therefore with the diatori cal strain rate $\tilde{\mathbf{D}}$ and the deviatoric stress $\tilde{\mathbf{S}}$. This implementation of CTI-RX model is distributed as the R³iCe code.

**Algorithm 1** Numerical algorithm, based on two nested fix loops: one for solving the viscosity, $\eta_n^*$, nested in a second one solving the temporal evolution of the orientation of **c**-axis.

1: **Inputs:**
   $\mathbf{c}^{\text{init}}$
2: **Initialize:**
   $\mathbf{c}^{-1} \leftarrow \mathbf{c}^{\text{init}}$
   $\mathbf{c}^{0} \leftarrow \mathbf{c}^{\text{init}}$
   $\mathbf{u}^{0} \leftarrow 0$
3: **for** $t$ from 0 to $N$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;$k_c \leftarrow 0$
5: &nbsp;&nbsp;&nbsp;&nbsp;$\text{res}_c \leftarrow +\infty$
6: &nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{u}^{t+1}, k_c=0, k_{nl}=0 \leftarrow \mathbf{u}^t$
7: &nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{c}^{t+1}, k_c=0 \leftarrow \mathbf{c}^t$
8: &nbsp;&nbsp;&nbsp;&nbsp;**while** $\text{res}_c > \text{tol}_c$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k_{nl} \leftarrow 0$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{res}_{nl} \leftarrow +\infty$
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{u}^{t+1}, k_c, k_{nl}=0 \leftarrow \mathbf{u}^{t+1}, k_c$
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**while** $\text{res}_{nl} > \text{tol}_{nl}$ **do**
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$H_n^* \leftarrow \eta_n^*(\mathbf{u}^{t+1}, k_c, k_{nl}, \mathbf{c}^{t+1}, k_c)$
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;find $\mathbf{u}^{t+1}, k_c, k_{nl}+1 \leftarrow \mathbf{u}^{t+1}, k_c, k_{nl}+1$ ($H_n^*, \mathbf{c}^{t+1}, k_c$) ▷ *problem_mixed* solver on Equation (31)
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{res}_{nl} \leftarrow \|\mathbf{S}(\eta_n^*(\mathbf{u}^{t+1}, k_c, k_{nl}+1, \mathbf{c}^{t+1}, k_c), \mathbf{u}^{t+1}, k_c, k_{nl}+1, \mathbf{c}^{t+1}, k_c) - \mathbf{S}(H_n^*, \mathbf{u}^{t+1}, k_c, k_{nl}+1, \mathbf{c}^{t+1}, k_c)\|$}
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k_{nl} \leftarrow k_{nl} + 1$
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end while**
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$H_n^* \leftarrow \eta_n^*(\mathbf{u}^{t+1}, k_c, k_{nl}, \mathbf{c}^{t+1}, k_c)$
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;compute $\mathbf{S} \leftarrow \mathbf{S}(H_n^*, \mathbf{u}^{t+1}, k_c, k_{nl})$ ▷ Equation (31)
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;find $\mathbf{c}_0 \leftarrow \mathbf{c}_0(\mathbf{S})$ ▷ Equation (26)
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;compute $\mathbf{c}^{t+1}, k_c+1 \leftarrow \mathbf{c}^{t+1}, k_c+1$ ($\mathbf{c}_0, \mathbf{u}^{t+1}, k_c, k_{nl}$) ▷ Equation (36)
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{res}_c \leftarrow \|\mathbf{S}(\eta_n^*(\mathbf{u}^{t+1}, k_c, k_{nl}, \mathbf{c}^{t+1}, k_c+1), \mathbf{u}^{t+1}, k_c, k_{nl}+1, \mathbf{c}^{t+1}, k_c+1) - \mathbf{S}(H_n^*, \mathbf{u}^{t+1}, k_c, k_{nl}, \mathbf{c}^{t+1}, k_c)\|$}
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k_c \leftarrow k_c + 1$
24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end while**
25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{u}^{t+1} \leftarrow \mathbf{u}^{t+1}, k_c, k_{nl}$
26: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{c}^{t+1} \leftarrow \mathbf{c}^{t+1}, k_c$
<table><caption>**Table 2.** Parameters of the CTI law for the ice single crystal.</caption><thead><tr><th>n</th><th>γ</th><th>β</th><th>η<sub>n</sub></th></tr></thead><tbody><tr><td>3</td><td>1</td><td>10<sup>-2</sup></td><td>7.5 MPa.s<sup>1/3</sup></td></tr></tbody></table>

### 3.4. CTI parameter values for the ice single crystal

As stated in Section 2, the values for the CTI parameters, $γ$, $η_n$, $n$ and $β$ (Equations (17) and (19)) are defined on the basis of experimental data.

$γ$ corresponds to the ratio of viscosities measured in compression and tension for an ice single crystal that is well oriented for basal slip. Laboratory experiments have shown that this ratio is very close to one. Therefore, in the present numerical simulations, we set $γ = 1$, as in previous studies [50, 51].

The following Duval et al. (1983) [9] and previous modeling work based on the CTI law [20, 48, 51], we set $n = 3$. While experiments in ice single crystals that are well oriented for basal slip are best fitted by $n = 2$, experiments on crystals that are poorly oriented for basal slip result in $n = 3$ [9]. Moreover, a value of 3 must be imposed in the CTI formulation to recover $n = 3$ for the polycrystal response.

$β$ is linked to the ratio between the basal and non-basal shearing viscosities. The strong viscoplastic anisotropy of ice is mainly carried by this parameter. In creep experiments at 1 MPa [9], the ratio between basal and non-basal shear viscosities is on the order of $10^4$. This ratio is related to $β$ through Equation (20), $β^{n+1/2} = 10^4$. For $n = 3$, $β = 10^{-2}$ is a value that is consistent with those considered in previous studies using the non-linear CTI formulation [20, 48, 51].

The viscosity $η_n$ can be adjusted by comparison with either the experimental response of a single crystal or that of isotropic polycrystals. Solving the CTI law for the uniaxial compression creep of a single crystal oriented at $45°$ from the compression axis and of a polycrystal with random orientations picked from a uniform texture results in a macroscopic strain rate of $\tilde{D}_{sc}^{sim} = 1.7 \times 10^{-2}(\frac{\eta_n}{\Sigma})^n$ for the single crystal and $\tilde{D}_{px}^{sim} = 5.3 \times 10^{-5}(\frac{\eta_n}{\Sigma})^n$ for the polycrystal. The viscosity $η_n$ is obtained by comparing the simulated values with the experimental data at 1 MPa [9] ($D_{sc}^{exp} = (2 \times 10^{-5}, 2 \times 10^{-4}) s^{-1}$ for the single crystal, $D_{px}^{exp} = (9 \times 10^{-8}, 3 \times 10^{-7}) s^{-1}$ for the polycrystal). $η_3^{sc} = (4.4, 9.4)$ MPa s$^{1/3}$ for the single crystal and $η_3^{px} = (5.6, 8.4)$ MPa s$^{1/3}$ for the polycrystal. The value of 7.5 MPa s$^{1/3}$ is taken for $η_3$ in order to reproduce the experiment duration as observed during uniaxial compression tests that are taken here as a reference (see Section 5.2). The parameters of the CTI law for the ice single crystal are summarized in Table 2.

### 4. Validation of the CTI-based implementation for predicting viscoplastic response

In this section, we verify that the CTI law implemented in R³iCe correctly reproduces the mechanical response of an ice polycrystal deforming by dislocation creep only (no DRX). This is a prerequisite for improving the CTI formulation with a representation of DRX. The mechanical response of polycrystalline ice is known to be characterized by strong strain and stress heterogeneities, which are very likely precursors of dynamic recrystallization [32]. We first compare the R³iCe predictions to those of the CraFT-EVP model [32], which has been shown to adequately reproduce the heterogeneous strain fields observed in experiments [32, 66]. We then investigate the reliability of a simplified mesh configuration for simulating the evolution of the texture and its impact on the macroscopic mechanical response in a numerically efficient manner.
<table><caption>**Table 2.** Parameters of the CTI law for the ice single crystal.</caption><thead><tr><td>n</td><td>γ</td><td>β</td><td>η<sub>n</sub></td></tr></thead><tbody><tr><td>3</td><td>1</td><td>10<sup>-2</sup></td><td>7.5 MPa.s<sup>1/3</sup></td></tr></tbody></table>

### 3.4. CTI parameter values for the ice single crystal

As stated in Section 2, the values for the CTI parameters, $\gamma$, $\eta_n$, $n$ and $\beta$ (Equations (17) and (19)) are defined on the basis of experimental data.

$\gamma$ corresponds to the ratio of viscosities measured in compression and tension for an ice single crystal that is well oriented for basal slip. Laboratory experiments have shown that this ratio is very close to one. Therefore, in the present numerical simulations, we set $\gamma = 1$, as in previous studies [50, 51].

The following Duval et al. (1983) [9] and previous modeling work based on the CTI law [20, 48, 51], we set $n = 3$. While experiments in ice single crystals that are well oriented for basal slip are best fitted by $n = 2$, experiments on crystals that are poorly oriented for basal slip result in $n = 3$ [9]. Moreover, a value of 3 must be imposed in the CTI formulation to recover $n = 3$ for the polycrystal response.

$\beta$ is linked to the ratio between the basal and non-basal shearing viscosities. The strong viscoplastic anisotropy of ice is mainly carried by this parameter. In creep experiments at 1 MPa [9], the ratio between basal and non-basal shear viscosities is on the order of $10^4$. This ratio is related to $\beta$ through Equation (20), $\beta^{n+1/2} = 10^4$. For $n = 3$, $\beta = 10^{-2}$ is a value that is consistent with those considered in previous studies using the non-linear CTI formulation [20, 48, 51].

The viscosity $\eta_n$ can be adjusted by comparison with either the experimental response of a single crystal or that of isotropic polycrystals. Solving the CTI law for the uniaxial compression creep of a single crystal oriented at $45^\circ$ from the compression axis and of a polycrystal with random orientations picked from a uniform texture results in a macroscopic strain rate of $\tilde{D}_{sc}^{sim} = 1.7 \times 10^{-2}(\frac{\eta_n}{\Sigma})^n$ for the single crystal and $\tilde{D}_{px}^{sim} = 5.3 \times 10^{-5}(\frac{\eta_n}{\Sigma})^n$ for the polycrystal. The viscosity $\eta_n$ is obtained by comparing the simulated values with the experimental data at 1 MPa [9] ($D_{sc}^{exp} = (2 \times 10^{-5}, 2 \times 10^{-4}) \text{ s}^{-1}$ for the single crystal, $D_{px}^{exp} = (9 \times 10^{-8}, 3 \times 10^{-7}) \text{ s}^{-1}$ for the polycrystal). $\eta_3^{sc} = (4.4, 9.4) \text{ MPa s}^{1/3}$ for the single crystal and $\eta_3^{px} = (5.6, 8.4) \text{ MPa s}^{1/3}$ for the polycrystal. The value of 7.5 MPa s$^{1/3}$ is taken for $\eta_3$ in order to reproduce the experiment duration as observed during uniaxial compression tests that are taken here as a reference (see Section 5.2). The parameters of the CTI law for the ice single crystal are summarized in Table 2.

### 4. Validation of the CTI-based implementation for predicting viscoplastic response

In this section, we verify that the CTI law implemented in R³iCe correctly reproduces the mechanical response of an ice polycrystal deforming by dislocation creep only (no DRX). This is a prerequisite for improving the CTI formulation with a representation of DRX. The mechanical response of polycrystalline ice is known to be characterized by strong strain and stress heterogeneities, which are very likely precursors of dynamic recrystallization [32]. We first compare the R³iCe predictions to those of the CraFT-EVP model [32], which has been shown to adequately reproduce the heterogeneous strain fields observed in experiments [32, 66]. We then investigate the reliability of a simplified mesh configuration for simulating the evolution of the texture and its impact on the macroscopic mechanical response in a numerically efficient manner.
4.1. *Stress and strain field predictions compared with CraFT-EVP*

CraFT-EVP has been used to provide full-field predictions of the mechanical response of ice during primary and secondary creep [32]. It is based on an elasto-viscoplastic formulation proposed by Suquet et al. 2012 [67] and solved using CraFT software. The predictions of CraFT-EVP for the mechanical response and strain fields were validated by direct comparisons with mechanical tests on ice polycrystals that included strain field measurements by Digital Image Correlation [32, 67]. The macroscopic strain fields predicted by CraFT-EVP agreed well with the experimental observations at 1% strain (secondary creep, before any texture evolution owing to DRX). With no texture evolution, CraFT-EVP and R³iCe both simulate a stationary secondary creep.

Reproducing accurate stress fields in R³iCe is key because the recrystallization attractor c₀ depends on the deviatoric stress S at the element scale (Section 2.3). Note that the Elasto-ViscoPlastic (EVP) law implemented in CraFT-EVP [67] is based on a full description of crystal plasticity, in which all possible slip systems contribute to the deformation as a function of their critical resolved shear stress and crystal orientation, whereas the CTI law in R³iCe is a simplified parameterization of viscoplastic anisotropy, which has been calibrated directly from experimental data (Section 3.4). The objective of the comparison performed below is to verify that the stress field predicted by R³iCe is, at minimum, statistically representative of a realistic stress field. Since stress fields have never been measured in polycrystalline ice, the stress field predicted by CraFT-EVP is here considered as a reference. R³iCe parameters were, of course, not fitted to the CraFT-EVP viscoplastic response, so the two predictions can be considered independent.

Figure 4 (a-b) shows the polycrystal microstructure generated by a grain growth algorithm using Neper [68] used for both CraFT-EVP and R³iCe simulations. In this configuration, grains are defined as a group of connected elements that share the same initial orientation. CraFT-EVP requires a regular cubic mesh, while the R³iCe polycrystalline microstructure used here was meshed with a tetrahedral geometry. The resolution of CraFT-EVP simulations was reduced to match the R³iCe one. Figure 5 shows the component εyy of the strain field and the component σyy of the stress field on a xy section at z = 0, for a uniaxial compression creep test and after 1% of macroscopic strain. The simulated strain fields do not perfectly compare between the CraFT-EVP and R³iCe predictions although the overall patterns do not show a strong mismatch. The locations of areas in compression (blue) and extension (red) are similarly prescribed (Figure 5 (a)). Significant differences between the responses of the models CraFT-EVP and R³iCe are restricted to the first layers of grains close to the top surface. These differences result from the different boundary conditions applied in the two models. In R³iCe, the Neumann boundary condition on the top surface results in a constant force applied on each element that can lead to heterogeneous displacement on this surface. As a result, some “weaker” grain deform faster, as represented by the dark blue areas on Figure 5 (a). In CraFT-EVP, periodic boundary conditions are imposed because equations are solved with a fast Fourier transform-based method. The simulated stress fields (Figure 5 (b)) compare satisfactorily in terms of amplitude and global pattern, although the precise stress localization shows more scattering. Some of this mismatch may also be explained by the differences in boundary conditions between the two models.

Similarly, the distributions of predicted Von Mises equivalent stress and strain by the two models show some discrepancies (Figure 6). These discrepancies are stronger between the two configurations in which grains are described by multiple elements (Figure 6) and where the impact of the boundary conditions mentioned above is stronger. Distributions for all components of both tensors (εij, σij) are given in Appendix B.
## 4.1. *Stress and strain field predictions compared with CraFT-EVP*

CraFT-EVP has been used to provide full-field predictions of the mechanical response of ice during primary and secondary creep [32]. It is based on an elasto-viscoplastic formulation proposed by Suquet et al. 2012 [67] and solved using CraFT software. The predictions of CraFT-EVP for the mechanical response and strain fields were validated by direct comparisons with mechanical tests on ice polycrystals that included strain field measurements by Digital Image Correlation [32, 67]. The macroscopic strain fields predicted by CraFT-EVP agreed well with the experimental observations at 1% strain (secondary creep, before any texture evolution owing to DRX). With no texture evolution, CraFT-EVP and R³iCe both simulate a stationary secondary creep.

Reproducing accurate stress fields in R³iCe is key because the recrystallization attractor $c_0$ depends on the deviatoric stress S at the element scale (Section 2.3). Note that the Elasto-ViscoPlastic (EVP) law implemented in CraFT-EVP [67] is based on a full description of crystal plasticity, in which all possible slip systems contribute to the deformation as a function of their critical resolved shear stress and crystal orientation, whereas the CTI law in R³iCe is a simplified parameterization of viscoplastic anisotropy, which has been calibrated directly from experimental data (Section 3.4). The objective of the comparison performed below is to verify that the stress field predicted by R³iCe is, at minimum, statistically representative of a realistic stress field. Since stress fields have never been measured in polycrystalline ice, the stress field predicted by CraFT-EVP is here considered as a reference. R³iCe parameters were, of course, not fitted to the CraFT-EVP viscoplastic response, so the two predictions can be considered independent.

Figure 4 (a-b) shows the polycrystal microstructure generated by a grain growth algorithm using Neper [68] used for both CraFT-EVP and R³iCe simulations. In this configuration, grains are defined as a group of connected elements that share the same initial orientation. CraFT-EVP requires a regular cubic mesh, while the R³iCe polycrystalline microstructure used here was meshed with a tetrahedral geometry. The resolution of CraFT-EVP simulations was reduced to match the R³iCe one. Figure 5 shows the component $\epsilon_{yy}$ of the strain field and the component $\sigma_{yy}$ of the stress field on a xy section at $z = 0$, for a uniaxial compression creep test and after 1% of macroscopic strain. The simulated strain fields do not perfectly compare between the CraFT-EVP and R³iCe predictions although the overall patterns do not show a strong mismatch. The locations of areas in compression (blue) and extension (red) are similarly prescribed (Figure 5 (a)). Significant differences between the responses of the models CraFT-EVP and R³iCe are restricted to the first layers of grains close to the top surface. These differences result from the different boundary conditions applied in the two models. In R³iCe, the Neumann boundary condition on the top surface results in a constant force applied on each element that can lead to heterogeneous displacement on this surface. As a result, some “weaker” grain deform faster, as represented by the dark blue areas on Figure 5 (a). In CraFT-EVP, periodic boundary conditions are imposed because equations are solved with a fast Fourier transform-based method. The simulated stress fields (Figure 5 (b)) compare satisfactorily in terms of amplitude and global pattern, although the precise stress localization shows more scattering. Some of this mismatch may also be explained by the differences in boundary conditions between the two models.

Similarly, the distributions of predicted Von Mises equivalent stress and strain by the two models show some discrepancies (Figure 6). These discrepancies are stronger between the two configurations in which grains are described by multiple elements (Figure 6) and where the impact of the boundary conditions mentioned above is stronger. Distributions for all components of both tensors ($\epsilon_{ij}, \sigma_{ij}$) are given in Appendix B.
In R³iCe, the stress tensor at the element scale is used to predict the recrystallization attrac-
tor **c**₀. The latter is based on the principal direction of the stress tensor **v**₁ and **v**₃ (Section 2.3).
Figure 7 shows the distribution of **v**₁ and **v**₃ orientations predicted using both CraFT-EVP (Fig-
ure 7 (a)), and R³iCe with a configuration with several elements per grain (Figure 7 (b)) and one
orientation per element (Figure 7 (c)). All three configurations predict similar **v**₁ and **v**₃ orienta-
tion distributions and therefore similar **c**₀ orientation distributions are expected.

In summary, if CraFT-EVP simulations are taken as a reference, the strain and stress fields at
secondary creep, prior to DRX activation, are qualitatively and statistically well predicted by the
CTI constitutive law implemented in R³iCe. In particular, the principal stress tensor directions
used in R³iCe to simulate the dynamic recrystallization impact on orientation evolution are
similarly predicted by the two models.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Microstructure generated via a grain growth using Neper meshed for (a) CraFT-EVP microstructure with 265302 grid points (b) R³iCe with 262688 tetrahedral elements. (c) Simplified microstructure in which orientations are defined at the element scale with a hexahedral mesh of 2740 elements.</figcaption></figure>

## 4.2. Stress and strain field prediction in the simplified mesh configuration

Resolving CTI-RX with sub-grain scale meshes is numerically very costly. In the following, we
test a microstructure representation composed of one orientation (or grain) per mesh element
(Figure 4 (c)) in order to keep a statistically reasonable number of orientations on a mesh size
(~ 2700 elements).

The distributions of the Von Mises equivalent stress and strain in these one-orientation-per-
mesh-element simulations are compared with those of the simulations with meshes at the sub-
grain scale using both the CraFT-EVP and R³iCe models (Figure 6). The distributions obtained in
this simplified microstructure representation compare well with those predicted using the CraFT-
EVP model, which is here considered as a reference. The fact that the CTI-RX fields predictions
In R³iCe, the stress tensor at the element scale is used to predict the recrystallization attrac-
tor **c**₀. The latter is based on the principal direction of the stress tensor **v**₁ and **v**₃ (Section 2.3).
Figure 7 shows the distribution of **v**₁ and **v**₃ orientations predicted using both CraFT-EVP (Fig-
ure 7 (a)), and R³iCe with a configuration with several elements per grain (Figure 7 (b)) and one
orientation per element (Figure 7 (c)). All three configurations predict similar **v**₁ and **v**₃ orienta-
tion distributions and therefore similar **c**₀ orientation distributions are expected.

In summary, if CraFT-EVP simulations are taken as a reference, the strain and stress fields at
secondary creep, prior to DRX activation, are qualitatively and statistically well predicted by the
CTI constitutive law implemented in R³iCe. In particular, the principal stress tensor directions
used in R³iCe to simulate the dynamic recrystallization impact on orientation evolution are
similarly predicted by the two models.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Microstructure generated via a grain growth using Neper meshed for (a) CraFT-EVP microstructure with 265302 grid points (b) R³iCe with 262688 tetrahedral elements. (c) Simplified microstructure in which orientations are defined at the element scale with a hexahedral mesh of 2740 elements.</figcaption></figure>

## 4.2. Stress and strain field prediction in the simplified mesh configuration

Resolving CTI-RX with sub-grain scale meshes is numerically very costly. In the following, we
test a microstructure representation composed of one orientation (or grain) per mesh element
(Figure 4 (c)) in order to keep a statistically reasonable number of orientations on a mesh size
(~ 2700 elements).

The distributions of the Von Mises equivalent stress and strain in these one-orientation-per-
mesh-element simulations are compared with those of the simulations with meshes at the sub-
grain scale using both the CraFT-EVP and R³iCe models (Figure 6). The distributions obtained in
this simplified microstructure representation compare well with those predicted using the CraFT-
EVP model, which is here considered as a reference. The fact that the CTI-RX fields predictions
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of predicted (a) strain ($\varepsilon_{yy}$) (b) stress ($\sigma_{yy} = S_{yy} + p$) fields for the grain microstructure using CraFT-EVP and R³iCe.</figcaption></figure>

obtained with a sub-grain scale mesh depart slightly from the Craft-EVP predictions is attributed, as mentioned before, to the discrepancy in the applied boundary conditions (periodic versus Neumann applied on the top surface). Please note that this discrepancy is reduced by a larger number of grains (orientations) in contact with the top surface in the simplified microstructure.

Figure 7 compares the distribution of the principal stress tensor directions $v_1$ and $v_3$ that are used in the calculation of the dynamic recrystallization attractor in R³iCe. There is a superb correspondence between the distributions predicted by CraFT-EVP and by R³iCe with the simplified microstructure. In the following, all simulations were performed using one independent orientation per element to have a large enough number of orientations while maintaining reasonable computation times. We assume that this choice has a limited impact on the simulated stress fields and therefore on the evaluation of the attractor $c_0$ at the sample scale.

**5. Texture evolution: comparison with experimental data**

This section tests the ability of the R³iCe model to predict the evolution of the c-axes orientation with and without dynamic recrystallization (Equation (32)). As mentioned in the introduction, deformation experiments with no or very limited DRX are not possible because of the fast strain rates and high temperature required. Deformation-only textures may only be sampled in deep ice cores from the central parts of Greenland or Antarctica [11]. These textures were well reproduced by uniaxial compression simulations performed using the mean-field ViscoPlastic Self Consistent model [13, 43]. R³iCe prediction for texture evolution without recrystallization are therefore compared with those of VPSC simulations extracted from [13]. In contrast, there is a large experimental dataset in which the texture evolved in response to both deformation and DRX. We use this dataset to constrain the predictions of R³iCe for the effect of DRX on texture evolution and mechanical behavior under uniaxial compression, uniaxial tension, and simple shear.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of predicted (a) strain ($\varepsilon_{yy}$) (b) stress ($\sigma_{yy} = S_{yy} + p$) fields for the grain microstructure using CraFT-EVP and R³iCe.</figcaption></figure>

obtained with a sub-grain scale mesh depart slightly from the Craft-EVP predictions is attributed, as mentioned before, to the discrepancy in the applied boundary conditions (periodic versus Neumann applied on the top surface). Please note that this discrepancy is reduced by a larger number of grains (orientations) in contact with the top surface in the simplified microstructure.

Figure 7 compares the distribution of the principal stress tensor directions $v_1$ and $v_3$ that are used in the calculation of the dynamic recrystallization attractor in R³iCe. There is a superb correspondence between the distributions predicted by CraFT-EVP and by R³iCe with the simplified microstructure. In the following, all simulations were performed using one independent orientation per element to have a large enough number of orientations while maintaining reasonable computation times. We assume that this choice has a limited impact on the simulated stress fields and therefore on the evaluation of the attractor $c_0$ at the sample scale.

**5. Texture evolution: comparison with experimental data**

This section tests the ability of the R³iCe model to predict the evolution of the c-axes orientation with and without dynamic recrystallization (Equation (32)). As mentioned in the introduction, deformation experiments with no or very limited DRX are not possible because of the fast strain rates and high temperature required. Deformation-only textures may only be sampled in deep ice cores from the central parts of Greenland or Antarctica [11]. These textures were well reproduced by uniaxial compression simulations performed using the mean-field ViscoPlastic Self Consistent model [13, 43]. R³iCe prediction for texture evolution without recrystallization are therefore compared with those of VPSC simulations extracted from [13]. In contrast, there is a large experimental dataset in which the texture evolved in response to both deformation and DRX. We use this dataset to constrain the predictions of R³iCe for the effect of DRX on texture evolution and mechanical behavior under uniaxial compression, uniaxial tension, and simple shear.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Kernel Density Estimation of the Von Mises equivalent strain (a) and stress (b) fields simulated with CraFT-EVP (EVP, black line) and R³iCe (CTI, cyan line). The dashed lines show the results obtained for simulations where grains are defined by multiple mesh elements (*grains* microstructure) (Figure 4 (a-b)). The full line shows a R³iCe simulation where a single orientation is attributed per element (Figure 4 (c)).</figcaption></figure>

The laboratory creep experiments selected to validate the model are:

* Uniaxial compression experiments from Montagnat et al. 2015 [29] ($\Sigma \sim 0.7$ MPa, $T = -5^\circ C$)

* Uniaxial tension experiments from Jacka and Maccagnan 1984 [27] ($\Sigma \sim 0.4$ MPa, $T = -3^\circ C$)

* Torsion experiments from Journaux et al. 2019 [31] ($\Sigma_{max} \sim 0.5$ MPa, $T = -7^\circ C$),

where $\Sigma$ stands for the macroscopic applied stress and $\Sigma_{max}$ corresponds to the maximum stress resulting from the applied torque.

The model predictions are also compared with strain-controlled uniaxial compression tests under confining pressure performed by Qi et al. 2017 [30]. These conditions enabled us to reach higher flow stresses without failure ($\Sigma \sim (1.3, 4.3)$ MPa, $T = -10^\circ C$).

For each configuration, simulations using a configuration with one independent orientation per element are performed, starting with initial orientations randomly selected in a uniform texture. The parameters for all simulations are given in Table 3.

For all configurations, the simulated texture evolution is represented by (1) the evolution of the eigenvalues $a_i$ of the second-order orientation tensor $\mathbf{a}^{(2)}$ ($\mathbf{a}^{(2)} = \frac{1}{N} \sum_{i=0}^{N} \mathbf{c}_i \otimes \mathbf{c}_i$, with N being the number of orientations), (2) pole figures representing the c-axes distribution at various strains, and (3) the evolution with strain of the texture in a section extracted from the pole figure. These
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Kernel Density Estimation of the Von Mises equivalent strain (a) and stress (b) fields simulated with CraFT-EVP (EVP, black line) and R³iCe (CTI, cyan line). The dashed lines show the results obtained for simulations where grains are defined by multiple mesh elements (*grains* microstructure) (Figure 4 (a-b)). The full line shows a R³iCe simulation where a single orientation is attributed per element (Figure 4 (c)).</figcaption></figure>

The laboratory creep experiments selected to validate the model are:

* Uniaxial compression experiments from Montagnat et al. 2015 [29] ($\Sigma \sim 0.7$ MPa, $T = -5^\circ C$)

* Uniaxial tension experiments from Jacka and Maccagnan 1984 [27] ($\Sigma \sim 0.4$ MPa, $T = -3^\circ C$)

* Torsion experiments from Journaux et al. 2019 [31] ($\Sigma_{max} \sim 0.5$ MPa, $T = -7^\circ C$),

where $\Sigma$ stands for the macroscopic applied stress and $\Sigma_{max}$ corresponds to the maximum stress resulting from the applied torque.

The model predictions are also compared with strain-controlled uniaxial compression tests under confining pressure performed by Qi et al. 2017 [30]. These conditions enabled us to reach higher flow stresses without failure ($\Sigma \sim (1.3, 4.3)$ MPa, $T = -10^\circ C$).

For each configuration, simulations using a configuration with one independent orientation per element are performed, starting with initial orientations randomly selected in a uniform texture. The parameters for all simulations are given in Table 3.

For all configurations, the simulated texture evolution is represented by (1) the evolution of the eigenvalues $a_i$ of the second-order orientation tensor $\mathbf{a}^{(2)}$ ($\mathbf{a}^{(2)} = \frac{1}{N} \sum_{i=0}^{N} \mathbf{c}_i \otimes \mathbf{c}_i$, with N being the number of orientations), (2) pole figures representing the c-axes distribution at various strains, and (3) the evolution with strain of the texture in a section extracted from the pole figure. These
outputs are compared with the pole figures obtained during the experiments at three selected strain values.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Distributions of the principal stress tensor directions **v**₁ and **v**₃ for uniaxial compression on (a) CraFT-EVP simulation (b) R³iCe simulation with grains (c) R³iCe simulation with one orientation per mesh element.</figcaption></figure>

<table><caption>**Table 3.** Values of the numerical and physical parameters used in the model R³iCe.</caption><thead><tr><th rowspan="2">References</th><th colspan="2">Model inputs</th><th colspan="3">Physical values</th></tr><tr><th>Stress</th><th>μ</th><th>Σ (MPa)</th><th>Γ<sub>RX</sub> (days)</th><th>T (°C)</th></tr></thead><tbody><tr><td>without recrystallization</td><td>Compression</td><td>0</td><td>1</td><td>+∞</td><td></td></tr><tr><td>Montagnat et al. 2015 [29]</td><td>Compression</td><td>1.03 × 10<sup>-3</sup></td><td>0.7</td><td>13.8</td><td>-5</td></tr><tr><td>Jacka et al. 1984 [27]</td><td>Tension</td><td>1.86 × 10<sup>-2</sup></td><td>0.4</td><td>4.1</td><td>-3</td></tr><tr><td>Journaux et al. 2019 [31]</td><td>Simple Shear</td><td>1.06 × 10<sup>-3</sup></td><td>0.5</td><td>36.8</td><td>-7</td></tr><tr><td>Qi et al. 2017 [30]</td><td>Compression</td><td>(1.34, 20.8) × 10<sup>-5</sup></td><td>(1.3, 4.3)</td><td>4.6</td><td>-10</td></tr></tbody></table>
outputs are compared with the pole figures obtained during the experiments at three selected strain values.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Distributions of the principal stress tensor directions **v**₁ and **v**₃ for uniaxial compression on (a) CraFT-EVP simulation (b) R³iCe simulation with grains (c) R³iCe simulation with one orientation per mesh element.</figcaption></figure>

<table><caption>**Table 3.** Values of the numerical and physical parameters used in the model R³iCe.</caption><thead><tr><th rowspan="2">References</th><th colspan="2">Model inputs</th><th colspan="3">Physical values</th></tr><tr><th>Stress</th><th>μ</th><th>Σ (MPa)</th><th>Γ<sub>RX</sub> (days)</th><th>T (°C)</th></tr></thead><tbody><tr><td>without recrystallization</td><td>Compression</td><td>0</td><td>1</td><td>+∞</td><td></td></tr><tr><td>Montagnat et al. 2015 [29]</td><td>Compression</td><td>1.03 × 10<sup>-3</sup></td><td>0.7</td><td>13.8</td><td>-5</td></tr><tr><td>Jacka et al. 1984 [27]</td><td>Tension</td><td>1.86 × 10<sup>-2</sup></td><td>0.4</td><td>4.1</td><td>-3</td></tr><tr><td>Journaux et al. 2019 [31]</td><td>Simple Shear</td><td>1.06 × 10<sup>-3</sup></td><td>0.5</td><td>36.8</td><td>-7</td></tr><tr><td>Qi et al. 2017 [30]</td><td>Compression</td><td>(1.34, 20.8) × 10<sup>-5</sup></td><td>(1.3, 4.3)</td><td>4.6</td><td>-10</td></tr></tbody></table>
The second-order orientation tensor $\mathbf{a}^{(2)}$ is a statistical representation of texture anisotropy [69]. The c-axes distribution lies within an ellipsoid whose axes are the eigenvalues of $\mathbf{a}^{(2)}$ with $1 \ge a_1 \ge a_2 \ge a_3 > 0$.

The simulated strain-rate evolution provides an evaluation of the evolution of the enhancement factor, which is the ratio between the maximum strain rate measured at the onset of DRX and the minimum strain rate measured at secondary creep in the experiments, which is the starting point of the simulations. In the simulation, this strain-rate enhancement results solely from the evolution of the texture, whereas in the experiments, changes in the microstructure might also contribute to it.

## 5.1. *Texture evolution in the absence of recrystallization*

To test the capability of R³iCe to predict texture evolution, we simulated uniaxial compression under a 1 MPa deviatoric stress using one orientation per mesh element with initial orientations randomly selected from a uniform texture.

Figure 8 (a-c) presents the evolution with strain (i.e. time) of the eigenvalues ($a_i, i = 1, 3$) of the second-order orientation tensor $\mathbf{a}^{(2)}$, the strain rate, and the texture in the Z plane (perpendicular to the compression direction). Figure 8 (d) shows textures predicted at 22%, 40% and 60% total strain. The strain-rate curve discontinuities (Figure 8 (b)) are technical artifacts imposed by the operating rules of the GRICAD clusters. The shared access restricts the simulation duration to 48 hours, which means that the code must be relaunched every 48 hours. To reduce the calculation time, the convergence criterion is set tougher for the first step and relaxed for the following steps. The steps observed at each restart reveal the convergence error during the simulation, which appears to have a limited impact.

A strong cluster (single maximum of c-axes) texture develops as strain increases. This texture evolution induces geometrical hardening of the polycrystal, which is expressed by a decrease in the strain rate by a factor ~ 1.5 between 20% and 60% total strain, when the texture is highly clustered ($a_1 \approx 0.71$).

Similar cluster textures are observed along deep cold ice cores, for instance, along the Talos Dome ice core [13] where a clustered texture with $a_1 \approx 0.7$ was measured at about 650 m depth, where an accumulated strain of 60% is estimated. The simulated texture evolution with strain is also in good agreement (although slightly weaker) with the one predicted by Montagnat et al. [13] using the VPSC model when departing from an isotropic texture (see red dots in Figure 8 (a)).

## 5.2. *Uniaxial compression*

The evolution with strain (i.e. time) of the eigenvalues ($a_i, i = 1, 3$) of the second-order orientation tensor $\mathbf{a}^{(2)}$, the strain rate, and the texture in the Z plane (perpendicular to the compression direction) predicted for uniaxial compression with DRX under an imposed compressive stress of 0.7 MPa is shown in Figure 9. This is compared with experimental measurements from Montagnat et al. 2015 [29]. These latter were obtained from initially isotropic granular ice compressed at -5°C under a constant load with initial stresses of 0.7 and 0.8 MPa, and up to various strains (from 2 to 17.8% bulk shortening).

Comparison of the c-axis pole figures shows good qualitative agreement between the simulated and experimental textures (Figure 9 (d-e)). Both are girdle textures that become more concentrated and show an evolution of the cone angle to the compression axis from 45° to 35° up to ~ 20% macroscopic strain and remain nearly constant thereafter (Figure 9 (d)). There is also good agreement between the predicted evolution of the texture eigenvalues and those measured
The second-order orientation tensor $\mathbf{a}^{(2)}$ is a statistical representation of texture anisotropy [69]. The c-axes distribution lies within an ellipsoid whose axes are the eigenvalues of $\mathbf{a}^{(2)}$ with $1 \ge a_1 \ge a_2 \ge a_3 > 0$.

The simulated strain-rate evolution provides an evaluation of the evolution of the enhancement factor, which is the ratio between the maximum strain rate measured at the onset of DRX and the minimum strain rate measured at secondary creep in the experiments, which is the starting point of the simulations. In the simulation, this strain-rate enhancement results solely from the evolution of the texture, whereas in the experiments, changes in the microstructure might also contribute to it.

## 5.1. *Texture evolution in the absence of recrystallization*

To test the capability of R³iCe to predict texture evolution, we simulated uniaxial compression under a 1 MPa deviatoric stress using one orientation per mesh element with initial orientations randomly selected from a uniform texture.

Figure 8 (a-c) presents the evolution with strain (i.e. time) of the eigenvalues ($a_i, i = 1, 3$) of the second-order orientation tensor $\mathbf{a}^{(2)}$, the strain rate, and the texture in the Z plane (perpendicular to the compression direction). Figure 8 (d) shows textures predicted at 22%, 40% and 60% total strain. The strain-rate curve discontinuities (Figure 8 (b)) are technical artifacts imposed by the operating rules of the GRICAD clusters. The shared access restricts the simulation duration to 48 hours, which means that the code must be relaunched every 48 hours. To reduce the calculation time, the convergence criterion is set tougher for the first step and relaxed for the following steps. The steps observed at each restart reveal the convergence error during the simulation, which appears to have a limited impact.

A strong cluster (single maximum of c-axes) texture develops as strain increases. This texture evolution induces geometrical hardening of the polycrystal, which is expressed by a decrease in the strain rate by a factor ~ 1.5 between 20% and 60% total strain, when the texture is highly clustered ($a_1 \approx 0.71$).

Similar cluster textures are observed along deep cold ice cores, for instance, along the Talos Dome ice core [13] where a clustered texture with $a_1 \approx 0.7$ was measured at about 650 m depth, where an accumulated strain of 60% is estimated. The simulated texture evolution with strain is also in good agreement (although slightly weaker) with the one predicted by Montagnat et al. [13] using the VPSC model when departing from an isotropic texture (see red dots in Figure 8 (a)).

## 5.2. *Uniaxial compression*

The evolution with strain (i.e. time) of the eigenvalues ($a_i, i = 1, 3$) of the second-order orientation tensor $\mathbf{a}^{(2)}$, the strain rate, and the texture in the Z plane (perpendicular to the compression direction) predicted for uniaxial compression with DRX under an imposed compressive stress of 0.7 MPa is shown in Figure 9. This is compared with experimental measurements from Montagnat et al. 2015 [29]. These latter were obtained from initially isotropic granular ice compressed at -5°C under a constant load with initial stresses of 0.7 and 0.8 MPa, and up to various strains (from 2 to 17.8% bulk shortening).

Comparison of the c-axis pole figures shows good qualitative agreement between the simulated and experimental textures (Figure 9 (d-e)). Both are girdle textures that become more concentrated and show an evolution of the cone angle to the compression axis from 45° to 35° up to ~ 20% macroscopic strain and remain nearly constant thereafter (Figure 9 (d)). There is also good agreement between the predicted evolution of the texture eigenvalues and those measured
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Simulation of uniaxial compression along *Y* with no recrystallization Γ<sub>RX</sub> = +∞, *M*<sub>o</sub> = 0. a) Evolution with strain of the eigenvalues of the second-order orientation tensor (*a*<sub>1</sub> > *a*<sub>2</sub> > *a*<sub>3</sub>). Red dots are predictions from the VPSC model taken from [13]; b) evolution with strain of the macroscopic strain rate; c) evolution of the intensity of the texture with the colatitude *θ* over a section defined on the *Z* plane (dashed blue lines in d)). Data at 22%, 40% and 60% macrostrain are shown. The white dashed lines define the angle at ±45° from the compression axis *Y*; d) textures presented as pole figures for 22% (*i*), 40% (*ii*), and 60% (*iii*) macroscopic strain. All figures presenting texture data are color-coded using the same colorbar, which is displayed on the right (scales in multiples of a uniform distribution).</figcaption></figure>

at 7%, 12% and 17.8% (Figure 9 (a)). This good agreement between the experimental and simulated textures was obtained for a value of Γ<sub>RX</sub> = 13.8 days. During the experiments, 17.8% strain was reached in 144 hours (~ 6 days). In the model simulations, this strain is reached in a shorter time (~ 130 hours), but the tertiary creep (maximum strain-rate) regime is reached between 10 to 20% strain, similar to the experimental observations (see also [27, 37]).

The macroscopic response predicted by R³iCe shows a clear softening associated with the
development of the texture. A maximum enhancement factor (EF) of 4.7 is obtained at ~ 19%
strain (black arrow, Figure 9), which corresponds to a girdle texture with a cone angle around
~ 45° to the compression axis. This softening is followed by a slight hardening as the girdle texture
slowly evolves toward a cone angle at ~ 35° to the compression axis.

The bulk strain associated with the maximum strain rate and the simulated enhancement
factor cannot be directly compared to that measured in Montagnat et al. 2015 [29] because
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Simulation of uniaxial compression along *Y* with no recrystallization Γ<sub>RX</sub> = +∞, *M*<sub>o</sub> = 0. a) Evolution with strain of the eigenvalues of the second-order orientation tensor (*a*<sub>1</sub> > *a*<sub>2</sub> > *a*<sub>3</sub>). Red dots are predictions from the VPSC model taken from [13]; b) evolution with strain of the macroscopic strain rate; c) evolution of the intensity of the texture with the colatitude *θ* over a section defined on the *Z* plane (dashed blue lines in d)). Data at 22%, 40% and 60% macrostrain are shown. The white dashed lines define the angle at ±45° from the compression axis *Y*; d) textures presented as pole figures for 22% (*i*), 40% (*ii*), and 60% (*iii*) macroscopic strain. All figures presenting texture data are color-coded using the same colorbar, which is displayed on the right (scales in multiples of a uniform distribution).</figcaption></figure>

at 7%, 12% and 17.8% (Figure 9 (a)). This good agreement between the experimental and simulated textures was obtained for a value of Γ<sub>RX</sub> = 13.8 days. During the experiments, 17.8% strain was reached in 144 hours (~ 6 days). In the model simulations, this strain is reached in a shorter time (~ 130 hours), but the tertiary creep (maximum strain-rate) regime is reached between 10 to 20% strain, similar to the experimental observations (see also [27, 37]).

The macroscopic response predicted by R³iCe shows a clear softening associated with the
development of the texture. A maximum enhancement factor (EF) of 4.7 is obtained at ~ 19%
strain (black arrow, Figure 9), which corresponds to a girdle texture with a cone angle around
~ 45° to the compression axis. This softening is followed by a slight hardening as the girdle texture
slowly evolves toward a cone angle at ~ 35° to the compression axis.

The bulk strain associated with the maximum strain rate and the simulated enhancement
factor cannot be directly compared to that measured in Montagnat et al. 2015 [29] because
these experiments were performed under constant load and the area of the sample section
evolved during the tests, whereas the model simulations were performed in a constant stress
configuration. The simulated strain rate may be corrected by assuming a constant sample volume
during the experiment, which leads to $\dot{\epsilon}_{sim}^{corr} = \dot{\epsilon}_{sim} \times (1 + \varepsilon)^3$. This corrected strain rate is plotted
in blue in Figure 9 (b). This brings the maximum strain rate to ~10% and the enhancement factor
to 3.0. These values are closer to the macroscopic response measured by Montagnat et al. [29],
which is characterized by an enhancement factor between 3.3 and 3.8.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Simulations for uniaxial compression creep along the *Y* direction with parameters from Table 3 and comparison with experimental data. a) Evolution with strain of the eigenvalues of the second-order orientation tensor (*a*₁ > *a*₂ > *a*₃). The blue dots correspond to experimental values from [29]; b) simulated macroscopic strain rate (black) and corrected simulated macroscopic strain rate (blue) to match the constant force boundary conditions instead of the constant stress one used in the simulation. EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude *θ* over a section defined on the *Z* plane (dashed blue lines in (d)). Data at 7%, 12% and 17.8% macrostrain are shown. The white dashed lines define the angle at ±45° from the compression axis *Y*; d) simulated textures (pole figures) for 7% (*i*), 12% (*ii*), and 17.8% (*iii*) macroscopic strains; e) measured textures for the same macroscopic strains [29]. All figures presenting texture data are color-coded using the same color bar (scales in multiples of a uniform distribution).</figcaption></figure>

5.3. *Uniaxial tension*

Experimental creep tests performed under uniaxial tension are scarce. Jacka et al. 1984 [27] pre-
sented one experiment performed at −3° C, up to 9.3% octahedral shear strain that corresponds
to the definition given in Jacka et al. 1984 [70], to 13.1% axial strain. At this strain level, they did
these experiments were performed under constant load and the area of the sample section
evolved during the tests, whereas the model simulations were performed in a constant stress
configuration. The simulated strain rate may be corrected by assuming a constant sample volume
during the experiment, which leads to $\dot{\epsilon}_{sim}^{corr} = \dot{\epsilon}_{sim} \times (1 + \varepsilon)^3$. This corrected strain rate is plotted
in blue in Figure 9 (b). This brings the maximum strain rate to ~10% and the enhancement factor
to 3.0. These values are closer to the macroscopic response measured by Montagnat et al. [29],
which is characterized by an enhancement factor between 3.3 and 3.8.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Simulations for uniaxial compression creep along the *Y* direction with parameters from Table 3 and comparison with experimental data. a) Evolution with strain of the eigenvalues of the second-order orientation tensor (*a*₁ > *a*₂ > *a*₃). The blue dots correspond to experimental values from [29]; b) simulated macroscopic strain rate (black) and corrected simulated macroscopic strain rate (blue) to match the constant force boundary conditions instead of the constant stress one used in the simulation. EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude *θ* over a section defined on the *Z* plane (dashed blue lines in (d)). Data at 7%, 12% and 17.8% macrostrain are shown. The white dashed lines define the angle at ±45° from the compression axis *Y*; d) simulated textures (pole figures) for 7% (*i*), 12% (*ii*), and 17.8% (*iii*) macroscopic strains; e) measured textures for the same macroscopic strains [29]. All figures presenting texture data are color-coded using the same color bar (scales in multiples of a uniform distribution).</figcaption></figure>

5.3. *Uniaxial tension*

Experimental creep tests performed under uniaxial tension are scarce. Jacka et al. 1984 [27] pre-
sented one experiment performed at −3° C, up to 9.3% octahedral shear strain that corresponds
to the definition given in Jacka et al. 1984 [70], to 13.1% axial strain. At this strain level, they did
not attain the quasi-constant strain rate typically of tertiary creep. The texture obtained, measured using the classical manual Rigsby stage technique with one orientation per grain, is presented in Figure 10. The texture is characterized by a small circle girdle with a mean half angle of 50.4° and a standard deviation of 15.3°. Although tertiary creep was not reached in this experiment, an enhancement factor of ~3, similar to that observed in compression under the same conditions, was measured.

The texture evolution and macroscopic response predicted by a simulation under uniaxial tension creep are shown in Figure 11. The texture evolves with increasing strain toward a girdle texture with a cone angle slightly around 45° to the extension direction. This texture appears to stabilize at ~5% axial strain. The distribution of the c-axis angle from the tension axis is given in Figure 12 for the simulated texture at 13% strain. It is similar to the one measured by Jacka et al. 1984 [27] (Figure 10). An enhancement factor of approximately 10 is simulated at 10% axial strain.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Texture obtained after 9% of octahedral tension strain (about 13% axial strain) performed at -3°C with an octahedral stress of 0.4 MPa. The star ⋆ in the center of the pole figure shows the direction of the tension axis. From Jacka et al. 1984 [27].</figcaption></figure>

### 5.4. Simple shear

The simulated texture evolution in simple shear is compared with observations made by Journaux et al. [31] for torsion experiments on isotropic granular ice performed at -7°C, under a constant torque corresponding to a maximum shear stress between 0.4 and 0.6 MPa.

The experimental textures are characterized by two submaxima, one almost perpendicular to the shear plane, called *M1*, and the second one, *M2*, initially at a low angle to the shear plane. Similar evolution was also observed in simple shear in [26, 41, 71]. As the shear strain increases, *M2* merges with *M1* to form a highly concentrated cluster texture characterized by a single maximum of c-axes normal to the shear plane. The simple shear simulations performed here reproduce this evolution relatively well, as shown in Figure 13.

By accounting for dynamic recrystallization, the R³iCe model is able to reproduce the evolution of both *M1* and *M2* maxima with strain. At low shear strains (γ ~ 0.2), *M1* concentrates normal to the shear plane, while *M2* follows the principal extension direction (white dashed line, Figure 13(c) bottom and Figure 13(d)). As strain increases, *M2* rotates toward *M1*. At a shear strain of 2, the texture is characterized by a slightly elongated single cluster normal to the shear plane as *M1* and *M2* merge. The evolution of the eigenvalues of the second-order orientation tensor, as observed in the experiments, is also well reproduced (Figure 13(a)).
not attain the quasi-constant strain rate typically of tertiary creep. The texture obtained, measured using the classical manual Rigsby stage technique with one orientation per grain, is presented in Figure 10. The texture is characterized by a small circle girdle with a mean half angle of 50.4° and a standard deviation of 15.3°. Although tertiary creep was not reached in this experiment, an enhancement factor of ~3, similar to that observed in compression under the same conditions, was measured.

The texture evolution and macroscopic response predicted by a simulation under uniaxial tension creep are shown in Figure 11. The texture evolves with increasing strain toward a girdle texture with a cone angle slightly around 45° to the extension direction. This texture appears to stabilize at ~5% axial strain. The distribution of the c-axis angle from the tension axis is given in Figure 12 for the simulated texture at 13% strain. It is similar to the one measured by Jacka et al. 1984 [27] (Figure 10). An enhancement factor of approximately 10 is simulated at 10% axial strain.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Texture obtained after 9% of octahedral tension strain (about 13% axial strain) performed at -3°C with an octahedral stress of 0.4 MPa. The star ⋆ in the center of the pole figure shows the direction of the tension axis. From Jacka et al. 1984 [27].</figcaption></figure>

### 5.4. Simple shear

The simulated texture evolution in simple shear is compared with observations made by Journaux et al. [31] for torsion experiments on isotropic granular ice performed at -7°C, under a constant torque corresponding to a maximum shear stress between 0.4 and 0.6 MPa.

The experimental textures are characterized by two submaxima, one almost perpendicular to the shear plane, called *M1*, and the second one, *M2*, initially at a low angle to the shear plane. Similar evolution was also observed in simple shear in [26, 41, 71]. As the shear strain increases, *M2* merges with *M1* to form a highly concentrated cluster texture characterized by a single maximum of c-axes normal to the shear plane. The simple shear simulations performed here reproduce this evolution relatively well, as shown in Figure 13.

By accounting for dynamic recrystallization, the R³iCe model is able to reproduce the evolution of both *M1* and *M2* maxima with strain. At low shear strains (γ ~ 0.2), *M1* concentrates normal to the shear plane, while *M2* follows the principal extension direction (white dashed line, Figure 13(c) bottom and Figure 13(d)). As strain increases, *M2* rotates toward *M1*. At a shear strain of 2, the texture is characterized by a slightly elongated single cluster normal to the shear plane as *M1* and *M2* merge. The evolution of the eigenvalues of the second-order orientation tensor, as observed in the experiments, is also well reproduced (Figure 13(a)).
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Simulations for uniaxial tension creep along the *Y* direction with parameters from Table 3 and comparison with experimental data. a) Evolution with strain of the eigenvalues of the second-order orientation tensor (*a*₁ > *a*₂ > *a*₃); b) simulated macroscopic strain rate. EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude *θ* over a section defined on the plane normal to *Z* (dashed blue lines in d)). Data at 5%, 13% and 15% axial macroscopic strains are shown. The white dashed lines define the angle at ±45° from the compression axis *Y*; d) simulated textures (represented as pole figures) for 5% (*i*), 13% (*ii*), and 15% (*iii*) axial macrostrain. All figures presenting texture data are color-coded using the same color bar (scales in multiples of a uniform distribution).</figcaption></figure>

In both experiments and simulations, the strain rate increased rapidly up to γ = 0.3 while the double-maxima texture develops (Figure 13(d)). It continues to increase, but at a progressively slowing rate as the texture evolves toward a single maximum normal to the shear plane, up to γ ~ 3 where it reaches a quasi-steady state. As explained in Section 5.2, the discontinuities in the strain-rate curve (Figure 13(b)) are technical artifacts imposed by the operating rules of GRICAD clusters.

In phase with Treverrow et al. 2012 [37], we define an enhancement factor for simple shear as the ratio between the strain rate during secondary creep (here it corresponds to the beginning of the run) and that obtained when the strongest simulated texture is reached. Using this definition, we obtain an enhancement factor (EF) of 6.7 that matches Treverrow et al. 2012 [37] data for a similar applied shear stress.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Simulations for uniaxial tension creep along the *Y* direction with parameters from Table 3 and comparison with experimental data. a) Evolution with strain of the eigenvalues of the second-order orientation tensor (*a*₁ > *a*₂ > *a*₃); b) simulated macroscopic strain rate. EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude *θ* over a section defined on the plane normal to *Z* (dashed blue lines in d)). Data at 5%, 13% and 15% axial macroscopic strains are shown. The white dashed lines define the angle at ±45° from the compression axis *Y*; d) simulated textures (represented as pole figures) for 5% (*i*), 13% (*ii*), and 15% (*iii*) axial macrostrain. All figures presenting texture data are color-coded using the same color bar (scales in multiples of a uniform distribution).</figcaption></figure>

In both experiments and simulations, the strain rate increased rapidly up to γ = 0.3 while the double-maxima texture develops (Figure 13(d)). It continues to increase, but at a progressively slowing rate as the texture evolves toward a single maximum normal to the shear plane, up to γ ~ 3 where it reaches a quasi-steady state. As explained in Section 5.2, the discontinuities in the strain-rate curve (Figure 13(b)) are technical artifacts imposed by the operating rules of GRICAD clusters.

In phase with Treverrow et al. 2012 [37], we define an enhancement factor for simple shear as the ratio between the strain rate during secondary creep (here it corresponds to the beginning of the run) and that obtained when the strongest simulated texture is reached. Using this definition, we obtain an enhancement factor (EF) of 6.7 that matches Treverrow et al. 2012 [37] data for a similar applied shear stress.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Distribution of the angle between the **c**-axis and the **y**-axis (tension axis) corresponding to pole figure (ii) in Figure 11 (d).</figcaption></figure>

5.5. *Testing the relative effects of deformation and recrystallization on texture evolution under high stresses*

The work of Qi et al. [30] provides some of the rare existing highly resolved textures obtained for temperatures and types of samples similar to those of the experiments simulated above, but under imposed displacement rates that induce higher stress conditions (see Table 3). The stress values reported in Table 3 correspond to the flow stress, i.e., the nearly constant stress reached after peak stress, at about 20% strain. To reach deviatoric stresses higher than 1 MPa without failure, a confining pressure of 10 MPa was applied.

To test the sensitivity of the dynamic recrystallization formulation to the level of stress, we performed simulations under creep conditions using the stress values from Qi et al. 2017 [30] experiments. As CTI-RX equations are relating the deviatoric stress to the deviatoric strain and are solved under the Stokes hypothesis (incompressibility), it is not necessary to numerically apply the confining pressure. Although boundary conditions differ between the simulations (macroscopic constant stress) and experiments (macroscopic constant displacement rate), we assumed the quasi-steady state behavior sampled by the flow stress measurements in the experiments to be equivalent to the quasi-steady state part of the tertiary creep in the creep simulations [72]. The stable textures expected at both stages can therefore be compared.

A comparison between the simulated and observed textures at different compressive flow stresses is shown in Figure 14. The experimental observations are well reproduced. Both measured and simulated textures show a transition from a girdle texture to a weak single maximum with increasing flow stress for a macroscopic strain of ~ 20%. To correctly fit the experimental data, a value of $\Gamma_{RX}$ of 4.6 days is used, which is lower than that used to reproduce the unconfinement compression creep tests of Montagnat et al. 2015 [29] (Table 3). This implies a higher contribution of recrystallization compared with deformation at higher stresses.

**6. Discussion**

The CTI-RX model presented here enables the simulation of texture evolution in polycrystalline ice, as shown in the simulations in compression, simple shear, and tension, including or not dynamic recrystallization. Without recrystallization, the texture evolution predicted by the model correctly reproduces the textures observed along deep ice cores where dynamic recrystallization is assumed to have a low impact.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Distribution of the angle between the **c**-axis and the **y**-axis (tension axis) corresponding to pole figure (ii) in Figure 11 (d).</figcaption></figure>

5.5. *Testing the relative effects of deformation and recrystallization on texture evolution under high stresses*

The work of Qi et al. [30] provides some of the rare existing highly resolved textures obtained for temperatures and types of samples similar to those of the experiments simulated above, but under imposed displacement rates that induce higher stress conditions (see Table 3). The stress values reported in Table 3 correspond to the flow stress, i.e., the nearly constant stress reached after peak stress, at about 20% strain. To reach deviatoric stresses higher than 1 MPa without failure, a confining pressure of 10 MPa was applied.

To test the sensitivity of the dynamic recrystallization formulation to the level of stress, we performed simulations under creep conditions using the stress values from Qi et al. 2017 [30] experiments. As CTI-RX equations are relating the deviatoric stress to the deviatoric strain and are solved under the Stokes hypothesis (incompressibility), it is not necessary to numerically apply the confining pressure. Although boundary conditions differ between the simulations (macroscopic constant stress) and experiments (macroscopic constant displacement rate), we assumed the quasi-steady state behavior sampled by the flow stress measurements in the experiments to be equivalent to the quasi-steady state part of the tertiary creep in the creep simulations [72]. The stable textures expected at both stages can therefore be compared.

A comparison between the simulated and observed textures at different compressive flow stresses is shown in Figure 14. The experimental observations are well reproduced. Both measured and simulated textures show a transition from a girdle texture to a weak single maximum with increasing flow stress for a macroscopic strain of ~ 20%. To correctly fit the experimental data, a value of $\Gamma_{RX}$ of 4.6 days is used, which is lower than that used to reproduce the unconfinement compression creep tests of Montagnat et al. 2015 [29] (Table 3). This implies a higher contribution of recrystallization compared with deformation at higher stresses.

**6. Discussion**

The CTI-RX model presented here enables the simulation of texture evolution in polycrystalline ice, as shown in the simulations in compression, simple shear, and tension, including or not dynamic recrystallization. Without recrystallization, the texture evolution predicted by the model correctly reproduces the textures observed along deep ice cores where dynamic recrystallization is assumed to have a low impact.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Simulations for a simple shear creep test (shear plane normal to the *Y*-axis and shear direction parallel to *X*, $\frac{\partial u_x}{\partial y} < 0$) with parameters from Table 3 and comparison with experimental data from [31]. a) Evolution with strain of the eigenvalues of the second-order orientation tensor ($a_1 > a_2 > a_3$). The blue dots correspond to the experimental values from Journaux et al. 2019 [31]; b) simulated macroscopic strain rate. EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude $\theta$ over a section defined on the *Z* plane (dashed blue lines in d)). Data at $\gamma = 0.2, 0.42$, and 1.96 are shown. The white dashed line follows the orientation of the principal extension axis; d) simulated textures (pole figures, upper hemisphere) at $\gamma = 0.2$ (i), $\gamma = 0.42$ (ii) and $\gamma = 1.96$ (iii); e) measured textures for the same macroscopic shear strains [31]. All figures presenting texture data are color-coded using the same colorbar (scales in multiples of a uniform distribution).</figcaption></figure>

The novelty of the present model and its implementation is based on the simulation of the evolution of the texture (orientations of c-axes) based on a new formulation that (i) introduces an attractor controlled by the local stress field to reproduce the impact of DRX mechanisms, (ii) is coupled with the continuous transverse isotropic (CTI) law to account for strain and stress field heterogeneities at the crystal scale arising from the crystal viscoplastic anisotropy, and (iii) solves these coupled equations using a semi-implicit numerical scheme in a full-field implementation using FEM, ensuring consistency between texture evolution and mechanical behavior. The first two ingredients are derived from experimentally based knowledge of dynamic recrystallization processes.

Indeed, recent experimental observations have (i) clearly documented the role of the crystal viscoplastic anisotropy in generating local stresses and strains markedly different from the macroscopic ones [3, 29, 66] and (ii) revealed that the formation of new orientations (nucleation) during dynamic recrystallization occurs through bulging (heterogeneous grain boundary migra-
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Simulations for a simple shear creep test (shear plane normal to the Y-axis and shear direction parallel to *X*, $\frac{\partial u_x}{\partial y} < 0$) with parameters from Table 3 and comparison with experimental data from [31]. a) Evolution with strain of the eigenvalues of the second-order orientation tensor ($a_1 > a_2 > a_3$). The blue dots correspond to the experimental values from Journaux et al. 2019 [31]; b) simulated macroscopic strain rate. EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude $\theta$ over a section defined on the Z plane (dashed blue lines in d)). Data at $\gamma = 0.2, 0.42$, and 1.96 are shown. The white dashed line follows the orientation of the principal extension axis; d) simulated textures (pole figures, upper hemisphere) at $\gamma = 0.2$ (i), $\gamma = 0.42$ (ii) and $\gamma = 1.96$ (iii); e) measured textures for the same macroscopic shear strains [31]. All figures presenting texture data are color-coded using the same colorbar (scales in multiples of a uniform distribution).</figcaption></figure>

The novelty of the present model and its implementation is based on the simulation of the evolution of the texture (orientations of c-axes) based on a new formulation that (i) introduces an attractor controlled by the local stress field to reproduce the impact of DRX mechanisms, (ii) is coupled with the continuous transverse isotropic (CTI) law to account for strain and stress field heterogeneities at the crystal scale arising from the crystal viscoplastic anisotropy, and (iii) solves these coupled equations using a semi-implicit numerical scheme in a full-field implementation using FEM, ensuring consistency between texture evolution and mechanical behavior. The first two ingredients are derived from experimentally based knowledge of dynamic recrystallization processes.

Indeed, recent experimental observations have (i) clearly documented the role of the crystal viscoplastic anisotropy in generating local stresses and strains markedly different from the macroscopic ones [3, 29, 66] and (ii) revealed that the formation of new orientations (nucleation) during dynamic recrystallization occurs through bulging (heterogeneous grain boundary migra-
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Simulation of the uniaxial compression creep test along *Y* at various stresses. Comparison between R³iCe predictions and experiments by Qi et al. [30]. (a) Textures from [30] obtained by EBSD for various strain rates leading to different quasi-steady state compressive flow stresses at macroscopic strains of ~ 0.2% (exact strain is indicated at the top left of each pole figure). (b) Simulated textures obtained by creep under the same compressive stresses and up to the same macroscopic strains as in (a).</figcaption></figure>

tion) and polygonization (formation of new grains boundaries by organization of dislocations within a grain) [3, 33, 63, 66, 73]. These two processes are controlled by the local strain and stress fields and result in recrystallized orientations that are strongly related to, but slightly different from, the parent grain orientations. The last important experimental observation is that the development of recrystallization textures leads to a weakening of the polycrystal mechanical response. All these observations justify the choice of a continuous formulation for the *c*-axis orientation evolution (Equation (22)), with a targeted recrystallized orientation that is calculated to maximize the local basal resolved shear stress (Equation (26)).

Comparison with tertiary creep experiments performed in compression [29, 37, 39], in simple shear [31, 37, 41] and in tension [27] supports the fact that the R³iCe model can reproduce texture evolution in conditions similar to those of laboratory experiments (T = -5 to -10°C, S ~ 1 MPa). For this, a single parameter is adjusted: the recrystallization parameter Γ<sub>RX</sub>, which weights the relative contribution of deformation and recrystallization to the rotation of the *c*-axes. As expected, under similar stress conditions, a lower temperature requires a higher Γ<sub>RX</sub> value (and therefore a lower weight of dynamic recrystallization) to correctly reproduce the experimental texture evolution. A noticeable result is that the two-maxima texture measured in simple shear (in natural and laboratory ice) is correctly reproduced here, together with the kinetics of its evolution with strain toward a single-maximum cluster.

For all conditions simulated with R³iCe, we show that the CTI-RX model reproduces accurately the softening at the transition from secondary to tertiary creep. The kinetics of the transition and
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Simulation of the uniaxial compression creep test along *Y* at various stresses. Comparison between R³iCe predictions and experiments by Qi et al. [30]. (a) Textures from [30] obtained by EBSD for various strain rates leading to different quasi-steady state compressive flow stresses at macroscopic strains of ~ 0.2% (exact strain is indicated at the top left of each pole figure). (b) Simulated textures obtained by creep under the same compressive stresses and up to the same macroscopic strains as in (a).</figcaption></figure>

tion) and polygonization (formation of new grains boundaries by organization of dislocations within a grain) [3, 33, 63, 66, 73]. These two processes are controlled by the local strain and stress fields and result in recrystallized orientations that are strongly related to, but slightly different from, the parent grain orientations. The last important experimental observation is that the development of recrystallization textures leads to a weakening of the polycrystal mechanical response. All these observations justify the choice of a continuous formulation for the *c*-axis orientation evolution (Equation (22)), with a targeted recrystallized orientation that is calculated to maximize the local basal resolved shear stress (Equation (26)).

Comparison with tertiary creep experiments performed in compression [29, 37, 39], in simple shear [31, 37, 41] and in tension [27] supports the fact that the R³iCe model can reproduce texture evolution in conditions similar to those of laboratory experiments (T = -5 to -10°C, S ~ 1 MPa). For this, a single parameter is adjusted: the recrystallization parameter Γ<sub>RX</sub>, which weights the relative contribution of deformation and recrystallization to the rotation of the *c*-axes. As expected, under similar stress conditions, a lower temperature requires a higher Γ<sub>RX</sub> value (and therefore a lower weight of dynamic recrystallization) to correctly reproduce the experimental texture evolution. A noticeable result is that the two-maxima texture measured in simple shear (in natural and laboratory ice) is correctly reproduced here, together with the kinetics of its evolution with strain toward a single-maximum cluster.

For all conditions simulated with R³iCe, we show that the CTI-RX model reproduces accurately the softening at the transition from secondary to tertiary creep. The kinetics of the transition and
the resulting enhancement factors, when departing from an isotropic texture, are also correctly simulated, compared with experimental results performed under similar conditions [22, 27, 29, 36, 37]. Because the present model does not consider any softening mechanism associated with dislocation interactions or grain boundary migration, this comparison seems to confirm that texture-induced viscoplastic anisotropy explains most of the mechanical softening measured during tertiary creep. Therefore, softening associated with dislocation annihilation and grain boundary migration during dynamic recrystallization appears to play a secondary role.

Simulations at higher uniaxial compressive stresses (from 1 to 4 MPa) have been performed to compare with the experiments of Qi et al. 2017 [30]. Contrary to the previously mentioned experiments, these experiments were performed under a range of imposed strain-rate conditions, which resulted in different flow stresses. A confining pressure of 10 MPa was applied to prevent failure at the relatively high strain rates considered. The observation of more clustered textures at higher flow stress implies that higher strain rates result in a larger impact of deformation relative to DRX on texture evolution. R³iCe run under imposed stress without confining pressure can correctly represent texture clustering under stress. Indeed, experimental textures measured at 20% strain are satisfactorily reproduced by R³iCe providing that Γ<sub>RX</sub> is given a relatively low value (4.6 days for T = -10°C) compared to the one determined for unconfined creep compression experiments (13.8 days for T = -5°C). This low value of Γ<sub>RX</sub> indicates that, in the simulations, a high amount of recrystallization is necessary to reproduce the experimental textures and counterbalance the impact of deformation on the c-axes rotation. We interpret this apparent contradiction as due to the activation of brittle-ductile processes in the experiments. This assumption is based on observations by Kalifa et al. [74], which performed triaxial tests under the same range of strain rate (≈ 10<sup>-5</sup> – 10<sup>4</sup> s<sup>-1</sup>) than Qi et al. 2017 [30], but under variable confining pressure from 0 to 10 MPa. By careful observations made immediately after the peak stress, they estimated the density of cracks in the microstructure and revealed that micro-cracks were visible up to 10 MPa confining pressure. Therefore, microcracking is also expected to have occurred in Qi et al. 2017 [30] experiments. Strong stress concentration occurs at crack tips, and Chauve et al. 2017 [73] have recently shown that, for ice, this stress concentration can be released through DRX. The R³iCe formulation does not account for microcracking or the effect of confining pressure. Therefore, the enhancement of recrystallization due to the plastic energy available around the crack tips needs to be reproduced by increasing the recrystallization rate (Γ<sub>RX</sub>).

Recent models of texture evolution during dynamic recrystallization have been proposed by Rathmann et al. 2021 [54] and Richard et al. 2021,2022 [47, 75], aiming to account for texture-induced anisotropy of ice in large-scale models. Their formulation is based on an orientation distribution function (ODF) evolution equation in which dynamic recrystallization is represented by two terms: one term that accounts for the production of orientations based on a given parameter, the *deformability*, that is calculated using the Sachs constant stress assumption, and a second term that accounts for the diffusion of orientations. This approach is not intended to predict the mechanical response evolution due to texture development, but the texture evolution equation can be compared with that in the CTI-RX formulation. The *deformability* term can be obtained from Equation (22) by assuming homogeneous stresses (Sachs model). The diffusive term, which is necessary to reproduce the texture dispersion effect owing to recrystallization, is spontaneously obtained in full-field approaches such as R³iCe, due to the heterogeneous stress field. Although these models accurately predict steady-state textures for different strain geometries, the predictions for transient texture evolution are less accurate than that from R³iCe. This highlights the fundamental role of stress heterogeneities in texture evolution during the transition between secondary and tertiary creep (from 1 to ~ 15% strain).

The major strength of [47, 54, 75] formulations is their numerical efficiency, which makes them suitable for large-scale modeling (glaciers, ice sheets), provided that transient behaviors
the resulting enhancement factors, when departing from an isotropic texture, are also correctly simulated, compared with experimental results performed under similar conditions [22, 27, 29, 36, 37]. Because the present model does not consider any softening mechanism associated with dislocation interactions or grain boundary migration, this comparison seems to confirm that texture-induced viscoplastic anisotropy explains most of the mechanical softening measured during tertiary creep. Therefore, softening associated with dislocation annihilation and grain boundary migration during dynamic recrystallization appears to play a secondary role.

Simulations at higher uniaxial compressive stresses (from 1 to 4 MPa) have been performed to compare with the experiments of Qi et al. 2017 [30]. Contrary to the previously mentioned experiments, these experiments were performed under a range of imposed strain-rate conditions, which resulted in different flow stresses. A confining pressure of 10 MPa was applied to prevent failure at the relatively high strain rates considered. The observation of more clustered textures at higher flow stress implies that higher strain rates result in a larger impact of deformation relative to DRX on texture evolution. R³iCe run under imposed stress without confining pressure can correctly represent texture clustering under stress. Indeed, experimental textures measured at 20% strain are satisfactorily reproduced by R³iCe providing that Γ<sub>RX</sub> is given a relatively low value (4.6 days for T = -10°C) compared to the one determined for unconfined creep compression experiments (13.8 days for T = -5°C). This low value of Γ<sub>RX</sub> indicates that, in the simulations, a high amount of recrystallization is necessary to reproduce the experimental textures and counterbalance the impact of deformation on the c-axes rotation. We interpret this apparent contradiction as due to the activation of brittle-ductile processes in the experiments. This assumption is based on observations by Kalifa et al. [74], which performed triaxial tests under the same range of strain rate (≈ 10<sup>-5</sup> – 10<sup>4</sup> s<sup>-1</sup>) than Qi et al. 2017 [30], but under variable confining pressure from 0 to 10 MPa. By careful observations made immediately after the peak stress, they estimated the density of cracks in the microstructure and revealed that micro-cracks were visible up to 10 MPa confining pressure. Therefore, microcracking is also expected to have occurred in Qi et al. 2017 [30] experiments. Strong stress concentration occurs at crack tips, and Chauve et al. 2017 [73] have recently shown that, for ice, this stress concentration can be released through DRX. The R³iCe formulation does not account for microcracking or the effect of confining pressure. Therefore, the enhancement of recrystallization due to the plastic energy available around the crack tips needs to be reproduced by increasing the recrystallization rate (Γ<sub>RX</sub>).

Recent models of texture evolution during dynamic recrystallization have been proposed by Rathmann et al. 2021 [54] and Richard et al. 2021,2022 [47, 75], aiming to account for texture-induced anisotropy of ice in large-scale models. Their formulation is based on an orientation distribution function (ODF) evolution equation in which dynamic recrystallization is represented by two terms: one term that accounts for the production of orientations based on a given parameter, the *deformability*, that is calculated using the Sachs constant stress assumption, and a second term that accounts for the diffusion of orientations. This approach is not intended to predict the mechanical response evolution due to texture development, but the texture evolution equation can be compared with that in the CTI-RX formulation. The *deformability* term can be obtained from Equation (22) by assuming homogeneous stresses (Sachs model). The diffusive term, which is necessary to reproduce the texture dispersion effect owing to recrystallization, is spontaneously obtained in full-field approaches such as R³iCe, due to the heterogeneous stress field. Although these models accurately predict steady-state textures for different strain geometries, the predictions for transient texture evolution are less accurate than that from R³iCe. This highlights the fundamental role of stress heterogeneities in texture evolution during the transition between secondary and tertiary creep (from 1 to ~ 15% strain).

The major strength of [47, 54, 75] formulations is their numerical efficiency, which makes them suitable for large-scale modeling (glaciers, ice sheets), provided that transient behaviors
are not of concern. In such a context, the R³iCe model could be a good candidate for providing
constraints to develop parameterizations for the temporal evolution, destined to inexpensive
formulations.

The CTI-RX model is not limited to the loading conditions prescribed here or to polycrystalline ice. Its formulation is generic and can be extended to other polycrystalline materials, provided that their single-crystal viscoplastic behavior can be adequately described using the continuous transverse isotropic law defined in Equation (17). For instance, the versatility of the CTI-RX model could enable the analysis of the mechanical behavior of materials such as e.g. magnesium and quartz, which also show strong viscoplastic anisotropy with hexagonal symmetry.

**7. Conclusions**

In this paper, we present a new formulation, CTI-RX, to model the effect of dynamic recrystal-
lization on the texture evolution of polycrystalline materials, provided that their viscoplastic be-
havior can be described using a continuous transverse isotropic (CTI) law. This formulation is
validated in the context of polycrystalline ice.

The integration of an orientation attractor, denoted as $c_0$ along with an anisotropic flow law in full-field resolution allows for the accurate replication of texture evolution during dynamic recrystallization in tertiary creep under compression, tension, and simple shear. For the latter, an accurate reproduction of the two-maxima texture and its evolution with strain is obtained. This is made possible by defining $c_0$ such as to maximize the local resolved shear stress, thus providing a physically based formulation for recrystallization-induced c-axes rotation. This formulation is as simple as possible, given our knowledge of the mechanisms.

The R³iCe model is the result of the implementation of the CTI-RX formulation in a finite
element framework. In the model, the texture results from a balance between c-axis rotation due
to viscoplastic deformation and dynamic recrystallization, which is controlled by Γ<sub>RX</sub>, the only
tuning parameter of the model. The R³iCe model can therefore be used to constrain the impact
of experimental conditions, such as pressure and temperature, on the recrystallization kinetics.

Accurate reproduction of the textures in the R³iCe model leads to a good prediction of the
mechanical softening associated with dynamic recrystallization. This confirms that texture-
induced viscoplastic anisotropy may explain most of the mechanical softening observed during
tertiary creep and suggests that recovery and grain boundary migration play a secondary role in
this softening.

Under conditions of high stress, where confining pressure is required in the laboratory to pre-
vent failure, we propose that local fracturing at the crystal scale likely enhances dynamic recrys-
tallization processes and accelerates texture evolution with strain. This is consistent with the
lower Γ<sub>RX</sub> required to fit the texture measured during the confined compression experiments [30]
relative to that used to fit the low stress unconfined compression creep experiments [29], despite
the higher experimental temperature of the latter (−5° C vs. −10° C, cf. Table 3). These results
provide an additional illustration of the ability of the R³iCe model to help resolve open questions
about ice deformation behavior.

The computational cost of the R³iCe full-field model is not adapted for direct implementation
in large-scale flow modeling frameworks. The model is nonetheless highly valuable for constrain-
ing the parameterization required to properly account for texture evolution and its impact on the
mechanical response in complex or changing boundary conditions such as those prevailing at
the bottom of deep ice cores, ice streams, and glaciers.
are not of concern. In such a context, the R³iCe model could be a good candidate for providing constraints to develop parameterizations for the temporal evolution, destined to inexpensive formulations.

The CTI-RX model is not limited to the loading conditions prescribed here or to polycrystalline ice. Its formulation is generic and can be extended to other polycrystalline materials, provided that their single-crystal viscoplastic behavior can be adequately described using the continuous transverse isotropic law defined in Equation (17). For instance, the versatility of the CTI-RX model could enable the analysis of the mechanical behavior of materials such as e.g. magnesium and quartz, which also show strong viscoplastic anisotropy with hexagonal symmetry.

## 7. Conclusions

In this paper, we present a new formulation, CTI-RX, to model the effect of dynamic recrystallization on the texture evolution of polycrystalline materials, provided that their viscoplastic behavior can be described using a continuous transverse isotropic (CTI) law. This formulation is validated in the context of polycrystalline ice.

The integration of an orientation attractor, denoted as $c_0$ along with an anisotropic flow law in full-field resolution allows for the accurate replication of texture evolution during dynamic recrystallization in tertiary creep under compression, tension, and simple shear. For the latter, an accurate reproduction of the two-maxima texture and its evolution with strain is obtained. This is made possible by defining $c_0$ such as to maximize the local resolved shear stress, thus providing a physically based formulation for recrystallization-induced c-axes rotation. This formulation is as simple as possible, given our knowledge of the mechanisms.

The R³iCe model is the result of the implementation of the CTI-RX formulation in a finite element framework. In the model, the texture results from a balance between c-axis rotation due to viscoplastic deformation and dynamic recrystallization, which is controlled by $\Gamma_{RX}$, the only tuning parameter of the model. The R³iCe model can therefore be used to constrain the impact of experimental conditions, such as pressure and temperature, on the recrystallization kinetics.

Accurate reproduction of the textures in the R³iCe model leads to a good prediction of the mechanical softening associated with dynamic recrystallization. This confirms that texture-induced viscoplastic anisotropy may explain most of the mechanical softening observed during tertiary creep and suggests that recovery and grain boundary migration play a secondary role in this softening.

Under conditions of high stress, where confining pressure is required in the laboratory to prevent failure, we propose that local fracturing at the crystal scale likely enhances dynamic recrystallization processes and accelerates texture evolution with strain. This is consistent with the lower $\Gamma_{RX}$ required to fit the texture measured during the confined compression experiments [30] relative to that used to fit the low stress unconfined compression creep experiments [29], despite the higher experimental temperature of the latter (−5° C vs. −10° C, cf. Table 3). These results provide an additional illustration of the ability of the R³iCe model to help resolve open questions about ice deformation behavior.

The computational cost of the R³iCe full-field model is not adapted for direct implementation in large-scale flow modeling frameworks. The model is nonetheless highly valuable for constraining the parameterization required to properly account for texture evolution and its impact on the mechanical response in complex or changing boundary conditions such as those prevailing at the bottom of deep ice cores, ice streams, and glaciers.
**Conflicts of interest**

The authors declare that the research was conducted in the absence of any commercial or
financial relationships that could be understood as a potential conflict of interest.

**Acknowledgements**

The authors are thankful to Chao Qi for quickly providing the EBSD data from [30] to test R³iCe abilities during deformation under confining pressure and high stresses.

All computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by the Grenoble research community.

The authors extend their sincere gratitude to the anonymous reviewers for their meticulous examination and valuable feedback, which significantly improved the quality and clarity of this manuscript.

**Authors' contributions**

Thomas Chauve conducted most of the study, developed the model, performed the numerical implementation of R³iCe, the simulations, their analyses and interpretations, and drafted the manuscript.

Maurine Montagnat elaborated on the scientific question, participated in the development of the equations, and followed the entire implementation of the model. She participated in the analyses and interpretation of the results. She is the co-principal author of the paper with TC.

Véronique Dansereau participated in the elaboration of the scientific question and co-developed the model with TC. She participated in writing the manuscript.

Pierre Saramito helped with the development of the model from the theoretical aspect to the technical solution.

Kevin Fourteau proposed, with TC and MM, the use of an attractor to model dynamic recrystallization.

Andréa Tommasi obtained funding and contributed to the evaluation of the model and writing of the final versions of the manuscript.

All authors have read and approved the final manuscript.

**Open Access**

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution, and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third-party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a full copy of this license, please visit http://creativecommons.org/licenses/by/4.0/.

The R³iCe code is available under the GPL-v3 license. Gitlab repository https://gricad-gitlab.univ-grenoble-alpes.fr/mecaiceige/tools/ice-polycrystal-models/rheolef_cti. The associated documentation is available at https://mecaiceige.gricad-pages.univ-grenoble-alpes.fr/tools/ice-polycrystal-models/ipms_documentation/Intro.html
**Conflicts of interest**

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be understood as a potential conflict of interest.

**Acknowledgements**

The authors are thankful to Chao Qi for quickly providing the EBSD data from [30] to test R³iCe abilities during deformation under confining pressure and high stresses.

All computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by the Grenoble research community.

The authors extend their sincere gratitude to the anonymous reviewers for their meticulous examination and valuable feedback, which significantly improved the quality and clarity of this manuscript.

**Authors' contributions**

Thomas Chauve conducted most of the study, developed the model, performed the numerical implementation of R³iCe, the simulations, their analyses and interpretations, and drafted the manuscript.

Maurine Montagnat elaborated on the scientific question, participated in the development of the equations, and followed the entire implementation of the model. She participated in the analyses and interpretation of the results. She is the co-principal author of the paper with TC.

Véronique Dansereau participated in the elaboration of the scientific question and co-developed the model with TC. She participated in writing the manuscript.

Pierre Saramito helped with the development of the model from the theoretical aspect to the technical solution.

Kevin Fourteau proposed, with TC and MM, the use of an attractor to model dynamic recrystallization.

Andréa Tommasi obtained funding and contributed to the evaluation of the model and writing of the final versions of the manuscript.

All authors have read and approved the final manuscript.

**Open Access**

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution, and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third-party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a full copy of this license, please visit http://creativecommons.org/licenses/by/4.0/.

The R³iCe code is available under the GPL-v3 license. Gitlab repository https://gricad-gitlab.univ-grenoble-alpes.fr/mecaiceige/tools/ice-polycrystal-models/rheolef_cti. The associated documentation is available at https://mecaiceige.gricad-pages.univ-grenoble-alpes.fr/tools/ice-polycrystal-models/ipms_documentation/Intro.html
# Appendix A. Demonstration of δᵢ identification

## A.1. Constitutive equation CTI and single crystal behavior

The constitutive equation is

$$ \mathbf{D} = \left( \frac{\partial \phi_S^{(1)}}{\partial \mathbf{S}} \right)^D = 2\delta_1 \mathbf{S} + 2\delta_2 \mathbf{M}^D \mathrm{Tr}(\mathbf{MS}) + \delta_3 (\mathbf{MS} + \mathbf{SM})^D \quad (A1) $$

Single crystal anisotropy can be described by three parameters derived from experimental data:

$\psi_1$ is the fluidity for parallel shearing to the basal plane:

$$ gD_{13} = \frac{\psi_1}{2} gS_{13}, \quad gD_{23} = \frac{\psi_1}{2} gS_{23} \quad (A2) $$

$\beta$ is the viscosity ratio between shear parallel to the basal plane $\eta$ and the viscosity for shear within the basal plane:

$$ gD_{12} = \beta \frac{\psi_1}{2} gS_{12} \quad (A3) $$

$\gamma$ is the viscosity ratio between compression (or traction) Along the c axis $\eta_c$ and compression (or traction) in one direction $g e_r$ within the basal plane:

$$ gS_{33} = 2\eta_c gD_{33}, \quad \gamma^g S_{rr} = 2\eta_c^g D_{rr} \quad (A4) $$

The same compressive (or extensive) stress parallel and normal to the c-axis $S = gS_{33} = gS_{rr}$ results therefore in different strain rates:

$$ gD_{rr} = \gamma^g D_{33} \quad (A5) $$

## A.2. δᵢ identification

To identify the parameters ($\delta_1, \delta_2, \delta_3$) in the constitutive Equation (A1), we developed it and compared it with the single crystal behavior.

Using $\mathbf{c} = (0, 0, 1)$

$$ \mathbf{M} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \mathbf{M}^D = \frac{1}{3} \begin{bmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 2 \end{bmatrix}, \quad \mathrm{Tr}(\mathbf{M}^g \mathbf{D}) = gD_{33} \quad (A6) $$

(1) $\delta_1$ is given by applying stress

$$ g\mathbf{S} = \begin{bmatrix} 0 & gS_{12} & 0 \\ gS_{12} & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$

Then

$$ \mathrm{Tr}(\mathbf{MS}) = 0, \quad (\mathbf{MS} + \mathbf{SM})^D = 0 \quad (A7) $$

The constitutive equation gives

$$ gD_{12} = 2\delta_1 gS_{12} \quad (A8) $$

Therefore:

$$ \delta_1 = \frac{\psi_1 \beta}{4} \quad (A9) $$
# Appendix A. Demonstration of δᵢ identification

## A.1. Constitutive equation CTI and single crystal behavior

The constitutive equation is

$$ \mathbf{D} = \left( \frac{\partial \phi_S^{(1)}}{\partial \mathbf{S}} \right)^D = 2\delta_1 \mathbf{S} + 2\delta_2 \mathbf{M}^D \mathrm{Tr}(\mathbf{MS}) + \delta_3 (\mathbf{MS} + \mathbf{SM})^D \quad (A1) $$

Single crystal anisotropy can be described by three parameters derived from experimental data:

$\psi_1$ is the fluidity for parallel shearing to the basal plane:

$$ gD_{13} = \frac{\psi_1}{2} gS_{13}, \quad gD_{23} = \frac{\psi_1}{2} gS_{23} \quad (A2) $$

$\beta$ is the viscosity ratio between shear parallel to the basal plane $\eta$ and the viscosity for shear within the basal plane:

$$ gD_{12} = \beta \frac{\psi_1}{2} gS_{12} \quad (A3) $$

$\gamma$ is the viscosity ratio between compression (or traction) Along the c axis $\eta_c$ and compression (or traction) in one direction $g e_r$ within the basal plane:

$$ gS_{33} = 2\eta_c gD_{33}, \quad \gamma^g S_{rr} = 2\eta_c^g D_{rr} \quad (A4) $$

The same compressive (or extensive) stress parallel and normal to the c-axis $S = gS_{33} = gS_{rr}$ results therefore in different strain rates:

$$ gD_{rr} = \gamma^g D_{33} \quad (A5) $$

## A.2. δᵢ identification

To identify the parameters ($\delta_1, \delta_2, \delta_3$) in the constitutive Equation (A1), we developed it and compared it with the single crystal behavior.

Using $\mathbf{c} = (0, 0, 1)$

$$ \mathbf{M} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \mathbf{M}^D = \frac{1}{3} \begin{bmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 2 \end{bmatrix}, \quad \mathrm{Tr}(\mathbf{M}^g \mathbf{D}) = gD_{33} \quad (A6) $$

(1) $\delta_1$ is given by applying stress

$$ g\mathbf{S} = \begin{bmatrix} 0 & gS_{12} & 0 \\ gS_{12} & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$

Then

$$ \mathrm{Tr}(\mathbf{MS}) = 0, \quad (\mathbf{MS} + \mathbf{SM})^D = 0 \quad (A7) $$

The constitutive equation gives

$$ gD_{12} = 2\delta_1 gS_{12} \quad (A8) $$

Therefore:

$$ \delta_1 = \frac{\psi_1 \beta}{4} \quad (A9) $$
(2) $\delta_3$ is given by applying stress

$$g\mathbf{S} = \begin{bmatrix} 0 & 0 & gS_{13} \\ 0 & 0 & 0 \\ gS_{13} & 0 & 0 \end{bmatrix}$$

Then

$$Tr(\mathbf{MS}) = 0, \quad (\mathbf{MS} + \mathbf{SM})^D = \mathbf{S} \qquad (A10)$$

The constitutive equation gives

$$gD_{13} = (2\delta_1 + \delta_3)gS_{13} \qquad (A11)$$

Therefore

$$\delta_3 = \frac{\psi_1}{2} (1 - \beta) \qquad (A12)$$

(3) $\delta_2$ identification is less trivial and presented in more detail.

* First, uniaxial stress $\tilde{\sigma}$ is applied along the c-axis. The stress tensor is $g\sigma = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \tilde{\sigma} \end{bmatrix}$ and the associated deviatoric stress tensor is

$$g\mathbf{S} = \begin{bmatrix} -\frac{1}{3}\tilde{\sigma} & 0 & 0 \\ 0 & -\frac{1}{3}\tilde{\sigma} & 0 \\ 0 & 0 & \frac{2}{3}\tilde{\sigma} \end{bmatrix}$$

$$- \mathrm{Tr}(\mathbf{M}^g \mathbf{S}) = \frac{2\tilde{\sigma}}{3}$$

$$-(\mathbf{M}^g \mathbf{S} + g\mathbf{SM})^D = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \frac{4}{3}\tilde{\sigma} \end{bmatrix}^D = \begin{bmatrix} -\frac{4}{9}\tilde{\sigma} & 0 & 0 \\ 0 & -\frac{4}{9}\tilde{\sigma} & 0 \\ 0 & 0 & \frac{8}{9}\tilde{\sigma} \end{bmatrix}$$

- using the constitutive equation, we obtain for $gD_{33}$:

$$gD_{33} = \delta_1 \frac{4}{3}\tilde{\sigma} + \delta_2 \frac{8}{9}\tilde{\sigma} + \delta_3 \frac{8}{9}\tilde{\sigma}$$

* Then a uniaxial stress $\tilde{\sigma}$ is applied normal to the c-axis, for instance in the direction 1. The stress tensor is

$$g\sigma = \begin{bmatrix} \tilde{\sigma} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}$$

and the associated diatoric stress tensor is

$$g\mathbf{S} = \begin{bmatrix} \frac{2}{3}\tilde{\sigma} & 0 & 0 \\ 0 & -\frac{1}{3}\tilde{\sigma} & 0 \\ 0 & 0 & -\frac{1}{3}\tilde{\sigma} \end{bmatrix}$$

$$- \mathrm{Tr}(\mathbf{M}^g \mathbf{S}) = -\frac{\tilde{\sigma}}{3}$$

$$-(\mathbf{M}^g \mathbf{S} + g\mathbf{SM})^D = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -\frac{2}{3}\tilde{\sigma} \end{bmatrix}^D = \begin{bmatrix} \frac{2}{9}\tilde{\sigma} & 0 & 0 \\ 0 & \frac{2}{9}\tilde{\sigma} & 0 \\ 0 & 0 & \frac{4}{9}\tilde{\sigma} \end{bmatrix}$$

- using the constitutive equation, we obtain for $gD_{11}$:

$$gD_{11} = \delta_1 \frac{4}{3}\tilde{\sigma} + \delta_2 \frac{2}{9}\tilde{\sigma} + \delta_3 \frac{2}{9}\tilde{\sigma}$$
(2) $\delta_3$ is given by applying stress

$$g\mathbf{S} = \begin{bmatrix} 0 & 0 & gS_{13} \\ 0 & 0 & 0 \\ gS_{13} & 0 & 0 \end{bmatrix}$$

Then

$$Tr(\mathbf{MS}) = 0, \quad (\mathbf{MS} + \mathbf{SM})^D = \mathbf{S} \qquad (A10)$$

The constitutive equation gives

$$gD_{13} = (2\delta_1 + \delta_3)gS_{13} \qquad (A11)$$

Therefore

$$\delta_3 = \frac{\psi_1}{2} (1 - \beta) \qquad (A12)$$

(3) $\delta_2$ identification is less trivial and presented in more detail.

* First, uniaxial stress $\tilde{\sigma}$ is applied along the c-axis. The stress tensor is $g\sigma = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \tilde{\sigma} \end{bmatrix}$ and the associated deviatoric stress tensor is

$$g\mathbf{S} = \begin{bmatrix} -\frac{1}{3}\tilde{\sigma} & 0 & 0 \\ 0 & -\frac{1}{3}\tilde{\sigma} & 0 \\ 0 & 0 & \frac{2}{3}\tilde{\sigma} \end{bmatrix}$$

$$- \mathrm{Tr}(\mathbf{M}^g \mathbf{S}) = \frac{2\tilde{\sigma}}{3}$$

$$-(\mathbf{M}^g \mathbf{S} + g\mathbf{SM})^D = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \frac{4}{3}\tilde{\sigma} \end{bmatrix}^D = \begin{bmatrix} -\frac{4}{9}\tilde{\sigma} & 0 & 0 \\ 0 & -\frac{4}{9}\tilde{\sigma} & 0 \\ 0 & 0 & \frac{8}{9}\tilde{\sigma} \end{bmatrix}$$

- using the constitutive equation, we obtain for $gD_{33}$:

$$gD_{33} = \delta_1 \frac{4}{3}\tilde{\sigma} + \delta_2 \frac{8}{9}\tilde{\sigma} + \delta_3 \frac{8}{9}\tilde{\sigma}$$

* Then a uniaxial stress $\tilde{\sigma}$ is applied normal to the c-axis, for instance in the direction 1. The stress tensor is

$$g\sigma = \begin{bmatrix} \tilde{\sigma} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}$$

and the associated diatoric stress tensor is

$$g\mathbf{S} = \begin{bmatrix} \frac{2}{3}\tilde{\sigma} & 0 & 0 \\ 0 & -\frac{1}{3}\tilde{\sigma} & 0 \\ 0 & 0 & -\frac{1}{3}\tilde{\sigma} \end{bmatrix}$$

$$- \mathrm{Tr}(\mathbf{M}^g \mathbf{S}) = -\frac{\tilde{\sigma}}{3}$$

$$-(\mathbf{M}^g \mathbf{S} + g\mathbf{SM})^D = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -\frac{2}{3}\tilde{\sigma} \end{bmatrix}^D = \begin{bmatrix} \frac{2}{9}\tilde{\sigma} & 0 & 0 \\ 0 & \frac{2}{9}\tilde{\sigma} & 0 \\ 0 & 0 & \frac{4}{9}\tilde{\sigma} \end{bmatrix}$$

- using the constitutive equation, we obtain for $gD_{11}$:

$$gD_{11} = \delta_1 \frac{4}{3}\tilde{\sigma} + \delta_2 \frac{2}{9}\tilde{\sigma} + \delta_3 \frac{2}{9}\tilde{\sigma}$$
* Finally from the grain behavior, we can write that $\gamma^g D_{33} =^g D_{11}$:

$$ \gamma \left( \frac{4}{3}\delta_1 + \frac{8}{9}\delta_2 + \frac{8}{9}\delta_3 \right) = \frac{4}{3}\delta_1 + \frac{2}{9}\delta_2 + \frac{2}{9}\delta_3 $$

* This gives:

$$ \delta_2 = 6\delta_1 \frac{1-\gamma}{4\gamma-1} - \delta_3 $$

* and using $\delta_1$ and $\delta_3$ as defined above:

$$ \delta_2 = \frac{\psi_1}{2} \left( \beta \frac{\gamma + 2}{4\gamma - 1} - 1 \right) $$

## Appendix B. Strain and stress field comparisons

The Figures 1 and 2 show the probability distribution functions for all components of the strain and stress tensors during uniaxial creep with c-axis evolution (Section 4.1).

<figure><img src="image_16.png" /><figcaption>**Figure 1.** Kernel Density Estimation of all components of the strain tensor predicted by an elasto-viscoplastic (EVP) simulation using CraFT and the continuous transverse isotropic (CTI) using R³iCe after 0.01 macroscopic strain. The dashed lines show the results for simulations where grains are well discretized. The full line shows the result for the simulation with one orientation per element.</figcaption></figure>
* Finally from the grain behavior, we can write that $\gamma^g D_{33} =^g D_{11}$:

$$ \gamma \left( \frac{4}{3} \delta_1 + \frac{8}{9} \delta_2 + \frac{8}{9} \delta_3 \right) = \frac{4}{3} \delta_1 + \frac{2}{9} \delta_2 + \frac{2}{9} \delta_3 $$

* This gives:

$$ \delta_2 = 6\delta_1 \frac{1-\gamma}{4\gamma-1} - \delta_3 $$

* and using $\delta_1$ and $\delta_3$ as defined above:

$$ \delta_2 = \frac{\psi_1}{2} \left( \beta \frac{\gamma + 2}{4\gamma - 1} - 1 \right) $$

## Appendix B. Strain and stress field comparisons

The Figures 1 and 2 show the probability distribution functions for all components of the strain and stress tensors during uniaxial creep with c-axis evolution (Section 4.1).

<figure><img src="image_17.png" /><figcaption>**Figure 1.** Kernel Density Estimation of all components of the strain tensor predicted by an elasto-viscoplastic (EVP) simulation using CraFT and the continuous transverse isotropic (CTI) using R³iCe after 0.01 macroscopic strain. The dashed lines show the results for simulations where grains are well discretized. The full line shows the result for the simulation with one orientation per element.</figcaption></figure>