%~Mouliné par MaN_auto v.0.32.0 2024-03-25 12:45:46
\documentclass[CRMECA, Unicode, XML, screen,biblatex,longabstracts]{cedram}
%


\addbibresource{crmeca20230896.bib}

\usepackage{algpseudocode}
%\usepackage{algorithmicx}
\usepackage{algorithm}
%\usepackage{algcompatible}
\usepackage{circuitikz}
\usepackage[capitalise,noabbrev]{cleveref}
\usepackage{mathtools}
\usepackage{amsbsy}

%\usepackage{marvosym}
%\usepackage{tikzsymbols}

\usepackage{float}
\usepackage{xspace}
%\usepackage{pstricks}
%\usepackage{auto-pst-pdf}
%\usepackage[inkscapelatex=true]{svg}


\newcommand{\noopsort}[1]{}


\DeclareRobustCommand{\PhaseOne}{[\tikz[baseline=-\the\dimexpr\fontdimen22\textfont2\relax,inner sep=0pt,thick]\draw[red](0,0)--(.5,0);]}
\DeclareRobustCommand{\PhaseTwo}{[\tikz[baseline=-\the\dimexpr\fontdimen22\textfont2\relax,inner sep=0pt,thick]\draw[blue](0,0)--(.5,0);]}


\newcommand{\grx}{\ensuremath{\Gamma_{RX}}\xspace}
\newcommand{\rrrice}{R$^{3}$iCe\xspace}
\newcommand{\ctirx}{CTI-RX\xspace}
\newcommand{\degree}{\ensuremath{^\circ}}


\newcommand{\bfa}{\mathbf{a}}
\newcommand{\bfb}{\mathbf{b}}
\newcommand{\bfc}{\mathbf{c}}
\newcommand{\bfy}{\mathbf{y}}
\newcommand{\bfu}{\mathbf{u}}
\newcommand{\bfx}{\mathbf{x}}
\newcommand{\bfA}{\mathbf{A}}
\newcommand{\bfB}{\mathbf{B}}
\newcommand{\bfD}{\mathbf{D}}
\newcommand{\bfM}{\mathbf{M}}
\newcommand{\bfS}{\mathbf{S}}
\newcommand{\bfW}{\mathbf{W}}
\newcommand{\bfX}{\mathbf{X}}
\newcommand{\clB}{\mathcal{B}}
\newcommand{\fMD}{\mathbf{MD}}
\newcommand{\fMS}{\mathbf{MS}}
\newcommand{\fSM}{\mathbf{SM}}
\newcommand{\fDM}{\mathbf{DM}}
\newcommand{\fWc}{\mathbf{Wc}}
\newcommand{\fDc}{\mathbf{Dc}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Colored tables and algorithms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\RequirePackage{adjustbox}
\RequirePackage{float}
\newfloat{Algorithm}{ht}{lop}
\RequirePackage{algpseudocode}
\RequirePackage{multicol}

\renewcommand{\arraystretch}{1.25}
\BeforeBeginEnvironment{tabular}{%
\centering\small\begin{adjustbox}{bgcolor=black!6}%
}%
\AfterEndEnvironment{tabular}{%
\end{adjustbox}%
}%
\BeforeBeginEnvironment{tabular*}{%
\centering\small\begin{adjustbox}{bgcolor=black!6}%
}%
\AfterEndEnvironment{tabular*}{%
\end{adjustbox}%
}%
\BeforeBeginEnvironment{algorithmic}{%
\centering\small%
\begin{adjustbox}{cfbox=black!6 .5pt 5pt,bgcolor=black!6}%%
\begin{minipage}[t]{\dimexpr\linewidth-11pt}%
}%
\AfterEndEnvironment{algorithmic}{%
\end{minipage}%
\end{adjustbox}%
}%compiler

%\RequirePackage[nameinlink,noabbrev]{cleveref}
\crefname{Algorithm}{Algorithm}{Algorithms}

%######## END ###########


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

%\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\TopicEN{Section information missing}
%\TopicFR{Information section manquante}

\title{A physically-based formulation for texture evolution during dynamic recrystallization. A~case study of ice}

\alttitle{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}

\author{\firstname{Thomas} \lastname{Chauve}\CDRorcid{0000-0003-2916-7472}\IsCorresp}
\address{Univ. Grenoble Alpes, CNRS, IRD, G-INP, IGE, Grenoble, France.}
\email[T. Chauve]{thomas.chauve@univ-grenoble-alpes.fr}

\thanks{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).}

\author{\firstname{Maurine} \lastname{Montagnat}\CDRorcid{0000-0001-9436-5163}}
\address[1]{Univ. Grenoble Alpes, CNRS, IRD, G-INP, IGE, Grenoble, France.}
\address{Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d’Études de la Neige, Grenoble, France.}
\email[M. Montagnat]{maurine.montagnat@univ-grenoble-alpes.fr}


\author{\firstname{Véronique} \lastname{Dansereau}\CDRorcid{0000-0002-6766-4595}}
\address{Institut des Sciences de la Terre, Université Grenoble Alpes, CNRS (5275), Grenoble, France.}
\email[V. Dansereau]{veronique.dansereau@univ-grenoble-alpes.fr}

\author{\firstname{Pierre} \lastname{Saramito}\CDRorcid{0000-0003-1989-2566}}
\address{Lab. Jean Kuntzmann, CNRS, Université Grenoble-Alpes, F-38041 Grenoble, France.}
\email[P. Saramito]{pierre.saramito@imag.fr}

\author{\firstname{Kévin} \lastname{Fourteau}\CDRorcid{0000-0002-9905-2446}}
\address[2]{Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d’Études de la Neige, Grenoble, France}
\email[K. Fourteau]{kevin.fourteau@meteo.fr}

\author{\firstname{Andréa} \lastname{Tommasi}\CDRorcid{0000-0002-6457-1852}}
\address{Géosciences Montpellier - CNRS, Université de Montpellier, France.}
\email[A. Tommasi]{andrea.tommasi@umontpellier.fr}
\shortrunauthors

\keywords{dynamic recrystallization, texture, viscoplastic anisotropy, finite-element method, ice}

\altkeywords{recristallisation dynamique, texture, anisotropie viscoplastique, méthode des éléments finis, glace}

\begin{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 \rrrice model is solved by a finite-element method with a semi-implicit scheme coded using the Rheolef library. The resulting open-source software
\rrrice 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 present formulation is too time-consuming for direct implementation in large-scale ice flow models, \rrrice can be used to adjust the parameterization used to implement texture-induced anisotropy in these models. The validation was performed for ice, but the \rrrice implementation is generic and applies to any material whose behavior may be adequately described using an anisotropic flow law.
\end{abstract}


\begin{altabstract}
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 attracteur d'orientation qui maximise la contrainte de cisaillement résolue sur le système de glissement le plus facile dans le cristal. Elle est mise en {\oe}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{\oe}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 {\oe}uvre directe dans des modèles d'écoulement de glace à grande échelle, R3iCe peut être utilisé pour ajuster la paramétrisation utilisée pour mettre en {\oe}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.
\end{altabstract}

\begin{DefTralics}
\newcommand{\rrrice}{R$^{3}$iCe\xspace}
\end{DefTralics}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
%\vspace*{-20pt}
\maketitle
\vspace*{-15pt}

\section*{Notation}
{\small
\begin{table}[!h]
\caption{Mathematical notations}
\vspace*{-10pt}
\centering
\begin{tabular}{llll}
$a$ & Scalar \\
$\bfa$ & Vector \\
$\bfA$ & Tensor \\
$\bfA^D$ & Deviatoric part of $\bfA$ & $\bfA^D=\bfA-\frac{1}{dim(\bfA)}Tr(\bfA)$ & $A^D_{ij}=A_{ij}-\frac{\delta_{ij}}{dim(\bfA)} Tr(\bfA)$\\
$\otimes$ & Tensorial (dyadic) product & $\bfA=\bfb\otimes \bfc$ & $A_{ij}=b_ic_j$ \\
$Tr()$ & Trace of a Tensor & $Tr(\bfA)$ & $A_{ii}$ \\
$dim()$ & dimension of a tensor & $dim(\bfA)$ & \\
$\nabla$ & $\nabla$ operator & $\left(\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z}\right)$ & $\frac{\partial}{\partial x_i}$ \\
$.$ & Scalar product & $\bfa.\bfb$ & $a_ib_i$ \\
$:$ & Double contracted product & $\bfA:\bfB$ & $A_{ij}B_{ij}$\\
$\|. \|$ & Norm & $\| \bfa \|$ & $a_i^2$ \\
$\bfA^n$ & matrix product $n$ times & $\underbrace{\bfA\times \bfA\dots \times \bfA}_n$
\end{tabular}
\label{tab:notation}
\end{table}
}

%\break
\section{Introduction}
Dynamic recrystallization (DRX) is a mechanism by which local stress incompatibilities and local strain heterogeneities are relaxed~\cite{Humphreys2004b, Sakai2014, Chauve2015}. By \emph{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) \cite{Kamb1972, Poirier1985, Humphreys1996, Doherty1997, Sakai2014}. 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 (\cref{fig:CTI_schema}) 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~\cite{Glen1954, Duval1983}. 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~\cite{Castelnau1998}.

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~\cite{De-la-chapelle1998}. 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 ``discontinuous'' dynamic recrystallization regime). Whether continuous or discontinuous DRX predominates depends on the temperature and strain conditions along the ice core~\cite{Thorsteinsson1997, Montagnat2012}. 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.~\cite{Samyn2008, Smith2017_E, Minchew2018, Thomas2021}). 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~\cite{Pettit2007, Gillet2006, Martin2009, Ma2010, Budd2013, McCormack2022}. 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.~\cite{Duval1979, Duval1981, Bouchez1982, Jacka1984b, Jacka1994, Montagnat2015, Chauve2015, Qi2017, Journaux2019}). Kamb (1972) \cite{Kamb1972} 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~\cite{Grennerat2012} 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~\cite{Chauve2017} 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.

\break
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~\cite{Duval1981, Jacka1984b, Budd2013}.
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~\cite{jacka1991, Jun1998, Li1996, Treverrow2012}.
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) \cite{Glen1952, Duval1983}. 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~\cite{Budd2013, Treverrow2012} 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 $\bfc$-axes oriented within a girdle at an angle between 30$^{\circ}$ and 50$^{\circ}$ 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\%$~\cite{Kamb1972, Duval1979, Jacka1984b, Treverrow2012, Piazolo2013, Montagnat2015, Qi2017}.
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 $\bfc$-axes aligned in a direction perpendicular to the shear plane~\cite{Hudleston1977, Duval1981, Bouchez1982, Journaux2019, Qi2019}.

Texture development in ice deforming by dislocation creep has been successfully simulated using approaches of variable complexity~\cite{Alley1988, Pettit2007, Gillet2006, Castelnau1996c}. 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~\cite{Wenk1997, Thorsteinsson2001, Llorens2016, Richards2021, Rathmann2021}. 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~\cite{Wenk1997, Thorsteinsson2001}, or from limitations in the deformation modeling frame despite the high complexity of microstructure evolution~\cite{Llorens2016}. In~\cite{Richards2021,Rathmann2021}, 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~\cite{Castelnau1996b}.

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 \ctirx, is solved in a full-field model using the finite-element method (FEM).

In \cref{sec:num_model}, we present the \ctirx model: its equations for the CTI law (\cref{subsec:CTI_law}) and the proposed formulation for the evolution of the $\bfc$-axis (\cref{subsec:c_evolution}). The numerical scheme designed to solve these equations is described in \cref{sec:r3ice_implementation}. \Cref{sec:validation_r3ice} presents the calibration and validation for cases without dynamic recrystallization. \Cref{sec:resultat} 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.


\section{\ctirx- a model for the ice single crystal behavior}\label{sec:num_model}
\subsection{ Continuous transverse isotropic behavior }\label{subsec:CTI_law}

Following~\cite{Meyssonnier1996,Mansuy2002,Gillet2005,Burr2017,Rathmann2021b}, 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 $\bfc$-axis, being isotopic in the basal plane (see Figure~\ref{fig:CTI_schema}). This formulation is based on theoretical work by Boehler (1987) \cite{Boehler1987}, which has shown, through geometrical arguments, that the relation between the deviatoric stress $\bfS$ and the deviatoric strain rate $\bfD$ in a transverse isotropic symmetry may be described by a potential $\phi_S$:
\begin{equation}
\bfD=\left(\frac{\partial \phi_S}{\partial\bfS}\right)^D
\end{equation}
where $\phi_S$ is expressed on the $\clB_{CTI}$ base, which is built from a combination of the following four invariants:
\begin{equation}\label{eq:CTI_invariant}
\clB_{CTI}=\left[Tr\left(\bfS^2\right),Tr\left(\bfS^3\right),Tr(\bfM\bfS),
Tr\left(\bfM\bfS^2\right)\right],
\end{equation}
With $\bfM =\mathbf{e_3} \otimes \mathbf{e_3}$ where $\mathbf{e_3}$ is the direction of the plane of isotropy (\cref{fig:CTI_schema}). The potential $\phi_S$ is described below in the linear ($\phi_S^{(1)}$) and nonlinear ($\phi_S^{(n)}$) cases.

\begin{figure}[!ht]
\begin{minipage}[t]{0.47\textwidth}\centering
\includegraphics[width=0.6\columnwidth]{figures/single_crsytal_sketch.png}
\end{minipage}
\qquad
\begin{minipage}[t]{0.47\textwidth}\centering
\resizebox{0.6\columnwidth}{!}{
\begin{circuitikz}
\tikzstyle{every node}=[font=\LARGE]
\draw [fill={rgb,255:red,31; green,0; blue,255}, line width=2pt] (7.75,22) ellipse (4cm and 1cm);
\draw [line width=4.7pt, ->, >=Stealth] (7.75,22).. controls (7.75,26) and (7.75,25.5).. (7.75,29);
\node [font=\Huge] at (8.5,29) {$e_3$};
\end{circuitikz}}
\end{minipage}
\caption{\label{fig:CTI_schema}
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 $\bfc$-axis (i.e. the normal to the basal $(0001)$ plane).}
\end{figure}

\subsubsection{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 $\clB_{CTI}$ base as:
\begin{equation}
\phi_S^{(1)}=\delta_1 Tr\left(\bfS^2\right)+\delta_2 Tr(\bfM\bfS)^2+\delta_3 Tr\left(\bfM\bfS^2\right).
\end{equation}
Considering the following derivation rules,
\begin{equation}\label{eq:derivation_rule}
\frac{\partial Tr(\fMS)}{\partial \bfS}=\bfM,\quad \frac{\partial Tr\left(\bfS^2\right)}{\partial \bfS}=2\bfS,\quad \frac{\partial Tr\left(\fMS^2\right)}{\partial \bfS}=\fMS+\fSM,
\end{equation}
the constitutive relation is given by
\begin{equation}\label{eq:constitutive_eq}
\bfD=\left(\frac{\partial \phi_S^{(1)}}{\partial \bfS}\right)^D=2\delta_1 \bfS + 2\delta_2 \bfM^D Tr(\fMS) + \delta_3 (\fMS+\fSM)^D.
\end{equation}
The three independent parameters $\delta_i$ in \cref{eq:constitutive_eq} 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 ${}^g\bfX$ indicates a tensor expressed in the grain reference frame:
\begin{itemize}
\item $\psi_1$ is the fluidity of the ice crystal subjected to shear on the basal plane:
\begin{equation}
{}^{g}D_{13}=\frac{\psi_1}{2}{}^{g}S_{13},\quad{}^{g}D_{23}=\frac{\psi_1}{2}{}^{g}S_{23},
\end{equation}
\item $\beta$ is the ratio between the fluidity for shear parallel to the basal plane and that for shear normal to the basal plane:
\begin{equation}
^{g}D_{12}=\beta\frac{\psi_1}{2}{}^{g}S_{12},
\end{equation}
\item $\gamma$ is the ratio of the viscosity for compression (or traction) Along the $\bfc$-axis $\eta_c$ and for compression (or traction) in any direction ${}^ge_r$ within the basal plane:
\begin{equation}
{}^{g}S_{33}=2\eta_c{}^{g}D_{33}, \quad \gamma {}^{g}S_{rr}=2\eta_c{}^{g}D_{rr}
\end{equation}
For solicitations of similar amplitude in different directions, for instance along $\mathbf{e_3}$ or perpendicular to $\mathbf{e_3}$ (i.e. along $\mathbf{e_{r}}$) $S={}^gS_{33}={}^gS_{rr}$, in the linear case, the strain rates are related by:
\begin{equation}
{}^gD_{rr}=\gamma{}^gD_{33}
\end{equation}
\end{itemize}
Using the above constraints on single crystal behavior, identification gives (see \cref{sec:demo_di} for details):
\begin{equation}
\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}\left(1-\beta\right)
\end{equation}


The constitutive equation~\eqref{eq:constitutive_eq} can also be expressed in an inverse manner as follows:
\begin{equation}
\bfS=2\alpha_1 \bfD + 2\alpha_2 \bfM^D Tr(\fMD) + \alpha_3 (\fMD+\fDM)^D
\end{equation}
and in terms of the potential
\begin{equation}
\phi_D^{(1)}=\alpha_1 Tr\left(\bfD^2\right)+\alpha_2 Tr(\bfM\bfD)^2+\alpha_3 Tr\left(\bfM\bfD^2\right)
\end{equation}
with : $\alpha_1= \frac{\eta_1}{\beta},\quad \alpha_2=2\eta_1(\frac{\gamma}{\beta}-1),\quad\alpha_3=2\eta_1(1-\frac{1}{\beta})$ and $\eta_1=\frac{1}{\psi_1}$ (see~\cite{Gillet2006} and appendix A for details).


\subsubsection{The non-linear CTI formulation}

A non-linear potential $\phi_D^{(n)}$ linking $\bfS$ to $\bfD^{\frac{1}{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$ :
\begin{equation}
\phi_D^{(n)}= \frac{1}{k} \left(\tilde{\phi}_D^{(1)}\right)^k.
\end{equation}
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~\cite{Gagliardini1999,Gillet2006,Rathmann2021}. Using the chain rule for deriving $\phi_D^{n}$ gives:
\begin{equation}\label{eq:non_linear_formulation}
\bfS=\frac{\partial \phi_D^{(n)}}{\partial \bfD}= \left(\tilde{\phi}_D^{(1)}\right)^{k-1} \frac{\partial \tilde{\phi}_D^{(1)}}{\partial \bfD}
\end{equation}
$k$ is further obtained by identification of the power exponent on \cref{eq:non_linear_formulation} in order to have $\bfS\sim \bfD^{\frac{1}{n}}$:
\begin{equation}
2(k-1)+2=\frac{1}{n} \implies k=\frac{n+1}{2n}
\end{equation}

Therefore, a non-linear CTI law can be written as:
\begin{equation}\label{eq:final_CTI_eq}
\phi_D^{(n)}= \frac{4n}{n+1} \left(\alpha_1 tr\left(\bfD^2\right) + \alpha_2 tr(\fMD)^2+\alpha_3 tr\left(\fMD^2\right) \right)^{\frac{n+1}{2n}},
\end{equation}
Leading to the constitutive equation:
\begin{equation}\label{eq:CTI_law}
\bfS=\eta_n^\star \left(2\alpha_1 \bfD + 2\alpha_2 \bfM^D Tr(\fMD) + \alpha_3 (\fMD+\fDM)^D \right),
\end{equation}
with effective viscosity
\begin{equation}\label{eq:eta}
\eta_n^\star = 2\eta_n\left(\alpha_1 tr\left(\bfD^2\right) + \alpha_2 tr(\fMD)^2+\alpha_3 tr\left(\fMD^2\right) \right)^{\frac{1-n}{2n}}
\end{equation}
and $\alpha_i$ parameters:
\begin{equation}\label{eq:alpha_i}
\alpha_1= \frac{1}{2\beta},\quad \alpha_2=\left(\frac{\gamma}{\beta}-1\right),\quad \alpha_3=\left(1-\frac{1}{\beta}\right).
\end{equation}
It is worth noting that $\beta$ links ${}^gD_{12}$ to ${}^gS_{12}$ through
\begin{equation}\label{eq:SC_shear_12}
{}^gD_{12}=\frac{\psi_n}{2}\beta^{\frac{n+1}{2}}{}^gS_{12}^n.
\end{equation}
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}}$.


\subsection{ Evolution of crystal rotation in response to deformation and recrystallization }\label{subsec:c_evolution}

Deformation resulting from dislocation glide in the different slip plane results in a rotation of the crystal whose rate is described as in~\cite{Gillet2006b,Gagliardini2009} by
\begin{equation}\label{eq:c_dot_old}
\frac{\partial \bfc}{\partial t}=\fWc-\lambda\left[\fDc-\left(\bfc^T\fDc\right)\bfc\right],
\end{equation}
where $\bfW(\bfu)=\frac{\mathbf{\nabla}(\bfu)-\mathbf{\nabla}^T(\bfu)}{2}$ and $\bfD(\bfu)=\frac{\mathbf{\nabla}(\bfu)+\mathbf{\nabla}^T(\bfu)}{2}$ are the spin- and strain-rate tensors. According to \cref{eq:c_dot_old}, the rotation of the $\bfc$-axis is due to the bulk spin, $\fWc$, and the viscoplastic spin expressed by $\lambda[\fDc-(\bfc^T\fDc)\bfc]$. $\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~\cite{Gillet2006b, Gagliardini2009}. This assumption is well funded for ice and has often been made in previous studies~\cite{Lliboutry1993,Azuma1995,Svendsen1996,Godert1998,Thorsteinsson2001}.

To account for the impact of dynamic recrystallization on the evolution of the $\bfc$-axis orientation, we add a new term to this equation based on the assumption that dynamic recrystallization induces a local re-orientation of the $\bfc$-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 $\bfc$-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^{\circ}$~\cite{Chauve2017}.


With the addition of the dynamic recrystallization term, the $\bfc$-axis rotation is described as follows:
\begin{equation}\label{eq:c_dot_new}
\frac{\partial \bfc}{\partial t}=\fWc-\lambda\left[\mathbf{Dc-\left(\bfc^T\fDc\right)\bfc}\right]+\frac{1}{\grx}\left(\mathbf{c_0}-\bfc\right),
\end{equation}
where \grx is a parameter that controls the rate of rotation of the $\bfc$-axis toward an attractor, $\mathbf{c_0}$, which represents the ``ideal'' orientation produced by recrystallization and whose formulation is described in \cref{subsec:c0_attractor}. Therefore, the evolution of the $\bfc$-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 \grx can be illustrated by considering the artificial case in which only dynamic recrystallization contributes to the $\bfc$-axis rotation. In this case, \cref{eq:c_dot_new} becomes
\begin{equation}
\frac{\partial \bfc}{\partial t}=\frac{1}{\grx}\left(\mathbf{c_0}-\bfc\right),
\end{equation}
The norm of $\frac{\partial\bfc}{\partial t}$ can be interpreted as the rotation rate due only to dynamic recrystallization. This rate increases as $\|\mathbf{c_0}-\bfc\|$. Considering a constant $\mathbf{c_0}$, the solution for the $\bfc$-axis rotation is:
\begin{equation}
\bfc(t)=\left(\bfc(0)-\mathbf{c_0} \right) e^{-\frac{t}{\grx}}+\bfc(0).
\end{equation}
In this case, $\bfc$ rotates towards $\mathbf{c_0}$ at a rate following an exponential decay with a characteristic time \grx, such that
$95\%$ ($99\%$) of the rotation is obtained after $4\grx$ ($5\grx$). In the \ctirx model, $\mathbf{c_0}$ is not constant but depends on the local diatoric stress tensor, $\bfS$.


\subsection{Formulation of the recrystallization attractor \texorpdfstring{$\mathbf{c_0}$}{c0}}\label{subsec:c0_attractor}

The formulation of $\mathbf{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~\cite{Chauve2015, Chauve2017, Chauve2017b}. 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 $\mathbf{a_i}$ (\cref{fig:CTI_schema}) gliding on the basal plane $\bfc$ in response to a diatoric stress $\bfS$ is given by
\begin{equation}
RSS(\mathbf{a_i},\bfc,\bfS)= \bfS:\boldsymbol{\mu},
\end{equation}
with $\boldsymbol{\mu}=\frac{1}{2} \frac{\bfc\otimes\mathbf{a_i}+\mathbf{a_i}\otimes\bfc}{\| \bfc \|.\| \mathbf{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 $\mathbf{v_i}$ ($i = 1-3$) as the corresponding eigenvectors of the diatoric stress tensor $\bfS$, the maximum resolved shear stress is $\| s_1-s_3\|$ and the associated orientations are~\cite{Wu2004}:
\begin{equation}\label{eq:c0_solution}
\mathbf{c_0}=\frac{1}{2}\left(\mathbf{v_1} \pm \mathbf{v_3}\right)
\end{equation}

Degenerated cases exist if at least two eigenvalues are equal:
\begin{itemize}
\item $s_1=s_2$, the possible attractors $\bfc_0$ are the orientations at $45^\circ$ from $\mathbf{v_3}$.
\item $s_2=s_3$, the possible attractors $\bfc_0$ are the orientations at $45^\circ$ form $\mathbf{v_1}$.
\item $s_1=s_2=s_3=0$, the possible attractors $\bfc_0$ are all orientations.
\end{itemize}
In all other cases, there are two possible orientations for the attractor $\mathbf{c_0}$. To solve \cref{eq:c_dot_new} we define $\mathbf{c_0}$ as which of the two solutions is closer to the actual $\bfc$ axis, in order to minimize the path between $\bfc$ and $\mathbf{c_0}$. In practice, $\mathbf{c_0}$ is chosen so that the scalar product between $\mathbf{c_0}$ and $\bfc$ is maximal.

\subsection{Interpretation of the orientation evolution equation}

\Cref{eq:c_dot_new} can be interpreted in terms of the relative impact of deformation versus recrystallization on the orientation of the $\bfc$-axis. In the following section,we describe three end-member cases.

For uniaxial compression along the $\bfy$-axis (Figure~\ref{fig:cdot_intper}$\MK$(a)), the stable solution for the deformation term $(\fWc-\lambda[\mathbf{Dc-(\bfc^T\fDc
)\bfc}])$ corresponds to $\bfc$ aligned with $\bfy$ (red dot, Figure~\ref{fig:cdot_intper}$\MK$(a)) \cite{Gillet2006}. For recrystallization $(\frac{1}{\grx}(\mathbf{c_0}-\bfc))$, the stable solutions are the orientations at $45^\circ$ to the $\bfy$-axis (blue line, Figure~\ref{fig:cdot_intper}$\MK$(a)). When both deformation and recrystallization are active, the stable orientations of the $\bfc$-axis are located on a girdle between these two end-members (black line, Figure~\ref{fig:cdot_intper}$\MK$(a)). The position of this girdle depends on the value of \grx (black arrows, Figure~\ref{fig:cdot_intper}$\MK$(a)). The higher is \grx the closer is the girdle to the compression axis.

For uniaxial tension along $\bfy$, the stable solutions for the deformation term are all the orientations located in the $xOz$ plane (red line, Figure~\ref{fig:cdot_intper}$\MK$(b)) \cite{Gillet2006}. When both deformation and recrystallization are active, the stable orientations (black circle, Figure~\ref{fig:cdot_intper}$\MK$(b)) form a girdle between $45^\circ$ and $90^\circ$ to the $\bfy$-axis depending on the value of \grx.

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

\begin{figure}[!ht]
\centering
\includegraphics[width=14cm]{figures/cdot_interp.png}
\caption{\label{fig:cdot_intper}
Pole figure schematic representation of the stable orientations for the evolution equation (eq.~\eqref{eq:c_dot_new}) for various stress configurations $\bfS$: (a) uniaxial compression along $\bfy$, (b) uniaxial tension along $\bfy$ and (c) simple shear in the $xOz$ plane, $\frac{\partial u_x}{ \partial y}<0$ (upper hemisphere).}
\end{figure}


\subsection{\ctirx flow model description}\label{sec:r3ice_model_ad}

The system of equations constituting the \ctirx 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
\begin{enumerate}
\item The momentum equation under the Stokes assumption:
\begin{equation}
-\nabla.\bfS+\nabla p=0
\end{equation}
\item The non-linear CTI law:
\begin{equation}
\begin{aligned}[t]
\bfS=\eta_n^\star \left(2\alpha_1 \bfD + 2\alpha_2 \bfM^D Tr(\fMD) + \alpha_3 (\fMD+\fDM)^D \right) \\
\eta_n^\star = 2\eta_n\left(\alpha_1 tr(\bfD^2) + \alpha_2 tr(\fMD)^2+\alpha_3 tr\left(\fMD^2\right) \right)^{\frac{1-n}{2n}}
\end{aligned}
\end{equation}
\item The equation for the evolution of the $\bfc$-axis:
\begin{equation}
\frac{\partial \bfc}{\partial t}=\fWc-\lambda\left[\mathbf{Dc-\left(\bfc^T\fDc\right)\bfc}\right]+\frac{1}{\grx}\left(\mathbf{c_0}-\bfc\right),
\end{equation}
With $\mathbf{c_0}=\frac{v_1 \pm v_3}{2}$ where $\mathbf{v_i}$ are the eigenvectors of $\bfS$ and $\|\bfc.\mathbf{c_0}\|$ is maximum.
\end{enumerate}

\section{Numerical scheme and implementation}\label{sec:r3ice_implementation}

The \ctirx equations (\cref{sec:r3ice_model_ad}) are solved on a cuboid of polycrystalline ice deformed by creep (stress boundary conditions) under either uniaxial compression, uniaxial tension, or simple shear (Figure~\ref{fig:BC_3D}). A full-field approach using FEM was chosen to describe the field of orientations $\bfc$ at the element size. This implementation results in the \rrrice model (\textbf{R}heology, \textbf{R}ecrystallization, \textbf{R}heolef, \textbf{i}n \textbf{C}ontinous Transverse Isotropic mat\textbf{e}rial). This design allows direct comparisons of texture evolution and macroscopic mechanical behavior with data from laboratory creep experiments~\cite{Montagnat2015,Journaux2019,Qi2017}.

\begin{figure}[!ht]
\centering
\resizebox{0.4\textwidth}{!}{
\begin{circuitikz}
\tikzstyle{every node}=[font=\normalsize]
\draw [, line width=1pt] (17,21.75) rectangle (21,25.75);
\draw [, line width=1pt, dashed] (18.5,23) rectangle (22.25,27);
\draw [line width=1pt] (17,25.75) -- (18.5,27);
\draw [line width=1pt] (18.5,27) -- (22.25,27);
\draw [line width=1pt] (22.25,27) -- (22.25,23);
\draw [line width=1pt] (21,21.75) -- (22.25,23);
\draw [line width=1pt] (21,25.75) -- (22.25,27);
\draw [line width=1pt, dashed] (17,21.75) -- (18.5,23);
\node [font=\Large] at (19.30,22.35) {\emph{bottom}};
\node [font=\Large] at (19.30,26.30) {\emph{top}};
\node [font=\Large, color={rgb,255:red,0; green,0; blue,255}] at (18.5,23.45) {left\_bottom\_back};
\node [font=\Large, color={rgb,255:red,0; green,0; blue,255}] at (17,21.2) {left\_bottom\_front};
\draw [color={rgb,255:red,136; green,138; blue,133}, line width=1pt, -Stealth] (22.25,21.75) -- (22.25,22.5);
\draw [color={rgb,255:red,136; green,138; blue,133}, line width=1pt, -Stealth] (22.25,21.75) -- (23,21.75);
\draw [color={rgb,255:red,136; green,138; blue,133}, line width=1pt, -Stealth] (22.25,21.75) -- (22.75,22.25);
\node [font=\normalsize, color={rgb,255:red,136; green,138; blue,133}] at (23,21.5) {x};
\node [font=\normalsize, color={rgb,255:red,136; green,138; blue,133}] at (22,22.25) {y};
\node [font=\normalsize, color={rgb,255:red,136; green,138; blue,133}] at (23,22.25) {z};
\draw [color={rgb,255:red,0; green,0; blue,255}, fill={rgb,255:red,0; green,0; blue,255}, line width=1pt] (17,21.75) circle (0.25cm);
\draw [color={rgb,255:red,0; green,0; blue,255}, fill={rgb,255:red,0; green,0; blue,255}, line width=1pt] (18.5,23) circle (0.25cm);
\draw [color={rgb,255:red,0; green,0; blue,255}, fill={rgb,255:red,0; green,0; blue,255}, line width=1pt] (18.5,23) circle (0.25cm);
\draw [color={rgb,255:red,0; green,0; blue,255}, fill={rgb,255:red,0; green,0; blue,255}, line width=1pt] (18.5,23) circle (0.25cm);
\draw [color={rgb,255:red,0; green,0; blue,255}, fill={rgb,255:red,0; green,0; blue,255}, line width=1pt] (18.5,23) circle (0.25cm);
\draw [color={rgb,255:red,0; green,0; blue,255}, fill={rgb,255:red,0; green,0; blue,255}, line width=1pt] (18.5,23) circle (0.25cm);
\draw [color={rgb,255:red,0; green,0; blue,255}, fill={rgb,255:red,0; green,0; blue,255}, line width=1pt] (18.5,23) circle (0.25cm);
\draw [Stealth-Stealth] (22.5,23) -- (22.5,27);
\draw [line width=1.1pt, -Stealth] (20,26.25) -- (20,28);
\draw [line width=1.1pt, -Stealth] (20,26.25) -- (21.5,26.25);
\node [font=\Large] at (20.2,27.5) {$F_y$};
\node [font=\Large] at (21.25,26.60) {$F_x$};
\node [font=\Large] at (23.,24.8) {$L$};
\end{circuitikz}}\label{fig:BC_3D}
\caption{Schematic representation of the rectangular cubic simulation domain, with surface and point boundaries labeled.}
\end{figure}

\subsection{Numerical experiment setup}
\subsubsection{Uniaxial creep tests}

In uniaxial creep, a zero vertical velocity ($u_y(\textit{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 ($\bfu(\textit{left\_bottom\_back})=(0,0,0)$). To avoid rotation around the $\bfy$-axis, a zero-velocity condition along the $\bfx$-axis is also applied on a second summit ($u_x(\textit{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(\textit{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 \cref{sec:EVPvsCTI}.


\subsubsection{Simple shear creep tests}

In the simple shear tests, forcing is applied as a constant tangential stress on the top boundary ($F_x(\textit{top})=F$) of the sample. Other boundary conditions are:
\begin{itemize}
\item $\bfu(\textit{bottom})=0$,
\item $u_y=u_z=0$ for the lateral surfaces (left, right, front, back) to avoid solid rotation,
\item $u_y(\textit{top})=u_z(\textit{top})=0$.
\end{itemize}

\subsection{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~\ref{fig:BC_3D}) and the applied macroscopic stress $\Sigma$ (in Pa), the system of equations becomes:
\begin{equation}
-\tilde{\nabla}\mathbf{\tilde{S}}+\tilde{\nabla} \tilde{p}=0
\end{equation}
\begin{equation}\label{eq:ad_cti_dot_new}
\begin{aligned}
\mathbf{\tilde{S}} &=\tilde{\eta}^\star_n \left(2\alpha_1 \mathbf{\tilde{D}}+2\alpha_2\bfM^DTr\left(\mathbf{M\tilde{D}}\right)+\alpha_3\left(\mathbf{M\tilde{D}}+\mathbf{\tilde{D}M}\right)^D \right) \\
\tilde{\eta}^\star_n &= 2 \left(\alpha_1 Tr\left(\mathbf{\tilde{D}}^2\right)+\alpha_2 Tr\left(\mathbf{M\tilde{D}}\right)^2+\alpha_3 Tr\left(\mathbf{M\tilde{D}}^2\right)\right)^\frac{1-n}{2n}
\end{aligned}
\end{equation}
\begin{equation}\label{eq:ad_c_dot_new}
\frac{\partial c}{\partial\tilde{t}}=\mathbf{\tilde{W}}\left(\mathbf{\tilde{u}}\right).\bfc-\lambda \left[\mathbf{\tilde{D}}\left(\mathbf{\tilde{u}}\right).\bfc-\left(\bfc^T.\mathbf{\tilde{D}}\left(\mathbf{\tilde{u}}\right).\bfc\right).\bfc \right] + \mathcal{M}_o \left(\mathbf{c_0}-\bfc\right)
\end{equation}
where the superscript $\tilde{(.)}$ is used for all non-dimensional variables and operators.

The dimensionless number that arises when scaling the $\bfc$-axis evolution equation is $\mathcal{M}_o$:
\begin{equation}
\mathcal{M}_o=\left(\frac{\eta_n}{\Sigma}\right)^n\frac{1}{\grx}
\end{equation}

$\mathcal{M}_o$ is the ratio between the deformation characteristic time $\tau=(\frac{\eta_n}{\Sigma})^n$ and the recrystallization characteristic time \grx, which controls the rate of rotation of $\bfc$ 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 $\bfc$-axis.

\subsection{Numerical algorithm}

Because the evolution of the $\bfc$-axis orientation depends on the strain rate ($\mathbf{\tilde{D}}(\mathbf{\tilde{u}})$) and the effective viscosity $\eta^\star$ depends on both the strain rate and the orientation of the $\bfc$-axis, the system of coupled \cref{eq:ad_c_dot_new,eq:ad_cti_dot_new} is non-linear. To solve this problem, we linearize the equation describing the evolution of the $\bfc$-axis using a second-order backward differentiation formula (BDF2). The time discretization of the \ctirx model results in the following system of equations:
\begin{equation}
-\tilde{\nabla}\mathbf{\tilde{S}}^{t+1}+\tilde{\nabla} \tilde{p}^{t+1}=0
\end{equation}
\begin{equation}\label{eq:CTI_td}
\begin{aligned}
\mathbf{\tilde{S}}^{t+1} &=\tilde{\eta}^{\star~t+1}_{n} \left(2\alpha_1 \mathbf{\tilde{D}}^{t+1}+2\alpha_2{\bfM^{t+1}}^{D}Tr\left(\bfM^{t+1}\mathbf{\tilde{D}}^{t+1}\right)+\alpha_3\left(\bfM^{t+1}\mathbf{\tilde{D}}^{t+1}+\mathbf{\tilde{D}}^{t+1}\bfM^{t+1}\right)^D \right) \\
\tilde{\eta}^{\star~t+1}_n &= 2 \left(\alpha_1 Tr\left(\left(\mathbf{\tilde{D}}^{t+1}\right)^{2}\right)+\alpha_2 Tr\left(\bfM^{t+1}{\mathbf{\tilde{D}}^{t+1}}\right)^2+\alpha_3 Tr\left(\bfM^{t+1}\left(\mathbf{\tilde{D}}^{t+1}\right)^{2}\right)\right)^\frac{1-n}{2n}
\end{aligned}
\end{equation}
\begin{multline}\label{eq:BDF2}
\frac{3}{2 \Delta \tilde{t}}\left(\bfc^{t+1}-\frac{4}{3}\bfc^t+\frac{1}{3}\bfc^{t-1}\right)\\
 =\mathbf{\tilde{W}}\left(\mathbf{\tilde{u}}^{t+1}\right).\bfc^{t+1}-\lambda \left[\mathbf{\tilde{D}}\left(\mathbf{\tilde{u}}^{t+1}\right).\bfc^{t+1}-\left({\bfc^{t+1}}^T.\mathbf{\tilde{D}}\left(\mathbf{\tilde{u}}^{t+1}\right).\bfc^{t+1}\right).\bfc^{t+1} \right]
+ \mathcal{M}_o \left(\mathbf{c_0}^{t+1}-\bfc^{t+1}\right)
\end{multline}

Finite elements and variational methods were used to solve the time-discretized problem on a Lagrangian grid within the C++ environment Rheolef~\cite{Saramito2022}. The domain is spatially discretize using hexahedral finite elements (unless otherwise mentioned), and an initial $\bfc$-axis orientation is prescribed for each element. The FE approximation for the tensorial mechanical parameter $\bfM$ and the scalar mechanical parameter $\bfc$ is computed with polynomials of the order $0$ ($P_0$). The FE approximation for the tensorial strain rate $\mathbf{\tilde{D}}$, spin rate $\mathbf{\tilde{W}}$ and stress $\mathbf{\tilde{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 \emph{problem\_mixed} solver is used to solve the displacement field $\tilde{\bfu}$ and the pressure field $\tilde{p}$ under the incompressibility hypothesis for the linearized version of \cref{eq:CTI_td}, which is obtained by fixing the apparent viscosity $\eta^\star$. A first fix loop is used to solve the non-linearity associated with $\eta^\star$. This loop is nested in a second fix loop that solves the temporal evolution of the $\bfc$-axis orientation using \cref{eq:BDF2}. In practice, it is recommended for numerical efficiency to flatten the two nested fix loops.

\Cref{alg:algo} presents the numerical algorithm of two nested fix points implemented in \rrrice.

This algorithm is a semi-implicit solver for the coupled \cref{eq:ad_cti_dot_new,eq:ad_c_dot_new}. The simulated orientation field $\bfc$ is thus consistent with the displacement field $\mathbf{\tilde{u}}$ and therefore with the diatoric strain rate $\mathbf{\tilde{D}}$ and the deviatoric stress $\mathbf{\tilde{S}}$. This implementation of \ctirx model is distributed as the \rrrice code.

\algnewcommand{\Inputs}[1]{
\State \textbf{Inputs:}
\Statex \hspace*{\algorithmicindent}\parbox[t]{.8\linewidth}{\raggedright #1} }
\algnewcommand{\Initialize}[1]{
\State \textbf{Initialize:}
\Statex \hspace*{\algorithmicindent}\parbox[t]{.8\linewidth}{\raggedright #1}}
%\algnewcommand{\Inputs}[1]{%
%\State \textbf{Inputs:}
%\Statex \hspace*{\algorithmicindent}\parbox[t]{.8\linewidth}{\raggedright #1}
%}
%\algnewcommand{\Initialize}[1]{%
%\State \textbf{Initialize:}
%\Statex \hspace*{\algorithmicindent}\parbox[t]{.8\linewidth}{\raggedright #1}
%}

\begin{Algorithm}
\begin{algorithmic}[1]
\Inputs {$\mathbf{c}^{init}$}
\Initialize {$\mathbf{c}^{-1} \gets \mathbf{c}^{init}$ \\ $\mathbf{c}^{0} \gets \mathbf{c}^{init}$ \\ $\mathbf{u}^0 \gets 0$}
\For {$t$ from $0$ to $N$} \Comment{Time evolution loop}
\State $k_c \gets 0$
\State $res_{c} \gets +\infty$
\State $\mathbf{u}^{t+1,k_c=0,k_{nl}=0} \gets \mathbf{u}^{t}$
\State $\mathbf{c}^{t+1,k_c=0} \gets \mathbf{c}^{t}$
\While{$res_c > tol_c$} \Comment{Fix loop to solved (\cref{eq:BDF2})}
\State $k_{nl} \gets 0$
\State $res_{nl} \gets +\infty$
\State $\mathbf{u}^{t+1,k_c,k_{nl}=0} \gets \mathbf{u}^{t+1,k_c}$
\While{$res_{nl}> tol_{nl}$} \Comment{Fix loop to solve non-linear CTI (\cref{eq:ad_cti_dot_new})}
\State $H^\star_n \gets \eta^\star_n(\mathbf{u}^{t+1,k_c,k_{nl}},\mathbf{c}^{t+1,k_c})$ \Comment{\cref{eq:ad_cti_dot_new}}
\State find $\mathbf{u}^{t+1,k_c,k_{nl}+1} \gets \mathbf{u}^{t+1,k_c,k_{nl}+1}(H^\star_n,\mathbf{c}^{t+1,k_c})$ \Comment{\emph{problem\_mixed} solver on \cref{eq:ad_cti_dot_new}}
\State $res_{nl} \gets \| \mathbf{S}(\eta^\star_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^\star_n,\mathbf{u}^{t+1,k_c,k_{nl}+1},\mathbf{c}^{t+1,k_c})\|$
\State $k_{nl} \gets k_{nl}+1$
\EndWhile
\State $H^\star_n \gets \eta^\star_n(\mathbf{u}^{t+1,k_c,k_{nl}},\mathbf{c}^{t+1,k_c})$
\State compute $\mathbf{S} \gets \mathbf{S}(H^\star_n,\mathbf{u}^{t+1,k_c,k_{nl}})$ \Comment{\cref{eq:ad_cti_dot_new}}
\State find $\mathbf{c_0} \gets \mathbf{c_0}(\mathbf{S})$ \Comment{\cref{eq:c0_solution}}
\State compute $\mathbf{c}^{t+1,k_c+1} \gets \mathbf{c}^{t+1,k_c+1}(\mathbf{c_0},\mathbf{u}^{t+1,k_c,k_{nl}})$ \Comment{\cref{eq:BDF2}}
\State $res_c \gets \| \mathbf{S}(\eta^\star_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^\star_n,\mathbf{u}^{t+1,k_c,k_{nl}},\mathbf{c}^{t+1,k_c})\|$
\State $k_c \gets k_c+1$
\EndWhile
\State $\mathbf{u}^{t+1} \gets \mathbf{u}^{t+1,k_c,k_{nl}}$
\State $\mathbf{c}^{t+1} \gets \mathbf{c}^{t+1,k_c}$
\EndFor
\end{algorithmic}
\caption{Numerical algorithm, based on two nested fix loops: one for solving the viscosity, $\eta^\star_n$, nested in a second one solving the temporal evolution of the orientation of $\bfc$-axis.}\label{alg:algo}
\end{Algorithm}


\subsection{ CTI parameter values for the ice single crystal }\label{subsec:cti_param}

As stated in \cref{sec:num_model}, the values for the CTI parameters, $\gamma$, $\eta_n, n$ and $\beta$ (\cref{eq:CTI_law,eq:alpha_i}) 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~\cite{Meyssonnier1996,Mansuy2002}.

The following Duval et al. (1983) \cite{Duval1983} and previous modeling work based on the CTI law~\cite{Mansuy2002,Martin2009,Rathmann2021}, 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$~\cite{Duval1983}. 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~\cite{Duval1983}, the ratio between basal and non-basal shear viscosities is on the order of $10^4$. This ratio is related to $\beta$ through~\cref{eq:SC_shear_12}, $\beta^{\frac{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~\cite{Mansuy2002,Martin2009,Rathmann2021}.

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}^{sim}_{sc}=1.7\times10^{-2}(\frac{\eta_n}{\Sigma})^n$ for the single crystal and $\tilde{D}^{sim}_{px}=5.3\times10^{-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~\cite{Duval1983} ($D^{exp}_{sc}=(2\times10^{-5},2\times10^{-4})$ s$^{-1}$ for the single crystal, $D^{exp}_{px}=(9\times 10^{-8},3\times 10^{-7})$ s$^{-1}$ for the polycrystal). $\eta_3^{sc}=(4.4,9.4)$ MPa s$^{1/3}$ for the single crystal and $\eta_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 $\eta_3$ in order to reproduce the experiment duration as observed during uniaxial compression tests that are taken here as a reference (see~\cref{sec:comp-Montagnat}). The parameters of the CTI law for the ice single crystal are summarized in~\cref{tab:cti_param}.

\begin{table}[]
\caption{Parameters of the CTI law for the ice single crystal.}
\centering
\begin{tabular}{c|c|c|c}
$n$ & $\gamma$ & $\beta$ & $\eta_n$ \\
\hline $3$ & $1$ & $10^{-2}$ & $7.5~MPa.s^{1/3}$
\end{tabular}
\label{tab:cti_param}
\end{table}


\section{Validation of the CTI-based implementation for predicting viscoplastic response}\label{sec:validation_r3ice}

In this section, we verify that the CTI law implemented in \rrrice 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~\cite{Grennerat2012}. We first compare the \rrrice predictions to those of the CraFT-EVP model~\cite{Grennerat2012}, which has been shown to adequately reproduce the heterogeneous strain fields observed in experiments~\cite{Grennerat2012, Piazolo2015}. 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.
\subsection{Stress and strain field predictionscompared with CraFT-EVP}\label{sec:EVPvsCTI}

CraFT-EVP has been used to provide full-field predictions of the mechanical response of ice during primary and secondary creep~\cite{Grennerat2012}. It is based on an elasto-viscoplastic formulation proposed by Suquet et al. 2012~\cite{Suquet2012} 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~\cite{Suquet2012, Grennerat2012}. 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 \rrrice both simulate a stationary secondary creep.

Reproducing accurate stress fields in \rrrice is key because the recrystallization attractor $\mathbf{c_0}$ depends on the deviatoric stress $\bfS$ at the element scale (\cref{subsec:c0_attractor}). Note that the Elasto-ViscoPlastic (EVP) law implemented in CraFT-EVP~\cite{Suquet2012} 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 \rrrice is a simplified parameterization of viscoplastic anisotropy, which has been calibrated directly from experimental data (\cref{subsec:cti_param}). The objective of the comparison performed below is to verify that the stress field predicted by \rrrice 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. \rrrice parameters were, of course, not fitted to the CraFT-EVP viscoplastic response, so the two predictions can be considered independent.

\Cref{fig:mesh_grain}$\MK$(a-b) shows the polycrystal microstructure generated by a grain growth algorithm using Neper~\cite{Quey2011} used for both CraFT-EVP and \rrrice 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 \rrrice polycrystalline microstructure used here was meshed with a tetrahedral geometry. The resolution of CraFT-EVP simulations was reduced to match the \rrrice one. \Cref{fig:EVPvsCTI_field} shows the component $\varepsilon_{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 \rrrice predictions although the overall patterns do not show a strong mismatch. The locations of areas in compression (blue) and extension (red) are similarly prescribed (\Cref{fig:EVPvsCTI_field}$\MK$(a)). Significant differences between the responses of the models CraFT-EVP and \rrrice 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 \rrrice, 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 \Cref{fig:EVPvsCTI_field}$\MK$(a). In CraFT-EVP, periodic boundary conditions are imposed because equations are solved with a fast Fourier transform-based method. The simulated stress fields (\Cref{fig:EVPvsCTI_field}$\MK$(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 (\cref{fig:EVPvsCTI}). These discrepancies are stronger between the two configurations in which grains are described by multiple elements (\cref{fig:EVPvsCTI}) and where the impact of the boundary conditions mentioned above is stronger. Distributions for all components of both tensors $\left(\varepsilon_{ij},\sigma_{ij}\right)$ are given in \cref{sec:dij_comp}.
\break

In \rrrice, the stress tensor at the element scale is used to predict the recrystallization attractor $\mathbf{c_0}$. The latter is based on the principal direction of the stress tensor $\mathbf{v_1}$ and $\mathbf{v_3}$ (\cref{subsec:c0_attractor}). \Cref{fig:S_e1_comp} shows the distribution of $\mathbf{v_1}$ and $\mathbf{v_3}$ orientations predicted using both CraFT-EVP (\cref{fig:S_e1_comp}$\MK$(a)), and \rrrice with a configuration with several elements per grain (\cref{fig:S_e1_comp}$\MK$(b)) and one orientation per element (\cref{fig:S_e1_comp}$\MK$(c)). All three configurations predict similar $\mathbf{v_1}$ and $\mathbf{v_3}$ orientation distributions and therefore similar $\mathbf{c_0}$ 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 \rrrice. In particular, the principal stress tensor directions used in \rrrice to simulate the dynamic recrystallization impact on orientation evolution are similarly predicted by the two models.


\begin{figure}[!ht]
\centering
\includegraphics[width=14cm]{figures/micro.png}
\vspace*{-10pt}
\caption{Microstructure generated via a grain growth using Neper meshed for (a) CraFT-EVP microstructure with 265302 grid points (b) \rrrice with 262688 tetrahedral elements. (c) Simplified microstructure in which orientations are defined at the element scale with a hexahedral mesh of 2740 elements.}\label{fig:mesh_grain}
\end{figure}

\begin{figure}[!ht]
\centering
\includegraphics[width=12cm]{figures/field.png}
\caption{Comparison of predicted (a) strain ($\varepsilon_{yy}$) (b) stress ($\sigma_{yy}=S_{yy}+p$) fields for the grain microstructure using CraFT-EVP and \rrrice.}\label{fig:EVPvsCTI_field}
\end{figure}

\begin{figure}[!ht]
\centering
\includegraphics[width=10cm]{figures/DSeq_comparaison.png}
\caption{Kernel Density Estimation of the Von Mises equivalent strain (a) and stress (b) fields simulated with CraFT-EVP (EVP, black line) and \rrrice (CTI, cyan line). The dashed lines show the results obtained for simulations where grains are defined by multiple mesh elements (\emph{grains} microstructure) (\cref{fig:mesh_grain}$\MK$(a-b)). The full line shows a \rrrice simulation where a single orientation is attributed per element (\cref{fig:mesh_grain}$\MK$(c)).}\label{fig:EVPvsCTI}
\end{figure}



\subsection{Stress and strain field prediction in the simplified mesh configuration}

Resolving \ctirx 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 (\cref{fig:mesh_grain}$\MK$(c)) in order to keep a statistically reasonable number of orientations on a mesh size ($\sim 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 \rrrice models (\cref{fig:EVPvsCTI}). 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 \ctirx fields predictions 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~\ref{fig:S_e1_comp} compares the distribution of the principal stress tensor directions $\mathbf{v_1}$ and $\mathbf{v_3}$ that are used in the calculation of the dynamic recrystallization attractor in \rrrice. There is a superb correspondence between the distributions predicted by CraFT-EVP and by \rrrice 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 $\mathbf{c_0}$ at the sample scale.



\section{Texture evolution: comparison with experimental data}\label{sec:resultat}

This section tests the ability of the \rrrice model to predict the evolution of the $\bfc$-axes orientation with and without dynamic recrystallization (\cref{eq:ad_c_dot_new}). 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~\cite{De-la-chapelle1998}. These textures were well reproduced by uniaxial compression simulations performed using the mean-field ViscoPlastic Self Consistent model~\cite{Castelnau1996c, Montagnat2012}.
\rrrice prediction for texture evolution without recrystallization are therefore compared with those of VPSC simulations extracted from~\cite{Montagnat2012}. 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 \rrrice for the effect of DRX on texture evolution and mechanical behavior under uniaxial compression, uniaxial tension, and simple shear.

\break
The laboratory creep experiments selected to validate the model are:
\begin{itemize}
\item Uniaxial compression experiments from Montagnat et al. 2015~\cite{Montagnat2015} ($\Sigma\sim 0.7$~MPa,\newline T$=-5^\circ$ C)
\item Uniaxial tension experiments from Jacka and Maccagnan 1984~\cite{Jacka1984b} ($\Sigma \sim 0.4$~MPa,\newline T$=-3^\circ$ C)
\item Torsion experiments from Journaux et al. 2019~\cite{Journaux2019} ($\Sigma_{max}\sim 0.5$~MPa, T=-7$^\circ$ C),
\end{itemize}
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~\cite{Qi2017}. 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 \cref{tab:simulation_param}.



For all configurations, the simulated texture evolution is represented by (1) the evolution of the eigenvalues $a_i$ of the second-order orientation tensor $\bfa^{(2)}$ ($\mathbf{a^{(2)}}=\frac{1}{N}\Sigma_{i=0}^N\mathbf{c_i}\otimes\mathbf{c_i}$, with $N$ being the number of orientations), (2) pole figures representing the $\bfc$-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.


\begin{figure}[!hb]
\centering
\includegraphics[width=12cm]{figures/e1_comparaison.png}
\caption{Distributions of the principal stress tensor directions $\mathbf{v_1}$ and $\mathbf{v_3}$ for uniaxial compression on (a) CraFT-EVP simulation (b) \rrrice simulation with grains (c) \rrrice simulation with one orientation per mesh element.}\label{fig:S_e1_comp}
\end{figure}

\begin{table}[!hb]
\caption{Values of the numerical and physical parameters used in the model \rrrice.}
\centering
\begin{tabular}{c|c|c|c|c|c}
References & \multicolumn{2}{c}{Model inputs} & \multicolumn{3}{c}{Physical values} \\
& Stress & $\mathcal{M}_o$ & $\Sigma~(MPa)$ & $\grx~(\emph{days})$ & T $({}^\circ C)$ \\
\hline without recrystallization & Compression & $0$ & $1$ & $+\infty$ \\
Montagnat et al. 2015~\cite{Montagnat2015} & Compression & $1.03\times 10^{-3}$ & $0.7$ & $13.8$ & $-5$\\
Jacka et al. 1984~\cite{Jacka1984b} & Tension & $1.86\times 10^{-2}$ & $0.4$ & $4.1$ & $-3$\\
Journaux et al. 2019~\cite{Journaux2019} & Simple Shear & $1.06\times 10^{-3}$ & $0.5$ & $36.8$ & $-7$ \\
Qi et al. 2017~\cite{Qi2017} & Compression & $\left(1.34,20.8\right)\times10^{-5}$ & $\left(1.3,~4.3\right)$ & $4.6$ & $-10$
\end{tabular}
\label{tab:simulation_param}
\end{table}


\break
The second-order orientation tensor $\bfa^{(2)}$ is a statistical representation of texture ani\-sotropy~\cite{Woodcock1977}. The $\bfc$-axes distribution lies within an ellipsoid whose axes are the eigenvalues of $\bfa^{(2)}$ with $1 \geq a_1 \geq a_2 \geq 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.

\subsection{Texture evolution in the absence of recrystallization}\label{sec:def-text}


To test the capability of \rrrice 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.

\Cref{fig:C_Mo0}$\MK$(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). \Cref{fig:C_Mo0}$\MK$(d) shows textures predicted at $22\%$, $40\%$ and $60\%$ total strain.
The strain-rate curve discontinuities (\cref{fig:C_Mo0}$\MK$(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 $\bfc$-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 $\sim 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~\cite{Montagnat2012} 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.~\cite{Montagnat2012} using the VPSC model when departing from an isotropic texture (see red dots in \cref{fig:C_Mo0}$\MK$(a)).

\begin{figure}[!ht]
\includegraphics[height=11cm]{figures/CS_0.png}
\caption{\label{fig:C_Mo0}
Simulation of uniaxial compression along $Y$ with no recrystallization $\grx=+\infty$, $\mathcal{M}_o=0$. a) Evolution with strain of the eigenvalues of the second-order orientation tensor ($a_1>a_2>a_3$). Red dots are predictions from the VPSC model taken from~\cite{Montagnat2012}; b) evolution with strain of the macroscopic strain rate; c) evolution 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 $22\%$, $40\%$ and $60\%$ macrostrain are shown. The white dashed lines define the angle at $\pm45^\circ$ 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). }
\end{figure}

\subsection{Uniaxial compression}\label{sec:comp-Montagnat}

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 conpression direction) predicted for uniaxial compression with DRX under an imposed compressive stress of $0.7$~MPa is shown in \cref{fig:C_MoW}. This is compared with experimental measurements from Montagnat et al. 2015~\cite{Montagnat2015}. These latter were obtained from initially isotropic granular ice compressed at $-5^\circ 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 $\bfc$-axis pole figures shows good qualitative agreement between the simulated and experimental textures (\cref{fig:C_MoW}$\MK$(d-e)). Both are girdle textures that become more concentrated and show an evolution of the cone angle to the compression axis from $45^\circ$ to $35^\circ$ up to $\sim20\%$ macroscopic strain and remain nearly constant thereafter (\cref{fig:C_MoW}$\MK$(d)). There is also good agreement between the predicted evolution of the texture eigenvalues and those measured at $7\%$, $12\%$ and $17.8\%$ (\cref{fig:C_MoW}$\MK$(a)).
This good agreement between the experimental and simulated textures was obtained for a value of $\grx=13.8$ days.
During the experiments, $17.8\%$ strain was reached in 144 hours ($\sim 6$ days). In the model simulations, this strain is reached in a shorter time ($\sim 130$ hours), but the tertiary creep (maximum strain-rate) regime is reached between 10 to 20\% strain, similar to the experimental observations (see also~\cite{Jacka1984b, Treverrow2012}).

The macroscopic response predicted by \rrrice shows a clear softening associated with the development of the texture. A maximum enhancement factor (EF) of $4.7$ is obtained at $\sim 19\%$ strain (black arrow, \cref{fig:C_MoW}), which corresponds to a girdle texture with a cone angle around $\sim 45^\circ$ to the compression axis.
This softening is followed by a slight hardening as the girdle texture slowly evolves toward a cone angle at $\sim 35^\circ$ 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~\cite{Montagnat2015} 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{\varepsilon}_{sim}^{corr}=\dot{\varepsilon}_{sim}\times\left(1+\varepsilon\right)^3$. This corrected strain rate is plotted in blue in \cref{fig:C_MoW}$\MK$(b).
This brings the maximum strain rate to $\sim 10\%$ and the enhancement factor to $3.0$.
These values are closer to the macroscopic response measured by Montagnat et al.~\cite{Montagnat2015}, which is characterized by an enhancement factor between $3.3$ and $3.8$.

\begin{figure}[!ht]
\includegraphics[height=9cm]{figures/Montagnat2015.png}
\vspace*{-5pt}
\caption{\label{fig:C_MoW}
Simulations for uniaxial compression creep along the $Y$ direction with parameters from \cref{tab:simulation_param} and comparison with experimental data. a) Evolution with strain of the eigenvalues of the second-order orientation tensor ($a_1>a_2>a_3$). The blue dots correspond to experimental values from~\cite{Montagnat2015}; 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 $\theta$ 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 $\pm45^\circ$ 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~\cite{Montagnat2015}. All figures presenting texture data are color-coded using the same color bar (scales in multiples of a uniform distribution).}
\end{figure}


\subsection{Uniaxial tension}


Experimental creep tests performed under uniaxial tension are scarce.
Jacka et al. 1984~\cite{Jacka1984b} presented one experiment performed at $-3^\circ$ C, up to $9.3\%$ octahedral shear strain that corresponds to the definition given in Jacka et al. 1984~\cite{Jacka1984a}, 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 \cref{fig:TextureJacka1984}. The texture is characterized by a small circle girdle with a mean half angle of $50.4^\circ$ and a standard deviation of $15.3^\circ$. Although tertiary creep was not reached in this experiment, an enhancement factor of $\sim 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 \cref{fig:T_MoW}.
The texture evolves with increasing strain toward a girdle texture with a cone angle slightly around $45^\circ$ to the extension direction. This texture appears to stabilize at $\sim 5\%$ axial strain. The distribution of the $\bfc$-axis angle from the tension axis is given in \cref{fig:jacka1984_hist_ii} for the simulated texture at $13\%$ strain. It is similar to the one measured by Jacka et al. 1984~\cite{Jacka1984b} (\cref{fig:TextureJacka1984}).
An enhancement factor of approximately $10$ is simulated at 10\% axial strain.

\begin{figure}[!ht]
\centering
\begin{minipage}[t]{0.47\textwidth}\centering
\includegraphics[width=0.7\columnwidth]{figures/Jacka1984_Texture.png}
\end{minipage}
\qquad
\begin{minipage}[t]{0.47\textwidth}\centering
\includegraphics[width=0.7\columnwidth]{figures/Jacka1984_distribution.png}
\end{minipage}
\caption{Texture obtained after $9\%$ of octahedral tension strain (about 13\% axial strain) performed at $-3^\circ C$ with an octahedral stress of $0.4~MPa$. The star $\star$ in the center of the pole figure shows the direction of the tension axis. From Jacka et al. 1984~\cite{Jacka1984b}.}\label{fig:TextureJacka1984}
\end{figure}


\begin{figure}[!ht]
\includegraphics[height=11cm]{figures/Jacka1984.png}
\caption{Simulations for uniaxial tension creep along the $Y$ direction with parameters from \cref{tab:simulation_param} and comparison with experimental data. a) Evolution with strain of the eigenvalues of the second-order orientation tensor ($a_1>a_2>a_3$); 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 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 $\pm45^\circ$ 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).}
\label{fig:T_MoW}
\end{figure}

\begin{figure}
\centering
\includegraphics[height=5cm]{figures/Jacka1984_hist.png}
\caption{Distribution of the angle between the $\bfc$-axis and the $\bfy$-axis (tension axis) corresponding to pole figure \emph{(ii)} in \cref{fig:T_MoW}$\MK$(d).}\label{fig:jacka1984_hist_ii}
\end{figure}

\subsection{Simple shear}

The simulated texture evolution in simple shear is compared with observations made by Journaux et al.~\cite{Journaux2019} for torsion experiments on isotropic granular ice performed at $-7^\circ$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~\cite{Bouchez1982, Burg1986, Qi2019}. As the shear strain increases, $M2$ merges with $M1$ to form a highly concentrated cluster texture characterized by a single maximum of $\bfc$-axes normal to the shear plane. The simple shear simulations performed here reproduce this evolution relatively well, as shown in \cref{fig:S_MoW}.

By accounting for dynamic recrystallization, the \rrrice model is able to reproduce the evolution of both $M1$ and $M2$ maxima with strain. At low shear strains ($\gamma \sim 0.2$), $M1$ concentrates normal to the shear plane, while $M2$ follows the principal extension direction (white dashed line, \cref{fig:S_MoW}$\MK$(c) bottom and \cref{fig:S_MoW}$\MK$(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 (\cref{fig:S_MoW}$\MK$(a)).

In both experiments and simulations, the strain rate increased rapidly up to $\gamma=0.3$ while the double-maxima texture develops (\cref{fig:S_MoW}$\MK$(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 $\gamma\sim3$ where is reaches a quasi-steady state. As explained in \cref{sec:comp-Montagnat}, the discontinuities in the strain-rate curve (\cref{fig:S_MoW}$\MK$(b)) are technical artifacts imposed by the operating rules of GRICAD clusters.

In phase with Treverrow et al. 2012~\cite{Treverrow2012}, 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~\cite{Treverrow2012} data for a similar applied shear stress.


\begin{figure}[!ht]
\includegraphics[scale=0.40]{figures/Journaux2019.png}
\caption{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 \cref{tab:simulation_param} and comparison with experimental data from~\cite{Journaux2019}. 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~\cite{Journaux2019}; 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~\cite{Journaux2019}. All figures presenting texture data are color-coded using the same colorbar (scales in multiples of a uniform distribution).}
\label{fig:S_MoW}
\end{figure}


\subsection{Testing the relative effects of deformation and recrystallization on texture evolution under high stresses}

The work of Qi et al.~\cite{Qi2017} 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 \cref{tab:simulation_param}).
The stress values reported in \cref{tab:simulation_param} 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~\cite{Qi2017} experiments. As \ctirx 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~\cite{Mellor1982}. 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 \cref{fig:S_QiW}.
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 $\sim 20\%$.
To correctly fit the experimental data, a value of \grx of $4.6$ days is used, which is lower than that used to reproduce the unconfined compression creep tests of Montagnat et al. 2015~\cite{Montagnat2015} (\cref{tab:simulation_param}). This implies a higher contribution of recrystallization compared with deformation at higher stresses.

\begin{figure}[!ht]
\includegraphics[width=13cm]{figures/Qi2017.png}
\caption{Simulation of the uniaxial compression creep test along $Y$ at various stresses. Comparison between \rrrice predictions and experiments by Qi et al.~\cite{Qi2017}. (a) Textures from~\cite{Qi2017} obtained by EBSD for various strain rates leading to different quasi-steady state compressive flow stresses at macroscopic strains of $\sim0.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).}
\label{fig:S_QiW}
\end{figure}


\section{Discussion}

The \ctirx 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.

The novelty of the present model and its implementation is based on the simulation of the evolution of the texture (orientations of $\bfc$-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~\cite{Grennerat2012, Piazolo2015, Chauve2015} and (ii) revealed that the formation of new orientations (nucleation) during dynamic recrystallization occurs through bulging (heterogeneous grain boundary migration) and polygonization (formation of new grains boundaries by organization of dislocations within a grain) \cite{Piazolo2015, Chauve2015, Chauve2017, Chauve2017b, Chauve2017c}. 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 $\bfc$-axis orientation evolution (\cref{eq:c_dot_new}), with a targeted recrystallized orientation that is calculated to maximize the local basal resolved shear stress (\cref{eq:c0_solution}).

Comparison with tertiary creep experiments performed in compression~\cite{Treverrow2012, Piazolo2013, Montagnat2015}, in simple shear~\cite{Treverrow2012, Qi2019, Journaux2019} and in tension~\cite{Jacka1984b} supports the fact that the \rrrice model can reproduce texture evolution in conditions similar to those of laboratory experiments ($T=-5$ to $-10^\circ$C, $S\sim 1$~MPa). For this, a single parameter is adjusted: the recrystallization parameter \grx, which weights the relative contribution of deformation and recrystallization to the rotation of the $\bfc$-axes. As expected, under similar stress conditions, a lower temperature requires a higher \grx 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 \rrrice, we show that the \ctirx 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~\cite{Jacka1984b, Li1996, Treverrow2012, Budd2013, Montagnat2015}. 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~\cite{Qi2017}. 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.
\rrrice 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 \rrrice providing that \grx is given a relatively low value ($4.6$ days for $T=-10\degree C$) compared to the one determined for unconfined creep compression experiments ($13.8$ days for $T=-5\degree$C).
This low value of \grx 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 $\bfc$-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.~\cite{Kalifa1992}, which performed triaxial tests under the same range of strain rate ($\sim 10^{-5}-10^{4}$~s$^{-1}$) than Qi et al. 2017~\cite{Qi2017}, 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\cite{Qi2017} experiments. Strong stress concentration occurs at crack tips, and Chauve et al. 2017~\cite{Chauve2017c} have recently shown that, for ice, this stress concentration can be released through DRX. The \rrrice 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 (\grx).



Recent models of texture evolution during dynamic recrystallization have been proposed by Rathmann et al. 2021~\cite{Rathmann2021b} and Richard et al. 2021,2022~\cite{Richards2021,Richards2022}, 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 \emph{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 \ctirx formulation. The \emph{deformability} term can be obtained from \cref{eq:c_dot_new} 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 \rrrice, 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 \rrrice.
This highlights the fundamental role of stress heterogeneities in texture evolution during the transition between secondary and tertiary creep (from $1$ to $\sim 15\%$ strain).

The major strength of~\cite{Rathmann2021b,Richards2021,Richards2022} 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 \rrrice model could be a good candidate for providing constraints to develop parameterizations for the temporal evolution, destined to inexpensive formulations.

The \ctirx 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 \cref{eq:CTI_law}. For instance, the versatility of the \ctirx 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.

\section{Conclusions}

In this paper, we present a new formulation, \ctirx, 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 $\mathbf{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 $\mathbf{c_0}$ such as to maximize the local resolved shear stress, thus providing a physically based formulation for recrystallization-induced $\bfc$-axes rotation. This formulation is as simple as possible, given our knowledge of the mechanisms.

The \rrrice model is the result of the implementation of the \ctirx formulation in a finite element framework. In the model, the texture results from a balance between $\bfc$-axis rotation due to viscoplastic deformation and dynamic recrystallization, which is controlled by \grx, the only tuning parameter of the model. The \rrrice 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 \rrrice 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 \grx required to fit the texture measured during the confined compression experiments~\cite{Qi2017} relative to that used to fit the low stress unconfined compression creep experiments~\cite{Montagnat2015}, despite the higher experimental temperature of the latter ($-5^\circ$ C vs. $-10^\circ$ C, cf. \cref{tab:simulation_param}). These results provide an additional illustration of the ability of the \rrrice model to help resolve open questions about ice deformation behavior.

The computational cost of the \rrrice 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.


\section*{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.




\section*{Acknowledgements}
The authors are thankful to Chao Qi for quickly providing the EBSD data from \cite{Qi2017} to test \rrrice abilities during deformation under confining pressure and high stresses.

All computations presented in this paper were performed using the GRICAD infrastructure (\url{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.

\section*{Authors' contributions} 
Thomas Chauve conducted most of the study, developed the model, performed the numerical implementation of \rrrice, 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.

\section*{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 \rrrice code is available under the GPL-v3 license. Gitlab repository \url{https://gricad-gitlab.univ-grenoble-alpes.fr/mecaiceige/tools/ice-polycrystal-models/rheolef_cti}. The associated documentation is available at \url{https://mecaiceige.gricad-pages.univ-grenoble-alpes.fr/tools/ice-polycrystal-models/ipms_documentation/Intro.html}

\CDRGrant[ERC RhEoVOLUTION]{882450}
\CDRGrant[ERC IVORI]{949516}


\appendix

\section{Demonstration of \texorpdfstring{$\delta_i$}{δi} identification}\label{sec:demo_di}
\subsection{Constitutive equation CTI and single crystal behavior}
\def\theequation{\Alph{section}\arabic{equation}}
\setcounter{equation}{0}
\def\thetable{\Alph{section}\arabic{table}}
\setcounter{figure}{0}
The constitutive equation is
\begin{equation}\label{Aeq:CTI}
\bfD=\left(\frac{\partial \phi_S^{(1)}}{\partial \bfS}\right)^D=2\delta_1 \bfS + 2\delta_2 \bfM^D Tr(\fMS) + \delta_3 (\fMS+\fSM)^D
\end{equation}

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:
\begin{equation}
{}^{g}D_{13}=\frac{\psi_1}{2}{}^{g}S_{13},\quad{}^{g}D_{23}=\frac{\psi_1}{2}{}^{g}S_{23}
\end{equation}

$\beta$ is the viscosity ratio between shear parallel to the basal plane $\eta$ and the viscosity for shear within the basal plane:
\begin{equation}
^{g}D_{12}=\beta\frac{\psi_1}{2}{}^{g}S_{12}
\end{equation}

$\gamma$ is the viscosity ratio between compression (or traction) Along the $c$ axis $\eta_c$ and compression (or traction) in one direction ${}^ge_r$ within the basal plane:
\begin{equation}
{}^{g}S_{33}=2\eta_c{}^{g}D_{33}, \quad \gamma {}^{g}S_{rr}=2\eta_c{}^{g}D_{rr}
\end{equation}
The same compressive (or extensive) stress parallel and normal to the $\bfc$-axis $S={}^gS_{33}={}^gS_{rr}$ results therefore in different strain rates:
\begin{equation}
{}^gD_{rr}=\gamma{}^gD_{33}
\end{equation}

\subsection{\texorpdfstring{$\delta_i$}{δi} identification}

To identify the parameters $(\delta_1,\delta_2,\delta_3)$ in the constitutive \cref{Aeq:CTI}, we developed it and compared it with the single crystal behavior.

Using $\bfc=(0,0,1)$
\begin{equation}
\bfM=
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 1
\end{bmatrix},
\quad \bfM^D=\frac{1}{3}
\begin{bmatrix}
1 & 0 & 0 \\
0 & -1 & 0 \\
0 & 0 & 2
\end{bmatrix}, \quad Tr\left(\bfM{}^g\bfD\right)={}^gD_{33}
\end{equation}
\begin{enumerate}
\item $\delta_1$ is given by applying stress
\[
{}^g\bfS=
\begin{bmatrix}
0 & {}^gS_{12} & 0 \\
{}^gS_{12} & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
\]
Then
\begin{equation}
Tr(\fMS)=0, \quad \left(\fMS+\fSM\right)^D=0
\end{equation}

The constitutive equation gives
\begin{equation}
{}^gD_{12}=2\delta_1{}^gS_{12}
\end{equation}


Therefore:
\begin{equation}
\delta_1=\frac{\psi_1\beta}{4}
\end{equation}


\item $\delta_3$ is given by applying stress
\[
{}^g\bfS=
\begin{bmatrix}
0 & 0 & {}^gS_{13} \\
 0 & 0 & 0 \\
{}^gS_{13} & 0 & 0
\end{bmatrix}
\]
Then
\begin{equation}
Tr(\fMS)=0, \quad \left(\fMS+\fSM\right)^D=\bfS
\end{equation}


The constitutive equation gives
\begin{equation}
{}^gD_{13}=(2\delta_1+\delta_3){}^gS_{13}
\end{equation}


Therefore
\begin{equation}
\delta_3=\frac{\psi_1}{2}\left(1-\beta\right)
\end{equation}



\item $\delta_2$ identification is less trivial and presented in more detail.
\begin{itemize}
\item First, uniaxial stress $\tilde{\sigma}$ is applied along the $\bfc$-axis. The stress tensor is ${}^g\mathbf{\sigma}=(\begin{smallmatrix}
0&0&0\\
0&0&0\\
0&0&\tilde{\sigma}
\end{smallmatrix})$
and the associated deviatoric stress tensor is
\[
{}^g\bfS=
\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}
\]
\begin{itemize}
\item $Tr(\bfM{}^g\bfS)=\frac{2\tilde{\sigma}}{3}$
\item $\left(\bfM{}^g\bfS+{}^g\fSM\right)^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}$
\item 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}
\]
\end{itemize}
\item Then a uniaxial stress $\tilde{\sigma}$ is applied normal to the $\bfc$-axis, for instance in the direction~$1$. The stress tensor is
\[
{}^g\mathbf{\sigma}=
\begin{bmatrix}
\tilde{\sigma}
&0&0\\0&0&0\\0&0&0
\end{bmatrix}
\]
and the associated diatoric stress tensor is
\[
{}^g\bfS=
\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}
\]
\begin{itemize}
\item $Tr(\bfM{}^g\bfS)=-\frac{\tilde{\sigma}}{3}$
\item $\left(\bfM{}^g\bfS+{}^g\fSM\right)^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}$
\item 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}
\]
\end{itemize}
\item Finally from the grain behavior, we can write that $\gamma {}^gD_{33}={}^gD_{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
\]
\item This gives:
\[
\delta_2=6\delta_1\frac{1-\gamma}{4\gamma-1}-\delta_3
\]
\item 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)
\]
\end{itemize}
\end{enumerate}


\section{Strain and stress field comparisons}\label{sec:dij_comp}

The \cref{fig:D_comp,fig:S_comp} show the probability distribution functions for all components of the strain and stress tensors during uniaxial creep with $\bfc$-axis evolution (\cref{sec:EVPvsCTI}).

\begin{figure}[ht!]
\centering
\includegraphics[width=11.5cm]{figures/Di_comparaison.png}
\caption{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 \rrrice 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.\label{fig:D_comp}}
\end{figure}

\begin{figure}[ht!]
\centering
\includegraphics[width=11.5cm]{figures/Si_comparaison.png}
\caption{Kernel Density Estimation of all components of the stress tensor predicted by an elasto-viscoplastic (EVP) simulation using CraFT and the continuous transverse isotropic (CTI) using \rrrice after $0.01$ macroscopic strain. The dashed lines show the results for simulations where grains are well discretize. The full line shows the result for the simulation with one orientation per element.\label{fig:S_comp}}
\end{figure}

% \bibliographystyle{crunsrt}
%\bibliography{Chauve}
 \CDRGrant[ERC RhEoVOLUTION]{882450}
\CDRGrant[ERC IVORI]{949516}

\break 
\printbibliography
\end{document}
