Thomas Chauve
email: thomas.chauve@univ-grenoble-alpes.fr

Maurine Montagnat

Véronique Dansereau

Pierre Saramito

Kévin Fourteau

Andrea Tommasi

Maurine Montagna

Véronique Danserea

Andréa Tommasi

Andréa Tommasi

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

Keywords: dynamic recrystallization, texture, viscoplastic anisotropy, finite element method, ice

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Dynamic recrystallization (DRX) is a mechanism by which local stress incompatibilities and local strain heterogeneities are relaxed [START_REF] Humphreys | Recrystallization and related annealing phenomena[END_REF][START_REF] Sakai | Dynamic and post-dynamic recrystallization under hot, cold and severe plastic deformation conditions[END_REF][START_REF] Chauve | Strain field evolution during dynamic recrystallization nucleation; A case study on ice[END_REF]. By local, we mean the scale of the grains that constitute crystalline materials such as ice, metals or most rocks. Dynamic recrystallization occurs in plastically deforming materials at high temperature. 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 to a drastic modification of the texture (crystallographic orientations of grains) [START_REF] Sakai | Dynamic and post-dynamic recrystallization under hot, cold and severe plastic deformation conditions[END_REF][START_REF] Kamb | Experimental recrystallization of ice under stress[END_REF][START_REF] Poirier | Creep of crystals[END_REF][START_REF] Humphreys | Recrystallization and related annealing phenomena[END_REF][START_REF] Doherty | Current issues in recrystallization: a review[END_REF]. These changes in the microstructure and texture may lead to modification of the rheological behavior, that is, the response of the material to imposed stresses or strains.

Ice is a hexagonal crystalline material (figure 1) that deforms, on Earth, at temperatures very close to its melting point. In the ductile regime, its deformation is characterized by a strong viscoplastic anisotropy, with dislocations gliding almost solely on the basal plane [START_REF] Glen | The Growth and Deformation of Ice Crystals[END_REF][START_REF] Duval | Rate controlling processes in the creep of polycrystalline ice[END_REF]. The orientation of an ice crystal relative to the direction of solicitation has therefore a strong control on its viscosity. Consequently, the distribution of crystallographic orientations in polycrystalline ice impacts the mechanical response in the form of a texture-induced viscoplastic anisotropy [START_REF] Castelnau | Anisotropic behavior of GRIP ices and flow in Central Greenland[END_REF].

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 the texture evolution along ice cores, a model has been proposed by De la Chapelle et al. 1998 [START_REF] De La Chapelle | Dynamic recrystallization and texture development in ice as revealed by the study of deep ice cores in Antarctica and Greenland[END_REF]. 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 predominate depends on the temperature and strain conditions along the ice core [START_REF] Thorsteinsson | Textures and fabrics in the GRIP ice core[END_REF][START_REF] Montagnat | Measurements and numerical simulation of fabric evolution along the Talos Dome ice core, Antarctica[END_REF]. 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, in particular when simple shear is dominant (see e.g. [START_REF] Samyn | Dynamic implications of discontinuous recrystallization in cold basal ice: Taylor Glacier, Antarctica[END_REF][START_REF] Smith | Ice fabric in an Antarctic ice stream interpreted from seismic anisotropy[END_REF][START_REF] Minchew | Processes controlling the downstream evolution of ice rheology in glacier shear margins: case study on Rutford Ice Stream, West Antarctica[END_REF][START_REF] Thomas | Microstructure and Crystallographic Preferred Orientations of an Azimuthally Oriented Ice Core from a Lateral Shear Margin: Priestley Glacier, Antarctica[END_REF]). 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 [START_REF] Pettit | The role of crystal fabric in flow near an ice divide[END_REF][START_REF] Gillet-Chaulet | Modélisation de l'écoulement de la glace polaire anisotrope et premières applications au forage de Dôme C[END_REF][START_REF] Martín | On the effects of anisotropic rheology on ice flow, internal structure, and the age-depth relationship at ice divides[END_REF][START_REF] Ma | Enhancement factors for grounded ice and ice shelves inferred from an anisotropic ice-flow model[END_REF][START_REF] Budd | Ice flow relations for stress and strain-rate components from combined shear and compression laboratory experiments[END_REF][START_REF] Mccormack | Modeling the Deformation Regime of Thwaites Glacier, West Antarctica, Using a Simple Flow Relation for Ice Anisotropy (ESTAR)[END_REF]. To model the deformation of an ice sheet, it is therefore essential to consider the effect of DRX on the texture evolution.

Dynamic recrystallization in ice has been extensively studied in the laboratory. These studies allowed to characterize the basic mechanisms at play, nucleation and grain-boundary migration, and the impact of strain heterogeneities on their kinetics (see e.g. [START_REF] Chauve | Strain field evolution during dynamic recrystallization nucleation; A case study on ice[END_REF][START_REF] Duval | Creep and recrystallization of polycrystalline ice[END_REF][START_REF] Duval | Creep and Fabrics of Polycrystalline Ice Under Shear and Compression[END_REF][START_REF] Bouchez | The fabric of polycrystalline ice deformed in simple shear : experiments in torsion, natural deformation and geometrical interpretation[END_REF][START_REF] Jacka | Ice crystallographic and strain rate changes with strain in compression and extension[END_REF][START_REF] Jacka | The steady-state crystal size of deforming ice[END_REF][START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF][START_REF] Qi | The down-stress transition from cluster to cone fabrics in experimentally deformed ice[END_REF][START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF]). [START_REF] Kamb | Experimental recrystallization of ice under stress[END_REF] [START_REF] Kamb | Experimental recrystallization of ice under stress[END_REF] was a pioneer in performing a large number of 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, characteristic of the tertiary creep stage could result from 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 [START_REF] Grennerat | Experimental characterization of the intragranular strain field in columnar ice during transient creep[END_REF] has revealed that the strain field during deformation in polycrystalline ice is highly heterogeneous, with no simple relation between intensity of the local strain and the Schmid factor, a parameter that quantifies the orientation relation 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 [START_REF] Chauve | Investigation of nucleation processes during dynamic recrystallization of ice using cryo-EBSD[END_REF] 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 with local microstructure changes.

Although this recent data show that DRX textures result from a complex interplay between the imposed stress or strain states, the local microstructure and the resulting redistribution of stresses, it also corroborates that development of DRX textures is systematically associated with the marked strain softening that signs the onset of tertiary creep [START_REF] Budd | Ice flow relations for stress and strain-rate components from combined shear and compression laboratory experiments[END_REF][START_REF] Duval | Creep and Fabrics of Polycrystalline Ice Under Shear and Compression[END_REF][START_REF] Jacka | Ice crystallographic and strain rate changes with strain in compression and extension[END_REF]. 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 on the initial texture and controlled by the imposed stress [START_REF] Jacka | The use of tertiary creep rates in ice at high strains in compression and shear[END_REF][START_REF] Jun | Horizontal shear rate of ice initially exhibiting vertical compression fabrics[END_REF][START_REF] Li | Deformation rates in combined compression and shear for ice which is initially isotropic and after the development of strong anisotropy[END_REF][START_REF] Treverrow | The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement[END_REF]. This contrasts with the minimum strain rate observed at the onset of secondary creep, which depends on the initial texture and that can be described by the classical Glen's law (a Norton-Hoff type law), which relates non-linearly the strain rate to the applied stress (through an exponent of 3 and an Arhenius type of dependence to the temperature) [START_REF] Duval | Rate controlling processes in the creep of polycrystalline ice[END_REF][START_REF] Glen | Experiments on the Deformation of Ice[END_REF]. 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 [START_REF] Budd | Ice flow relations for stress and strain-rate components from combined shear and compression laboratory experiments[END_REF][START_REF] Treverrow | The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement[END_REF] accounts for the impact of the DRX texture on the strain softening during tertiary creep.

Textures that develop in ice during DRX are well characterized. Under uniaxial compression, the textures measured at tertiary creep are characterized by c-axes oriented within a girdle at an angle between 30 • and 50 • to the compression axis depending on the experimental temperature T and stress. These textures are stable, that is, they evolve very slowly for strains higher than 10% [START_REF] Kamb | Experimental recrystallization of ice under stress[END_REF][START_REF] Duval | Creep and recrystallization of polycrystalline ice[END_REF][START_REF] Jacka | Ice crystallographic and strain rate changes with strain in compression and extension[END_REF][START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF][START_REF] Qi | The down-stress transition from cluster to cone fabrics in experimentally deformed ice[END_REF][START_REF] Treverrow | The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement[END_REF][START_REF] Piazolo | Dynamics of ice mass deformation: Linking processes to rheology, texture, and microstructure[END_REF]. During simple shear, tertiary creep textures seem to evolve with increasing strain from a two-maxima pattern, which is symmetric relative to the principal finite stretching direction, toward a strongly clustered texture, characterized by c-axes aligned in a direction perpendicular to the shear plane [START_REF] Duval | Creep and Fabrics of Polycrystalline Ice Under Shear and Compression[END_REF][START_REF] Bouchez | The fabric of polycrystalline ice deformed in simple shear : experiments in torsion, natural deformation and geometrical interpretation[END_REF][START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF][START_REF] Hudleston | Progressive development of fabrics across zones of shear in glacial ice[END_REF][START_REF] Qi | Crystallographic preferred orientations of ice deformed in direct-shear experiments at low temperatures[END_REF].

Texture development in ice deforming by dislocation creep has been successfully simulated using approaches of variable complexity [START_REF] Pettit | The role of crystal fabric in flow near an ice divide[END_REF][START_REF] Gillet-Chaulet | Modélisation de l'écoulement de la glace polaire anisotrope et premières applications au forage de Dôme C[END_REF][START_REF] Alley | Fabrics in polar ice sheets -Development and prediction[END_REF][START_REF] Castelnau | Modelling fabric development along the GRIP ice core, central Greenland[END_REF]. 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 [START_REF] Wenk | A deformation-based model for recrystallization of anisotropic materials[END_REF][START_REF] Thorsteinsson | An analytical approach to deformation of anisotropic ice-crystal aggregates[END_REF][START_REF] Llorens | Full-field predictions of ice dynamic recrystallisation under simple shear conditions[END_REF][START_REF] Richards | The evolution of ice fabrics: A continuum modelling approach validated against laboratory experiments[END_REF][START_REF] Rathmann | Effect of an orientation-dependent nonlinear grain fluidity on bulk directional enhancement factors[END_REF]. 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 [START_REF] Wenk | A deformation-based model for recrystallization of anisotropic materials[END_REF][START_REF] Thorsteinsson | An analytical approach to deformation of anisotropic ice-crystal aggregates[END_REF], or from limitations in the deformation modeling frame despite considering high complexity in the microstructure evolution [START_REF] Llorens | Full-field predictions of ice dynamic recrystallisation under simple shear conditions[END_REF]. In [START_REF] Richards | The evolution of ice fabrics: A continuum modelling approach validated against laboratory experiments[END_REF][START_REF] Rathmann | Effect of an orientation-dependent nonlinear grain fluidity on bulk directional enhancement factors[END_REF], the simulated textures are very close to the observed ones, but the kinetics of evolution is not. This discrepancy may result from the use of a Sachs approximation for modeling the mechanical response of the polycrystal, which results in a too soft mechanical response and incoherence in the grain rotation calculation [START_REF] Castelnau | Viscoplastic modeling of texture development in polycrystalline ice with a self-consistent approach : Comparison with bound estimates[END_REF].

In the present work, we propose a new formulation to predict the impact of DRX on the texture evolution of polycrystalline materials, like 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 the texture evolution is modeled by a physically-based orientation attractor that maximizes the resolved shear stress on the easiest slip system in the crystal (basal slip for ice). This set of equations, named CTI-RX, is solved in a full-field model using the finite-elements method (FEM).

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

CTI-RX-a model for the ice single crystal behavior

The Continuous Transverse Isotropic behavior

Following [START_REF] Meyssonnier | A model for the tangent viscous behaviour of anisotropic polar ice[END_REF][START_REF] Mansuy | Localization of deformation in polycrystalline ice: experiments and numerical simulations with a simple grain model[END_REF][START_REF] Gillet-Chaulet | A user-friendly anisotropic flow law for ice-sheet modelling[END_REF][START_REF] Burr | The anisotropic contact response of viscoplastic monocrystalline ice particles[END_REF][START_REF] Rathmann | Inferred basal friction and mass flux affected by crystal-orientation fabrics[END_REF], we account for the viscoplastic anisotropy of the ice single crystal by modelling its mechanical behavior using the continuous transverse isotropic (CTI) law. In this framework, the viscoplastic response of the ice single crystal depends only on the orientation of the c-axis, being isotopic in the basal plane (see Figure 1). This formulation is based on theoretical work by [START_REF] Boehler | On Experimental Testing Methods for Anisotropic Materials[END_REF] [START_REF] Boehler | On Experimental Testing Methods for Anisotropic Materials[END_REF], which has shown, through geometrical arguments, that the relation between the deviatoric stress S and the deviatoric strain rate D in a transverse isotropic symmetry may be described by a potential φ S :

D = ∂φ S ∂S D (1)
where φ S is expressed on the B C T I base, which is built out of a combination of the following four invariants:

B C T I = Tr (S 2 ), Tr (S 3 ), Tr (MS), Tr (MS 2 ) , (2)
with M = e 3 ⊗ e 3 where e 3 is the direction of the plane of isotropy (figure 1). The potential φ S is described below in the linear (φ (1) S ) and non linear (φ (n) S ) cases.

The linear CTI formulation

To obtain a linear CTI formulation, the potential, φ (1) S , must be quadratic. It is therefore constructed from the quadratic invariants taken from the B C T I base as:

φ (1) S = δ 1 Tr (S 2 ) + δ 2 Tr (MS) 2 + δ 3 Tr (MS 2 ). (3)
Considering the following derivation rules, the constitutive relation is given by

∂Tr (MS) ∂S = M, ∂Tr (S 2 ) ∂S = 2S, ∂Tr (MS 2 ) ∂S = MS + SM, (4)
D = ∂φ (1) S ∂S D = 2δ 1 S + 2δ 2 M D Tr (MS) + δ 3 (MS + SM) D . ( 5
)
The three independent parameters δ i in equation ( 5) can be fitted from 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 X indicates a tensor expressed in the grain reference frame:

• ψ 1 is the fluidity of the ice crystal subjected to shear on the basal plane:

g D 13 = ψ 1 2 g S 13 , g D 23 = ψ 1 2 g S 23 , (6)
• β is the ratio between the fluidity for shear parallel to the basal plane and that for shear normal to the basal plane:

g D 12 = β ψ 1 2 g S 12 , (7)
• γ is the ratio between the viscosity for compression (or traction) along the c-axis η c and for compression (or traction) in any direction g e r within the basal plane:

g S 33 = 2η c g D 33 , γ g S r r = 2η c g D r r ( 8
)
For solicitations of similar amplitude in different directions, for instance along e 3 or perpendicular to e 3 (i.e. along e r ) S = g S 33 = g S r r , in the linear case the strain rates are related by: g D r r = γ g D 33 [START_REF] Duval | Rate controlling processes in the creep of polycrystalline ice[END_REF] Using the above constraints on the single crystal behavior, identification gives (see appendix A for details):

δ 1 = ψ 1 β 4 , δ 2 = ψ 1 2 β γ + 2 4γ -1 -1 , δ 3 = ψ 1 2 1 -β ( 10
)
The constitutive equation (5) can also be expressed in an inverse manner as

S = 2α 1 D + 2α 2 M D Tr (MD) + α 3 (MD + DM) D (11)
and in terms of the potential

φ (1) D = α 1 Tr (D 2 ) + α 2 Tr (MD) 2 + α 3 Tr (MD 2 ) (12)
with : [START_REF] Gillet-Chaulet | Modélisation de l'écoulement de la glace polaire anisotrope et premières applications au forage de Dôme C[END_REF] and appendix A for details).

α 1 = η 1 β , α 2 = 2η 1 γ β -1 , α 3 = 2η 1 1 -1 β and η 1 = 1 ψ 1 (see

The non-linear CTI formulation

A non-linear potential φ (n) D linking S to D 1 n is obtained by taking the linear potential φ(1)

D = 1 2η 1 φ (1)
D to a higher order, k :

φ (n) D = 1 k φ(1) D k . ( 13
)
This hypothesis limits the number of independent parameters to 3 instead of 7 in the general case. A similar hypothesis was made in previous studies [START_REF] Gillet-Chaulet | Modélisation de l'écoulement de la glace polaire anisotrope et premières applications au forage de Dôme C[END_REF][START_REF] Rathmann | Effect of an orientation-dependent nonlinear grain fluidity on bulk directional enhancement factors[END_REF][START_REF] Gagliardini | Simulation numérique d'un écoulement bidimensionnel de glace polaire présentant une anisotropie induite évolutive[END_REF]. Using the chain rule for deriving φ n D gives:

S = ∂φ (n) D ∂D = φ(1) D k-1 ∂ φ(1) D ∂D ( 14
)
k is further obtained by identification of the power exponent on equation ( 14) in order to have

S ∼ D 1 n : 2(k -1) + 2 = 1 n =⇒ k = n + 1 2n (15)
Therefore, a non-linear CTI law can be written as:

φ (n) D = 4n n + 1 α 1 t r (D 2 ) + α 2 t r (MD) 2 + α 3 t r (MD 2 ) n+1 2n , (16)
leading to the following constitutive equation:

S = η ⋆ n 2α 1 D + 2α 2 M D Tr (MD) + α 3 (MD + DM) D , (17)
with the effective viscosity

η ⋆ n = 2η n α 1 t r (D 2 ) + α 2 t r (MD) 2 + α 3 t r (MD 2 )
1-n 2n [START_REF] Pettit | The role of crystal fabric in flow near an ice divide[END_REF] and α i parameters:

α 1 = 1 2β , α 2 = γ β -1 , α 3 = 1 - 1 β . ( 19
)
It is worth noting that β links g D 12 to g S 12 through:

g D 12 = ψ n 2 β n+1 2 g S n 12 . (20)
Therefore, β depends on the non-linearity, i.e. on the value of n. A similar enhancement will be obtained for a linear behavior n = 1 and a non-linear behavior with n = 3 by setting β n=3 = β n=1 .

Evolution of the crystal rotation in response to deformation and recrystallization

Deformation resulting from dislocation glide in the different slip planes results into a rotation of the crystal whose rate is described as in [START_REF] Gillet-Chaulet | Flow-induced anisotropy in polar ice and related ice-sheet flow modelling[END_REF][START_REF] Gagliardini | A review of anisotropic polar ice models: from crystal to ice-sheet flow models[END_REF], by:

∂c ∂t = Wc -λ Dc -c T Dc c , (21)
where W(u) = ∇(u)-∇ T (u) 2 and D(u) = ∇(u)+∇ T (u) 2 are the spin-and strain-rate tensors. According to equation [START_REF] Ma | Enhancement factors for grounded ice and ice shelves inferred from an anisotropic ice-flow model[END_REF], the rotation of the c-axis is due to the bulk spin, Wc, and to the viscoplastic spin expressed by λ Dcc T Dc c . λ is a parameter that depends on the plasticity model chosen for the ice crystal. In the present frame we will assume that deformation is only due to the contribution of dislocations gliding on the basal place, λ is set to 1. A value of λ < 1 is required to match a plasticity model that include prismatic and pyramidal slip [START_REF] Gillet-Chaulet | Flow-induced anisotropy in polar ice and related ice-sheet flow modelling[END_REF][START_REF] Gagliardini | A review of anisotropic polar ice models: from crystal to ice-sheet flow models[END_REF]. This assumption is well funded for ice and has been ofen made in previous studies [START_REF] Thorsteinsson | An analytical approach to deformation of anisotropic ice-crystal aggregates[END_REF][START_REF] Lliboutry | Anisotropic, transversely isotropic nonlinear viscosity of rock ice and rheolocical parameters inferred from homogenization[END_REF][START_REF] Azuma | A flow law for anisotropic polycrystalline ice under uniaxial compressive deformation[END_REF][START_REF] Svendsen | A continuum approach for modelling induced anisotropy in glaciers and ice sheets[END_REF][START_REF] Gödert | Induced anisotropy in large ice shields: Theory and its homogenization[END_REF].

To account for the impact of dynamic recrystallization on the evolution of the c-axis orientation, we add a new term to this equation based on the assumption that dynamic recrystallization induces a local re-orientation of the c-axis, either through nucleation of a grain with a new orientation, or by migration of a grain boundary, and that this re-orientation leads to local softening. This continuous description of the impact of DRX on the evolution of the c-axis orientation is supported by experimental observations of DRX nucleation in ice. These observations show that nucleation by bulging and polygonisation results in progressive misorientation of the recrystallized volumes relative to the parent grains, with misorientations most ofen < 20 • [START_REF] Chauve | Investigation of nucleation processes during dynamic recrystallization of ice using cryo-EBSD[END_REF].

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

∂c ∂t = Wc -λ Dc -c T Dc c + 1 Γ R X (c 0 -c) , (22)
where Γ R X is a parameter controlling the rate of rotation of the c-axis toward an attractor, c 0 , which represents the "ideal" orientation produced by recrystallization and whose formulation is described in section 2.3. The evolution of the c-axis orientation is therefore controlled by the balance between (i) the deformation-induced rotation and (ii) the rotation toward a recrystallization-driven attractor, c 0 . The meaning of Γ R X can be illustrated by considering the artificial case where only dynamic recrystallization contributes to the c-axis rotation. In this case, equation [START_REF] Budd | Ice flow relations for stress and strain-rate components from combined shear and compression laboratory experiments[END_REF] becomes

∂c ∂t = 1 Γ R X (c 0 -c) , (23)
The norm of ∂c ∂t can be interpreted as the rotation rate due to dynamic recrystallization only. This rate increases as ||c 0 -c||. Considering a constant c 0 , the solution for the c-axis rotation is:

c(t ) = (c(0) -c 0 ) e -t Γ R X + c(0). (24)
In this case, c rotates towards c 0 at a rate following an exponential decay with a characteristic time Γ R X , such that 95% (99%) of the rotation is obtained after 4Γ R X (5Γ R X ). In the CTI-RX model, c 0 is not constant but depends on the local deviatoric stress tensor, S.

Formulation of the recrystallization attractor c 0

The formulation of c 0 stems from the assumption that a grain (or part of a grain) rotates towards, or/and nucleates in, an orientation that facilitates the local viscoplastic deformation (dislocation glide). This assumption is supported by experimental observations which show that nucleation and grain boundary migration occur in areas of high strain incompatibility and produce orientations that tend to reduce this incompatibility [START_REF] Chauve | Strain field evolution during dynamic recrystallization nucleation; A case study on ice[END_REF][START_REF] Chauve | Investigation of nucleation processes during dynamic recrystallization of ice using cryo-EBSD[END_REF][START_REF] Chauve | Non-basal dislocations should be accounted for in simulating ice mass flow[END_REF]. Considering that the ice single crystal deforms mainly by dislocation glide on the basal plane, the recrystallizationinduced rotation should maximize the resolved shear stress (RSS) on the basal plane. The RSS for a dislocation with a Burgers vector a i (figure 1) gliding on the basal plane c in response to a deviatoric stress S is given by RSS(a i , c, S) = S : µ, [START_REF] Duval | Creep and Fabrics of Polycrystalline Ice Under Shear and Compression[END_REF] with µ = 1 2 c⊗a i +a i ⊗c ∥c∥.∥a i ∥ .

In the CTI formulation, the orientation for the maximum RSS is the orientation where the shear stress is maximum in the basal plane and the behavior is isotropic in this plane. Considering s 1 > s 2 > s 3 the eigenvalues and v i (i = 1 -3) the corresponding eigenvectors of the deviatoric stress tensor S, the maximum resolved shear stress is |s 1s 3 | and the associated orientations are [START_REF] Wu | Continuum mechanics and plasticity[END_REF]:

c 0 = 1 2 (v 1 ± v 3 ) (26)
Degenerated cases exit if at least two eigenvalues are equal:

• s 1 = s 2 , the possible attractors c 0 are the orientations at 45 • from v 3 .

• s 2 = s 3 , the possible attractors c 0 are the orientations at 45 • form v 1 .

• s 1 = s 2 = s 3 = 0, the possible attractors c 0 are all the orientations.

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

Interpretation of the orientation evolution equation

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

For uniaxial compression along the y-axis (figure 2a), the stable solution for the deformation term Wc -λ Dcc T Dc c corresponds to c aligned with y (red dot, figure 2a) [START_REF] Gillet-Chaulet | Modélisation de l'écoulement de la glace polaire anisotrope et premières applications au forage de Dôme C[END_REF]. For recrystallization 1 Γ R X (c 0c) , the stable solutions are the orientations at 45 • to the y-axis (blue line, figure 2a). When both deformation and recrystallization are active, the stable orientations of the c-axis are located on a girdle between these two end-members (black line, figure 2a). The position of this girdle depends on the value of Γ R X (black arrows, figure 2a). The higher is Γ R X the closer is the girdle to the compression axis.

For uniaxial tension along y, the stable solutions for the deformation term are all the orientations located in the xOz plane (red line, figure 2b) [START_REF] Gillet-Chaulet | Modélisation de l'écoulement de la glace polaire anisotrope et premières applications au forage de Dôme C[END_REF]. When both deformation and recrystallization are active, the stable orientations (black circle, figure 2b) form a girdle between 45 • and 90 • to the y-axis depending on the value of Γ R X .

For simple shear parallel to the xOz plane (with ∂u x ∂y < 0, figure 2c), there are four stable solutions of the recrystallization term that are ±y and ±x (figure 2c, blue dots), and two for the deformation term, ±y (figure 2c, red dot). The resulting stable orientations of the c-axis when both deformation and recrystallization are active are located at ±y and between x and x + y (and the symmetrical equivalent) (black dashed line dots, figure 2c). This latter stable orientation, whose exact position depends on the value of Γ R X , only exists for low Γ R X and vanishes for large Γ R X (Γ R X → +∞).

CTI-RX flow model description

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

(1) The momentum equation under the Stokes assumption:

-∇.S + ∇p = 0 (27)
(2) The non-linear CTI law:

S = η ⋆ n 2α 1 D + 2α 2 M D Tr (MD) + α 3 (MD + DM) D η ⋆ n = 2η n α 1 t r (D 2 ) + α 2 t r (MD) 2 + α 3 t r (MD 2 ) 1-n 2n (28)
(3) The equation for the evolution of the c-axis:

∂c ∂t = Wc -λ Dc -c T Dc c + 1 Γ R X (c 0 -c) , (29)
with

c 0 = v 1 ±v 3 2
where v i are the eigenvectors of S and |c.c 0 | is maximum.

Numerical scheme and implementation

The CTI-RX equations (section 2.5) are solved on a cuboid of polycrystalline ice deformed in creep (stress boundary conditions) in either uniaxial compression, uniaxial tension, or simple shear (Figure 3). A full-field approach using FEM is chosen to describe the field of orientations c at the element size. This implementation results in the R 3 iCe model (Rheology, Recrystallization, Rheolef, in Continous Transverse Isotropic material). This design allows direct comparisons of texture evolution and macroscopic mechanical behavior with data from laboratory creep experiments [START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF][START_REF] Qi | The down-stress transition from cluster to cone fabrics in experimentally deformed ice[END_REF][START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF].

. Uniaxial creep tests

In uniaxial creep, a zero vertical velocity (u y (bot t om) = 0) condition is applied at the bottom surface to simulate a fixed plateau, as in laboratory creep tests. To avoid translation of the sample, the position of one bottom summit is fully fixed (u(l e f t _bot t om_back) = (0, 0, 0)). To avoid rotation around the y-axis, a zero-velocity condition along the x-axis is also applied on a second summit (u x (l e f t _bot t om_ f r ont ) = 0). The sample is deformed in uniaxial compression (or tension) by applying a constant normal stress on the top surface (F y (t op) = 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 may deform faster. This may hinder the direct comparisons with laboratory experiments or other numerical models, in which homogeneous displacement conditions are imposed on the surface to which the compression is applied. This issue is addressed in section 4.1.

Simple shear creep tests

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

• u(bot t om) = 0, • u y = u z = 0 for the lateral surfaces (left, right, front, back) to avoid solid rotation, • u y (t op) = u z (t op) = 0.

Characteristic numbers and time

When adimensionalizied with respect to the material viscosity η n in Pa s 1 n , the vertical extent of the cuboid L in meter (Figure 3) and the applied macroscopic stress Σ (in Pa), the system of equations becomes:

-

∇S + ∇ p = 0 (30) S = η⋆ n 2α 1 D + 2α 2 M D Tr (M D) + α 3 (M D + DM) D η⋆ n = 2 α 1 Tr ( D2 ) + α 2 Tr (M D) 2 + α 3 Tr (M D2 ) 1-n 2n (31) ∂c ∂ t = W( ũ).c -λ D( ũ).c -(c T . D( ũ).c).c + M o(c 0 -c) (32)
where the superscript ( .) is used for all non-dimensional variables and operators. The dimensionless number that arises when scaling the c-axis evolution equation is M o:

M o = η n Σ n 1 Γ R X ( 33
)
M o is the ratio between a deformation characteristic time τ = η n Σ n and the recrystallization characteristic time Γ R X , which controls the rate of rotation of c toward the attractor c 0 . This dimensionless number quantifies the relative weight of the contributions of recrystallization versus deformation to the rotation of the c-axis.

Numerical algorithm

As the evolution of the c-axis orientation depends on the strain rate ( D( ũ)) and the effective viscosity η ⋆ depends on both the strain rate and the orientation of the c-axis, the system of coupled equations ( 31) and ( 32) is non-linear. To solve it, we linearize the equation describing the evolution of the c-axis using a second order backward differentiation formula (BDF2). The time discretization of the CTI-RX model results in the following system of equations:

-∇S t +1 + ∇ pt+1 = 0 (34) St+1 = η⋆ t +1 n 2α 1 Dt+1 + 2α 2 M t +1 D Tr (M t +1 Dt+1 ) + α 3 (M t +1 Dt+1 + Dt+1 M t +1 ) D η⋆ t +1 n = 2 α 1 Tr (( Dt+1 ) 2 ) + α 2 Tr (M t +1 Dt+1 ) 2 + α 3 Tr (M t +1 ( Dt+1 ) 2 ) 1-n 2n (35) 3 2∆ t c t +1 - 4 3 c t + 1 3 c t -1 = W( ũt+1 ).c t +1 -λ D( ũt+1 ).c t +1 -(c t +1 T . D( ũt+1 ).c t +1 ).c t +1 + M o(c 0 t +1 -c t +1 ) (36)
Finite elements and variational methods are used to solve the time-discretized problem on a Lagrangian grid within the C++ environment Rheolef [START_REF] Saramito | Efficient C++ finite element computing with Rheolef[END_REF]. The domain is spatially discretized using hexahedral finite elements (unless otherwise mentioned) and an initial c-axis orientation is prescribed for each element. The FE approximation for the tensorial mechanical parameter M and the scalar mechanical parameter c is computed with polynomes of order 0 (P 0 ). The FE approximation for the tensorial strain rate D, spin rate W and stress S fields is computed with polynomes of order 1 (P 1 ). The FE approximation for the displacement field is computed with polynomes of order 2 (P 2 ).

The Rheolef problem_mixed solver is used to solve the displacement field ũ and the pressure field p under the incompressibility hypothesis for the linearized version of equation [START_REF] Jun | Horizontal shear rate of ice initially exhibiting vertical compression fabrics[END_REF], which is obtained by fixing the apparent viscosity η ⋆ . A first fix loop is used to solve the non-linearity associated with η ⋆ . This loop is nested in a second fix loop that solves the temporal evolution of the c-axis orientation using equation [START_REF] Li | Deformation rates in combined compression and shear for ice which is initially isotropic and after the development of strong anisotropy[END_REF]. In practice, it is recommended for numerical efficiency to flatten the two nested fix loops. Algorithm 1 presents the numerical algorithm of two nested fix points implemented in R 3 iCe. This algorithm is a semi-implicit solver for the coupled equations [START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF] and [START_REF] Grennerat | Experimental characterization of the intragranular strain field in columnar ice during transient creep[END_REF]. The simulated orientation field c is thus consistent with the displacement field ũ and therefore with the deviatoric strain rate D and the deviatoric stress S. This implementation of CTI-RX model is distributed as the R 3 iCe code.

CTI parameters values for the ice single crystal

As stated in section 2, the values for the CTI parameters, γ, η n , n and β (equations ( 17) and ( 19)) are defined based on experimental data.

γ corresponds to the ratio of viscosities measured in compression and in tension for an ice single crystal that is well-oriented for basal slip. Laboratory experiments have shown that this ratio is very close to one. Therefore in the present numerical simulations, we set γ = 1, as previous studies [START_REF] Meyssonnier | A model for the tangent viscous behaviour of anisotropic polar ice[END_REF][START_REF] Mansuy | Localization of deformation in polycrystalline ice: experiments and numerical simulations with a simple grain model[END_REF].

Following [START_REF] Duval | Rate controlling processes in the creep of polycrystalline ice[END_REF] [START_REF] Duval | Rate controlling processes in the creep of polycrystalline ice[END_REF] and previous modeling work based on the CTI law [START_REF] Martín | On the effects of anisotropic rheology on ice flow, internal structure, and the age-depth relationship at ice divides[END_REF][START_REF] Rathmann | Effect of an orientation-dependent nonlinear grain fluidity on bulk directional enhancement factors[END_REF][START_REF] Mansuy | Localization of deformation in polycrystalline ice: experiments and numerical simulations with a simple grain model[END_REF], 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 [START_REF] Duval | Rate controlling processes in the creep of polycrystalline ice[END_REF]. Moreover, a value of 3 must be imposed in the CTI formulation in order to recover n = 3 for the polycrystal response.

β is linked to the ratio between basal and non-basal shearing viscosities. The strong viscoplastic anisotropy of ice is mostly carried by this parameter. In creep experiments at 1 MPa [START_REF] Duval | Rate controlling processes in the creep of polycrystalline ice[END_REF], the ratio between basal and non-basal shear viscosities is on the order of 10 4 . This ratio is related to β through equation [START_REF] Martín | On the effects of anisotropic rheology on ice flow, internal structure, and the age-depth relationship at ice divides[END_REF], β n+1 2 = 10 4 . For n = 3, β = 10 -2 is a value coherent with those considered in previous studies using the non-linear CTI formulation [START_REF] Martín | On the effects of anisotropic rheology on ice flow, internal structure, and the age-depth relationship at ice divides[END_REF][START_REF] Rathmann | Effect of an orientation-dependent nonlinear grain fluidity on bulk directional enhancement factors[END_REF][START_REF] Mansuy | Localization of deformation in polycrystalline ice: experiments and numerical simulations with a simple grain model[END_REF].

u t +1,k c =0,k nl =0 ← u t 7: c t +1,k c =0 ← c t 8:
while r es c > t ol c do ▷ Fix loop to solved (equation ( 36))

9:

k nl ← 0 10:
r es nl ← +∞ 11:

u t +1,k c ,k nl =0 ← u t +1,k c 12:
while r es nl > t ol nl do ▷ Fix loop to solve non linear CTI (equation ( 31))

13:

H ⋆ n ← η ⋆ n (u t +1,k c ,k nl , c t +1,k c ) ▷ equation (31)
14:

find u t +1,k c ,k nl +1 ← u t +1,k c ,k nl +1 (H ⋆ n , c t +1,k c
) ▷ problem_mixed solver on equation [START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF] 15:

r es nl ←||S(η ⋆ n (u t +1,k c ,k nl +1 , c t +1,k c ), u t +1,k c ,k nl +1 , c t +1,k c ) -S(H ⋆ n , u t +1,k c ,k nl +1 , c t +1,k c )|| 16: k nl ← k nl + 1 17:
end while 18:

H ⋆ n ← η ⋆ n (u t +1,k c ,k nl , c t +1,k c ) 19: compute S ← S(H ⋆ n , u t +1,k c ,k nl ) ▷ equation (31) 20:
find c 0 ← c 0 (S) ▷ equation ( 26)

21: compute c t +1,k c +1 ← c t +1,k c +1 (c 0 , u t +1,k c ,k nl ) ▷ equation (36) 22: r es c ←||S(η ⋆ n (u t +1,k c ,k nl , c t +1,k c +1 ), u t +1,k c ,k nl +1 , c t +1,k c +1 ) -S(H ⋆ n , u t +1,k c ,k nl , c t +1,k c )|| 23: k c ← k c + 1 24:
end while 25: ) MPa s 1/3 for the polycrystal are obtained. The value of 7.5 MPa s 1/3 is taken for η 3 in order to reproduce the experiment duration as observed during uniaxial compression tests that are taken, here, as a reference (see section 5.2). The parameters of the CTI law for the ice single crystal are summarized in table 2.

u t +1 ← u t +1,k c ,k nl

Validation of the CTI-based implementation for the prediction of the viscoplastic response

In this section, we verify that the CTI law implemented in R 3 iCe correctly reproduces the mechanical response of an ice polycrystal deforming by dislocation creep only (no DRX). This is a prerequisite before improving the CTI formulation with a representation of DRX. The ice mechanical response of the polycrystalline ice is known to be characterized by strong strain and stress heterogeneities, which are very likely precursors of dynamic recrystallization [START_REF] Grennerat | Experimental characterization of the intragranular strain field in columnar ice during transient creep[END_REF]. We first compare the R 3 iCe predictions to that of the CraFT-EVP model [START_REF] Grennerat | Experimental characterization of the intragranular strain field in columnar ice during transient creep[END_REF], which has been shown to adequately reproduce these heterogeneous strain fields observed in experiments [START_REF] Grennerat | Experimental characterization of the intragranular strain field in columnar ice during transient creep[END_REF][START_REF] Piazolo | Effect of local stress heterogeneities on dislocation fields: Examples from transient creep in polycrystalline ice[END_REF]. Then we 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 way.

Stress and strain field predictions: comparison with CraFT-EVP

CraFT-EVP has been used to provide full-field predictions of the mechanical response of ice during primary and secondary creep [START_REF] Grennerat | Experimental characterization of the intragranular strain field in columnar ice during transient creep[END_REF]. It is based on an elasto-viscoplastic formulation proposed by Suquet et al. 2012 [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF] and solved using the CraFT software. The predictions of CraFT-EVP for the mechanical response as well as strain fields were validated by direct comparisons with mechanical tests on ice polycrystals that included strain fields measurements by Digital Image Correlation [START_REF] Grennerat | Experimental characterization of the intragranular strain field in columnar ice during transient creep[END_REF][START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF]. The macroscopic strain fields predicted by CraFT-EVP were in good agreement with the experimental observations at 1% strain (secondary creep, prior to any texture evolution owing to DRX). With no texture evolution, CraFT-EVP and R 3 iCe both simulate a stationary secondary creep.

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

Figure 4(a-b) show the polycrystal microstructure generated by a grain growth algorithm using Neper [START_REF] Quey | Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing[END_REF] used for both CraFT-EVP and R 3 iCe simulations. In this configuration, grains are defined as a group of connected elements that share the same initial orientation. CraFT-EVP requires a regular cubic mesh, while the R 3 iCe polycrystalline microstructure used here was meshed with a tetrahedral geometry. The resolution of CraFT-EVP simulations were reduced to match the R 3 iCe one. Figure 5 shows the component ε y y of the strain field and the component σ y y of the stress field on a x y section at z = 0, for a uniaxial compression creep test and after 1% of macroscopic strain. The simulated strain fields do not compare perfectly between CraFT-EVP and R 3 iCe predictions although the overall patterns do not show strong mismatch. Location of areas in compression (blue) and in extension (red) are similarly prescribed (Figure 5a). Significant differences between the responses of the models CraFT-EVP and R 3 iCe are restricted to the first layers of grains close to the top surface. These differences result from the different boundary conditions applied in the two models. In R 3 iCe, the Neumann boundary condition on the top 14Thomas Chauve, Maurine Montagnat, Véronique Dansereau, Pierre Saramito, Kévin Fourteau and Andréa Tommasi surface results on a constant force applied on each element that can lead to heterogeneous displacement on this surface. As a result, some "weaker" grains deform faster, as represented by the dark blue areas on Figure 5a. In CraFT-EVP, periodic boundary conditions are imposed owing to the fact that equations are solved with a Fast Fourier Transform-based method. The simulated stress fields (Figure 5b) compare satisfactorily in terms of amplitude and global pattern although the precise stress localisation shows more scattering. Some of this mismatch may also be explained by the differences in boundary conditions between the two models.

Similarly, the distributions of predicted Von Mises equivalent stress and strain by the two models show some discrepancies (figure 6). These discrepancies are stronger between the two configuration in which grains are described by multiple elements (figure 6) and where the impact of boundary conditions just mentioned is the stronger. Distributions for all components of both tensors ε i j , σ i j are given in appendix B.

In R 3 iCe, the stress tensor at the element scale is used to predict the recrystallization attractor c 0 . This latter is based on the principal direction of the stress tensor v 1 and v 3 (section 2.3). Figure 7 shows the distribution of v 1 and v 3 orientations predicted using both CraFT-EVP (figure 7a), and R 3 iCe with a configuration with several elements per grain (figure 7b) and for one orientation per element (figure 7c). All three configurations predict similar v 1 and v 3 orientation distributions and therefore similar c 0 orientation distributions are expected.

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

Stress and strain field prediction in a simplified mesh configuration

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

The distributions of the Von Mises equivalent stress and strain in this one-orientation-permesh-element simulations are compared to those of the simulations with meshes at sub-grain scale using both the CraFT-EVP and R 3 iCe models (figure 6). The distributions obtained in this simplified microstructure representation compare well with those predicted using the CraFT-EVP model, here considered as a reference. The fact that the CTI-RX fields predictions obtained with a sub-grain scale mesh departs slightly from the Craft-EVP predictions is attributed, as mentioned before, to the discrepancy in 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 case of the simplified microstructure.

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

Texture evolution: comparison with experimental data

The aim of this section is to test the ability of the R 3 iCe model to predict the evolution of the c-axes orientation with and without dynamic recrystallization (equation ( 32)). As mentioned in the introduction, deformation experiments with no or very limited DRX are not possible due to the fast strain rates and high temperature conditions required. Deformation-only textures may 4 a-b). The full line shows a R 3 iCe simulation where a single orientation is attributed per element (figure 4 c) only be sampled in deep ice cores from the central parts of Greenland or Antarctica [START_REF] De La Chapelle | Dynamic recrystallization and texture development in ice as revealed by the study of deep ice cores in Antarctica and Greenland[END_REF]. These textures were well reproduced by uniaxial compression simulations performed using the meanfield ViscoPlastic Self Consistent model [START_REF] Montagnat | Measurements and numerical simulation of fabric evolution along the Talos Dome ice core, Antarctica[END_REF][START_REF] Castelnau | Modelling fabric development along the GRIP ice core, central Greenland[END_REF]. R 3 iCe prediction for texture evolution without recrystallization are therefore compared to those of VPSC simulations extracted from [START_REF] Montagnat | Measurements and numerical simulation of fabric evolution along the Talos Dome ice core, Antarctica[END_REF]. In contrast, there is a large experimental dataset in which texture evolved in response to both deformation and DRX. We use this dataset to constrain the predictions of R 3 iCe for the effect of DRX on the texture evolution and mechanical behavior in uniaxial compression, uniaxial tension and simple shear.

The laboratory creep experiments selected to validate the model are:

• where Σ stands for the macroscopic applied stress and Σ 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 [START_REF] Qi | The down-stress transition from cluster to cone fabrics in experimentally deformed ice[END_REF]. These conditions enabled to reach higher flow stresses without failure (Σ ∼ (1.3, 4.3) MPa, T=-10 • C ).

For each configuration, simulations using a configuration with one independent orientation

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.

Texture evolution in absence of recrystallization

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

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

A strong cluster (single maximum of c-axes) texture develops as strain increases. This texture evolution induces geometrical hardening of the polycrystal, expressed by a decrease of the strainrate by a factor ∼ 1.5 between 20% and 60% total strain, when the texture is highly clustered (a 1 ≈ 0.71).

Similar cluster textures are observed along deep cold ice cores, for instance along the Talos Dome ice core [START_REF] Montagnat | Measurements and numerical simulation of fabric evolution along the Talos Dome ice core, Antarctica[END_REF] where a clustered texture with a 1 ≈ 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. [START_REF] Montagnat | Measurements and numerical simulation of fabric evolution along the Talos Dome ice core, Antarctica[END_REF] using the VPSC model, when departing from an isotropic texture (see red dots in figure 8a).

Uniaxial compression

The evolution with strain (i.e. time) of the eigenvalues (a i , i = 1, 3) of the second order orientation tensor 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 figure 9. It is compared with experimental measurements from Montagnat et al. 2015 [START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF]. These latter were obtained from initially isotropic granular ice compressed at -5 • C under a constant load with initial stresses of 0.7 and 0.8 MPa, and up to various strains (from 2 to 17.8% bulk shortening).

Comparison of the c-axis pole figures shows a good qualitative agreement between the simulated and the experimental textures (figure 9(d-e)). Both are girdle textures that become more concentrated and show an evolution of the cone angle to the compression axis from 45 • to 35 • up to ∼ 20% macroscopic strain and remain nearly constant thereafter (figure 9(d)). There is also a good agreement between the predicted evolution of the eigenvalues of the texture and those measured at 7%, 12% and 17.8% (figure 9(a)). This good agreement between experimental and simulated textures was obtained for a value of Γ R X = 13.8 days. During the experiments, 17.8% strain was reached in 144 hours (∼ 6 days). In the model simulations, this strain is reached in a ). Red dots are predictions from the VPSC model taken from [START_REF] Montagnat | Measurements and numerical simulation of fabric evolution along the Talos Dome ice core, Antarctica[END_REF]; b) evolution with strain of the macroscopic strain rate; c) evolution of the intensity of the texture with the colatitude θ over a section defined on the Z plane (dashed blue lines in d)). Data at 22%, 40% and 60% macro strain are shown. The white dashed lines define the angle at ±45 • from the compression axis Y ; d) textures presented as pole figures for 22% (i ), 40% (i i ), and 60% (i i i ) 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).

shorter time (∼ 130 hours) but the tertiary creep (maximum strain-rate) regime is reached between 10 to 20% strain, similarly to the experimental observations (see also [START_REF] Jacka | Ice crystallographic and strain rate changes with strain in compression and extension[END_REF][START_REF] Treverrow | The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement[END_REF]).

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

The bulk strain associated with the maximum strain rate and the simulated enhancement factor cannot be directly compared to that measured in Montagnat et al. 2015 [29] because these experiments were performed under constant load and the area of the sample section evolved during the tests, whereas the model simulations are performed in a constant stress configuration. The simulated strain rate may be corrected by assuming a constant sample volume during the experiment, what leads to εcor r si m = εsim × (1 + ε) 3 . This corrected strain rate is plotted in blue in figure 9(b). It brings the maximum strain rate to ∼ 10% and the enhancement factor to 3.0. These values are closer to the macroscopic response measured by Montagnat et al. [START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF], which is characterized by an enhancement factor between 3.3 to 3.8. 3 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 [START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF]; b) simulated macroscopic strain rate (black) and corrected simulated macroscopic strain rate (blue) to match a constant force boundary conditions instead of the constant stress one used in the simulation. EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude θ over a section defined on the Z plane (dashed blue lines in (d)). Data at 7%, 12% and 17.8% macro strain are shown. The white dashed lines define the angle at ±45 • from the compression axis Y ; d) simulated textures (pole figures) for 7% (i ), 12% (i i ), and 17.8% (i i i ) macroscopic strain; e) measured textures for the same macroscopic strains [START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF]. All figures presenting texture data are color-coded using the same colorbar (scales in multiples of a uniform distribution).

Uniaxial tension

Experimental creep tests performed in uniaxial tension are scarce. Jacka et al. 1984 [27] presented one experiment performed at -3 • C, up to 9.3% octahedral shear strain that corresponds, following the definition given in Jacka et al. 1984 [START_REF] Jacka | The time and strain required for development of minimum strain rates in ice[END_REF], to 13.1% axial strain. At this strain level, they did not attain the quasi-constant strain rate typical of tertiary creep. The texture obtained, measured by the classical manual Rigsby stage technique with one orientation per grain, is presented in figure 10. The texture is characterized by a small circle girdle with a mean half angle of 50.4 • and a standard deviation of 15.3 • . Although tertiary creep was not reached in this experiment, an enhancement factor of ∼ 3, similar to that observed in compression in the same conditions, is measured.

The texture evolution and macroscopic response predicted by a simulation under uniaxial tension creep are shown in figure 11. The texture evolves with increasing strain towards a girdle texture with a cone angle slightly around 45 • to the extension direction. This texture seems to stabilize at ∼ 5% axial strain. The distribution of the c-axis angle from the tension axis is given in figure 12 for the simulated texture at 13% strain. It is similar to the one measured by Jacka et al. 1984 [START_REF] Jacka | Ice crystallographic and strain rate changes with strain in compression and extension[END_REF] (figure 10). An enhancement factor of about 10 is simulated at 10% axial strain.

Simple shear

The simulated texture evolution in simple shear is compared with observations made by Journaux et al. [START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF] for torsion experiments on isotropic granular ice performed at -7 • C, under a constant torque corresponding to a maximum shear stress between 0.4 and 0.6 MPa.

The experimental textures are characterized by two sub-maxima, one almost perpendicular to the shear plane, called M 1 and the second one, M 2, initially at low angle to the shear plane. Similar evolution was also observed in simple shear in [START_REF] Bouchez | The fabric of polycrystalline ice deformed in simple shear : experiments in torsion, natural deformation and geometrical interpretation[END_REF][START_REF] Qi | Crystallographic preferred orientations of ice deformed in direct-shear experiments at low temperatures[END_REF][START_REF] Burg | Dynamic recrystallization and fabric development during the simple shear deformation of ice[END_REF]. As the shear strain increases, M 2 merges with M 1 to form a highly concentrated cluster texture characterized by a single maximum of c-axes normal to the shear plane. The simple shear simulations performed here reproduce this evolution relatively well, as shown in figure 13.

By accounting for dynamic recrystallization, the R 3 iCe model is able to reproduce the evolution of both M 1 and M 2 maxima with strain. At low shear strains (γ ∼ 0.2), M 1 concentrates normal to the shear plane while M 2 follows the principal extension direction (white dashed line, figure 13(c) bottom and figure 13(d)). As strain increases, M 2 rotates towards M 1. At a shear strain of 2, the texture is characterized by a slightly elongated single cluster normal to the shear plane as M 1 and M 2 merge. The evolution of the eigenvalues of the second order orientation tensor, as observed in the experiments are also well reproduced (figure 13(a)).

In both experiments and simulations, the strain rate increases rapidly up to γ = 0.3 while the double-maxima texture develops (figure 13(d)). It continues to increase, but at a progressively slowing rate as the texture evolves towards a single maximum normal to the shear plane, up to 3 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 θ over a section defined on the plane normal to Z (dashed blue lines in d)). Data at 5%, 13% and 15% axial macroscopic strain are shown. The white dashed lines define the angle at ±45 • from the compression axis Y ; d) simulated textures (represented as pole figures) for 5% (i ), 13% (i i ), and 15% (i i i ) axial macro strain. All figures presenting texture data are color-coded using the same colorbar (scales in multiples of a uniform distribution).

γ ∼ 3 where is reaches a quasi-steady state. As explained in section 5.2, the discontinuities in the strain-rate curve (figure 13(b)) are technical artefacts imposed by the operating rules of GRICAD clusters.

In phase with Treverrow et al. 2012 [START_REF] Treverrow | The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement[END_REF], 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 (E F ) of 6.7 that matches Treverrow et al. 2012 [START_REF] Treverrow | The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement[END_REF] data for a similar applied shear stress.

Testing the relative effects of deformation and recrystallization on texture evolution under higher stresses

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

To test the sensitivity of the dynamic recrystallization formulation to the level of stress, we performed simulations under creep conditions with the stress values from Qi et al. 2017 [START_REF] Qi | The down-stress transition from cluster to cone fabrics in experimentally deformed ice[END_REF] experiments. As CTI-RX equations are relating the deviatoric stress to the deviatoric strain and are solved under 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 [START_REF] Mellor | Deformation and failure of ice under constant stress or constant strain-rate[END_REF]. Stable textures expected at both stages can therefore be compared.

The comparison between simulated and observed textures at different compressive flow stresses is shown in figure 14. The experimental observations are well reproduced. Both measured and simulated textures show a transition from a girdle texture to a weak single maximum with increasing flow stress for a macroscopic strain of ∼ 20%. To correctly fit the experimental data, a value of Γ R X of 4.6 days is used, which is lower than that used to reproduce the unconfined compression creep tests of Montagnat et al. 2015 [START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF] (table 3). This implies a higher contribution of recrystallization compared to deformation at higher stresses.

Discussion

The CTI-RX model presented here enables the simulation of texture evolution in polycrystalline ice, as shown in the simulations in compression, simple shear and tension, including or not dynamic recrystallization. Without recrystallization, the texture evolution predicted by the model reproduces correctly the textures observed along deep ice cores where dynamic recrystallization is supposed to have a low impact. 3 and comparison with experimental data from [START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF]. 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 Journaux et al. 2019 [START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF]; b) simulated macroscopic strain rate. EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude θ over a section defined on the Z plane (dashed blue lines in d)). Data at γ = 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 γ = 0.2 (i ), γ = 0.42 (i i ) and γ = 1.96 (i i i ); e) measured textures for the same macroscopic shear strains [START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF]. All figures presenting texture data are color-coded using the same colorbar (scales in multiples of a uniform distribution).

The novelty of the present model, and its implementation, stands on the simulation of the evolution of the texture (orientations of c-axes) based on a new formulation that (i) introduces an attractor controlled by the local stress field to reproduce the impact of DRX mechanisms (ii) is coupled with the continuous transverse isotropic (CTI) law to account for strain and stress field heterogeneities at the crystal scale arising from the crystal viscoplastic anisotropy, (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 derive 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 [START_REF] Chauve | Strain field evolution during dynamic recrystallization nucleation; A case study on ice[END_REF][START_REF] Grennerat | Experimental characterization of the intragranular strain field in columnar ice during transient creep[END_REF][START_REF] Piazolo | Effect of local stress heterogeneities on dislocation fields: Examples from transient creep in polycrystalline ice[END_REF] 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) [START_REF] Chauve | Strain field evolution during dynamic recrystallization nucleation; A case study on ice[END_REF][START_REF] Chauve | Investigation of nucleation processes during dynamic recrystallization of ice using cryo-EBSD[END_REF][START_REF] Chauve | Non-basal dislocations should be accounted for in simulating ice mass flow[END_REF][START_REF] Piazolo | Effect of local stress heterogeneities on dislocation fields: Examples from transient creep in polycrystalline ice[END_REF][START_REF] Chauve | Strain field evolution at the ductile-to-brittle transition: a case study on ice[END_REF]. 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 grains orientations. The last important experimental observation is that the development of recrystallization textures leads to the weakening of the polycrystal mechanical response. All these observations justify the choice of a continuous formulation for the c-axis orientation evolution (equation ( 22)), with a targeted recrystallized orientation that is calculated so that to maximize the local basal resolved shear stress (equation ( 26)).

Comparison with tertiary creep experiments made in compression [START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF][START_REF] Treverrow | The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement[END_REF][START_REF] Piazolo | Dynamics of ice mass deformation: Linking processes to rheology, texture, and microstructure[END_REF], in simple shear [START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF][START_REF] Treverrow | The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement[END_REF][START_REF] Qi | Crystallographic preferred orientations of ice deformed in direct-shear experiments at low temperatures[END_REF] and in tension [START_REF] Jacka | Ice crystallographic and strain rate changes with strain in compression and extension[END_REF] supports the fact that the R 3 iCe model is able to reproduce texture evolution in conditions similar to those of laboratory experiments (T = -5 to -10 • C, S ∼ 1 MPa).

For that, a single parameter is adjusted: the recrystallization parameter Γ R X , which weights the relative contribution of deformation and recrystallization to the rotation of the c-axes. As expected, at e.g. similar stress conditions, a lower temperature requires a higher Γ R X value (and therefore a lower weight of dynamic recrystallization) to correctly reproduce the experimental texture evolution. A noticeable result is that the two-maxima texture measured in simple shear (in natural and laboratory ice) is correctly reproduced here, together with the kinetics of its evolution with strain toward a single-maximum cluster.

For all conditions simulated with R 3 iCe, we show that the CTI-RX model reproduces accurately the softening at the transition from secondary to tertiary creep. The kinetics of the transition and the resulting enhancement factors, when departing from an isotropic texture, are correctly simulated as well, if compared with experimental results performed in similar conditions [START_REF] Budd | Ice flow relations for stress and strain-rate components from combined shear and compression laboratory experiments[END_REF][START_REF] Jacka | Ice crystallographic and strain rate changes with strain in compression and extension[END_REF][START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF][START_REF] Li | Deformation rates in combined compression and shear for ice which is initially isotropic and after the development of strong anisotropy[END_REF][START_REF] Treverrow | The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement[END_REF]. Since 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, the softening associated with dislocations annihilation and grain boundary migration during dynamic recrystallization seems to play a secondary role.

Simulations at higher uni-axial compressive stresses (from 1 to 4 MPa) have been performed to compare with the experiments of Qi et al. 2017 [START_REF] Qi | The down-stress transition from cluster to cone fabrics in experimentally deformed ice[END_REF]. Contrary to previously mentioned experiments, these latter 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 from 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 the texture evolution. R 3 iCe run under imposed stress without confining pressure is able to correctly represent this texture clustering with stress. Indeed, experimental textures measured at 20% strain are satisfyingly reproduced by R 3 iCe providing that Γ R X is given a relatively low value (4.6 days for T = -10 • C ) compared to the one determined for unconfined creep compression experiments (13.8 days for T = -5 • C). This low value of Γ R X indicates that, in the simulations, a high amount of recrystallization is necessary to reproduce the experimental textures and to counterbalance the impact of deformation on the c-axes rotation. We interpret this apparent contradiction as due to activation of brittle-ductile processes in the experiments. This assumption is based on observations by Kalifa et al. [START_REF] Kalifa | Microcracking and the failure of polycrystalline ice under triaxial compression[END_REF], which performed triaxial tests under the same range of strain rate (∼ 10 -5 -10 4 s -1 ) than Qi et al. 2017 [START_REF] Qi | The down-stress transition from cluster to cone fabrics in experimentally deformed ice[END_REF], but under variable confining pressure from 0 to 10 MPa. By careful observations made right 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. Micro-cracking is therefore also expected to have occurred in Qi et al. 2017 [START_REF] Qi | The down-stress transition from cluster to cone fabrics in experimentally deformed ice[END_REF] experiments. Strong stress concentration happens at crack tips and Chauve et al. 2017 [START_REF] Chauve | Strain field evolution at the ductile-to-brittle transition: a case study on ice[END_REF] have recently shown that, for ice, this stress concentration can be released through DRX. The R 3 iCe formulation does not account for neither micro-cracking nor an effect of the confining pressure. The enhancement of recrystallization due to plastic energy available around crack tips needs therefore to be reproduced by increasing the recrystallization rate (Γ R X ).

Recent models of texture evolution during dynamic recrystallization have been proposed by Rathmann et al. 2021 [START_REF] Rathmann | Inferred basal friction and mass flux affected by crystal-orientation fabrics[END_REF] and Richard et al. 2021,2022 [START_REF] Richards | The evolution of ice fabrics: A continuum modelling approach validated against laboratory experiments[END_REF][START_REF] Richards | Ice fabrics in two-dimensional flows: beyond pure and simple shear[END_REF], aiming to account for textureinduced anisotropy of ice in large-scale models. Their formulation is based on an orientation distribution function (ODF) evolution equation in which dynamic recrystallization is represented by two terms: one term that accounts for the production of orientations based on a given parameter, the deformability, that is calculated using the Sachs constant stress assumption, and a second term that accounts for the diffusion of orientations. This approach is not intended to predict the mechanical response evolution due to texture development, but the texture evolution equation can be compared to that in the CTI-RX formulation. The deformability term can be obtained from equation ( 22) by assuming homogeneous stresses (Sachs model). The diffusive term, which is necessary to reproduce the texture dispersion effect owing to recrystallization is spontaneously obtained in full-field approaches such as R 3 iCe, due to the heterogeneous stress field. Although these models accurately predict steady-state textures for different strain geometries, the predictions for transient textures evolution are less accurate than that from R 3 iCe. This highlights the fundamental role of stress heterogeneities in texture evolution during the transition between secondary and tertiary creep (from 1 to ∼ 15% strain).

The major strength of [START_REF] Richards | The evolution of ice fabrics: A continuum modelling approach validated against laboratory experiments[END_REF][START_REF] Rathmann | Inferred basal friction and mass flux affected by crystal-orientation fabrics[END_REF][START_REF] Richards | Ice fabrics in two-dimensional flows: beyond pure and simple shear[END_REF] formulations is the numerical efficiency that makes it suitable for large-scale modelling (glaciers, ice sheets), providing that transient behaviors are not of concern. In such context, the R 3 iCe model could be a good candidate to provide constraints to formulate parameterizations for the temporal evolution, destined to less costly formulations.

The CTI-RX model is not limited to the loading conditions prescribed here or to polycrystalline ice. Its formulation is generic and can be extended to other polycrystalline materials, provided that their single crystal viscoplastic behavior can be adequately described using the continuous transverse isotropic law as defined in equation [START_REF] Thomas | Microstructure and Crystallographic Preferred Orientations of an Azimuthally Oriented Ice Core from a Lateral Shear Margin: Priestley Glacier, Antarctica[END_REF]. For instance, the versatility of the CTI-RX model could enable analyzing the mechanical behavior of materials such as e.g. magnesium and quartz, which also show a strong viscoplastic anisotropy with hexagonal symmetry.

Conclusions

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

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

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

The accurate reproduction of the textures in the R 3 iCe model leads to a good prediction of the mechanical softening associated with dynamic recrystallization. It 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.

In conditions of high stresses, where confining pressure is requested in the laboratory to prevent failure, we propose that the local fracturing at the crystal scale likely enhances dynamic recrystallization processes and fasten texture evolution with strain. This is consistent with the lower Γ R X required to fit the texture measured during the confined compression experiments [START_REF] Qi | The down-stress transition from cluster to cone fabrics in experimentally deformed ice[END_REF] relative to that used to fit the low stress unconfined compression creep experiments [START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF], despite the higher experimental temperature of the latter (-5 • C vs. -10 • C , cf. table 3). These results are an additional illustration of the ability of the R 3 iCe model to help resolve open questions about ice deformation behavior.

The computational cost of the R 3 iCe full-field model is not adapted to a direct implementation in large-scale flow modeling frameworks. The model is nonetheless highly valuable to constrain the parameterization required to properly account for the texture evolution and its impact on the mechanical response in complex or changing boundary conditions as those prevailing in the bottom of deep ice cores, ice streams, and glaciers. The constitutive equation gives :

g D 12 = 2δ 1 g S 12 (44)
Therefore :

δ 1 = ψ 1 β 4 ( 45
)
(2) δ 3 is given by applying a stress g S = The constitutive equation gives :

g D 13 = (2δ 1 + δ 3 ) g S 13 [START_REF] Richards | The evolution of ice fabrics: A continuum modelling approach validated against laboratory experiments[END_REF] Therefore

δ 3 = ψ 1 2 1 -β (48)
(3) δ 2 identification is less trivial and presented in more details.

• First, a uniaxial stress σ is applied along the c-axis. The stress tensor is g σ = • Finally from the grain behavior we can write that γ g D 33 = g D 11 : • This gives :

γ 4 3 δ 1 + 8
δ 2 = 6δ 1 1 -γ 4γ -1 -δ 3
• and using δ 1 and δ 3 as defined above:

δ 2 = ψ 1 2 β γ + 2 4γ -1 -1

Appendix B. Strain and stress field comparison

The figures 15 and 16 show the probability distribution functions for all components of the strain and stress tensors during uniaxial creep with c-axis evolution (section 4.1).

1 d

 1 i m(A) Tr (A) A D i j = A i j -δ i j d i m(A) Tr (A) ⊗ Tensorial (dyadic) product A = b ⊗ c A i j = b i c j

e 3 Figure 1 .

 31 Figure 1.The transverse isotropic geometry is fully described by one direction, e 3 , that defines the plane of isotropy. For an ice single crystal, this direction corresponds to that of the c-axis (i.e. the normal to the basal (0001) plane).

Figure 2 .

 2 Figure 2. Pole figure schematic representation of the stable orientations for the evolution equation (eq. 22) for various stress configurations S: (a) uniaxial compression along y, (b) uniaxial tension along y and (c) simple shear in the xOz plane, ∂u x ∂y < 0 (upper hemisphere).

Figure 3 . 1 .

 31 Figure 3. Schematic representation of the rectangular cubic simulation domain, with surface and point boundaries labeled.

  12Thomas Chauve, Maurine Montagnat, Véronique Dansereau, Pierre Saramito, Kévin Fourteau and Andréa Tommasi 1: Inputs: c i ni t 2: Initialize: c -1 ← c i ni t c 0 ← c i ni t u 0 ← 0 3: for t from 0 to N do ▷ Time evolution loop 4:

26 : 3 Table 2 .

 2632 c t +1 ← c t +1,k c 27: end for Algorithm 1: Numerical algorithm, based on two nested fix loops: one for solving the viscosity, η ⋆ n , nested in a second one solving the temporal evolution of the orientation of c-axis. -2 7.5 M P a.s 1/Parameters of the CTI law for the ice single crystal. The viscosity η n can be adjusted by comparison with either the experimental response of the single crystal or that of isotropic polycrystals. Solving the CTI law for uniaxial compression creep of a single crystal oriented at 45 • from the compression axis and of a polycrystal with random orientations picked from a uniform texture results in macroscopic strain rate of Dsim sc = 1.7×10 -2 η n Σ n for the single crystal and Dsim px = 5.3×10 -5 η n Σ n for the polycrystal. The viscosity η n is obtained by comparing the simulated values to the experimental data at 1 MPa [9] (D exp sc = 2 × 10 -5 , 2 × 10 -4 s -1 for the single crystal, D exp px = 9 × 10 -8 , 3 × 10 -7 s -1 for the polycrystal). η sc 3 = (4.4, 9.4) MPa s 1/3 for the single crystal and η px 3 = (5.6, 8.4

Figure 4 .

 4 Figure 4. Microstructure generated via a grain growth using Neper meshed for (a) CraFT-EVP microstructure with 265302 grid points (b) R 3 iCe with 262688 tetrahedral elements. (c) Simplified microstructure in which orientations are defined at the element scale with hexahedral mesh of 2740 elements.

Figure 5 .

 5 Figure 5. Comparison of predicted (a) strain (ε y y ) (b) stress (σ y y = S y y + p) fields for the "grain" microstructure using CraFT-EVP and R 3 iCe.

Figure 6 .

 6 Figure 6. Kernel Density Estimation of the Von Mises equivalent strain (a) and stress (b) fields simulated with CraFT-EVP (EVP, black line) and R 3 iCe (CTI, cyan line). The dashed lines show results obtained for simulations where grains are defined by multiple mesh elements (grains microstructure) (figure 4 a-b). The full line shows a R 3 iCe simulation where a single orientation is attributed per element (figure 4 c)

  Uniaxial compression experiments from Montagnat et al. 2015 [29] (Σ ∼ 0.7 MPa, T=-5 • C ) • Uniaxial tension experiments from Jacka and Maccagnan 1984 [27] (Σ ∼ 0.4 MPa, T=-3 • C ) • Torsion experiments from Journaux et al. 2019 [31] (Σ max ∼ 0.5 MPa, T=-7 • C ),

Figure 8 .

 8 Figure 8. Simulation of uniaxial compression along Y with no recrystallization Γ R X = +∞, 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[START_REF] Montagnat | Measurements and numerical simulation of fabric evolution along the Talos Dome ice core, Antarctica[END_REF]; b) evolution with strain of the macroscopic strain rate; c) evolution of the intensity of the texture with the colatitude θ over a section defined on the Z plane (dashed blue lines in d)). Data at 22%, 40% and 60% macro strain are shown. The white dashed lines define the angle at ±45 • from the compression axis Y ; d) textures presented as pole figures for 22% (i ), 40% (i i ), and 60% (i i i ) 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).

Figure 9 .

 9 Figure 9. Simulations for uniaxial compression creep along the Y direction with parameters from table3and 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[START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF]; b) simulated macroscopic strain rate (black) and corrected simulated macroscopic strain rate (blue) to match a constant force boundary conditions instead of the constant stress one used in the simulation. EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude θ over a section defined on the Z plane (dashed blue lines in (d)). Data at 7%, 12% and 17.8% macro strain are shown. The white dashed lines define the angle at ±45 • from the compression axis Y ; d) simulated textures (pole figures) for 7% (i ), 12% (i i ), and 17.8% (i i i ) macroscopic strain; e) measured textures for the same macroscopic strains[START_REF] Montagnat | Analysis of dynamic recrystallization of ice from EBSD orientation mapping[END_REF]. All figures presenting texture data are color-coded using the same colorbar (scales in multiples of a uniform distribution).

Figure 10 .

 10 Figure 10. Texture obtained after 9% of octahedral tension strain (about 13% axial strain) performed at -3 • C with an octahedral stress of 0.4 M P a. The star ⋆ in the centre of the pole figure shows the direction of the tension axis. From Jacka et al. 1984 [27].

Figure 11 .

 11 Figure 11. Simulations for uniaxial tension creep along the Y direction with parameters from table3and 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 θ over a section defined on the plane normal to Z (dashed blue lines in d)). Data at 5%, 13% and 15% axial macroscopic strain are shown. The white dashed lines define the angle at ±45 • from the compression axis Y ; d) simulated textures (represented as pole figures) for 5% (i ), 13% (i i ), and 15% (i i i ) axial macro strain. All figures presenting texture data are color-coded using the same colorbar (scales in multiples of a uniform distribution).

Figure 12 .

 12 Figure 12. Distribution of the angle between the c-axis and the y-axis (tension axis) corresponding to the pole figure (ii) in figure 11d.

Figure 13 .

 13 Figure 13. Simulations for simple shear creep test (shear plane normal to the Y -axis and shear direction parallel to X , ∂u x ∂y < 0) with parameters from table3and comparison with experimental data from[START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF]. 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 Journaux et al. 2019[START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF]; b) simulated macroscopic strain rate. EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude θ over a section defined on the Z plane (dashed blue lines in d)). Data at γ = 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 γ = 0.2 (i ), γ = 0.42 (i i ) and γ = 1.96 (i i i ); e) measured textures for the same macroscopic shear strains[START_REF] Journaux | Recrystallization processes, microstructure and crystallographic preferred orientation evolution in polycrystalline ice during high-temperature simple shear[END_REF]. All figures presenting texture data are color-coded using the same colorbar (scales in multiples of a uniform distribution).

Figure 14 .

 14 Figure 14. Simulation of uniaxial compression creep test along Y at various stresses. Comparison between R 3 iCe predictions and experiments by Qi et al. [30]. (a) Textures from [30] obtained by EBSD for various strain rates leading to different quasi-steady state compressive flow stresses at macroscopic strains of ∼ 0.2% (exact strain is indicated at the top left of each pole figure). (b) Simulated textures obtained by creep under the same compressive stresses and up to the same macroscopic strains as in (a).

  Tr (M g D) = g D 33[START_REF] Alley | Fabrics in polar ice sheets -Development and prediction[END_REF] (1) δ 1 is given by applying a stress g S = ) = 0, (MS + SM) D = 0[START_REF] Castelnau | Modelling fabric development along the GRIP ice core, central Greenland[END_REF]

  ) = 0, (MS + SM) D = S[START_REF] Llorens | Full-field predictions of ice dynamic recrystallisation under simple shear conditions[END_REF]

-•-

  using the constitutive equation, we obtain for g D 33 :g D 33 = δ 1 Then a uniaxial stress σ is applied normal to the c-axis, for instance in the direction 1. The stress tensor is g σ = using the constitutive equation, we obtain for g D 11 :g D 11 = δ 1

Figure 15 .

 15 Figure 15. Kernel Density Estimation of all components of the strain tensor predicted by an elasto-viscoplastic (EVP) simulation using CraFT and the Continous Transverse Isotropic (CTI) using R 3 iCe after 0.01 macroscopic strain. The dashed lines show results for simulations where grains are well discretized. The full line shows result for the simulation with one orientation per element.

Figure 16 .

 16 Figure 16. Kernel Density Estimation of all components of the stress tensor predicted by an elasto-viscoplastic (EVP) simulation using CraFT and the Continous Transverse Isotropic (CTI) using R 3 iCe after 0.01 macroscopic strain. The dashed lines show results for simulations where grains are well discretized. The full line shows result for the simulation with one orientation per element.

Table 1 .

 1 Mathematical notations

	Tr ()	Trace of a tensor	Tr (A)	A i i
	d i m() dimension of a tensor	d i m(A)
	∇	∇ operator	∂ ∂x , ∂ ∂y , ∂ ∂z	∂ ∂x i
	.	Scalar product	a.b	a i b i
	:	Double contracted product A : B	A i j B i j
	∥.∥	Norm	∥a∥	a 2

i A n matrix product n times A × A... × A n

Acknowledgements

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

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

Funding

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

per element are performed, starting with initial orientations randomly picked in a uniform texture. The parameters for all simulations are given in table 3. For all configurations, the simulated texture evolution is represented by (1) the evolution of the eigenvalues a i of the second order orientation tensor a (2) (a

with N the number of orientations), (2) pole figures representing the c-axes distribution at various strains, and (3) the evolution with strain of the texture in a section extracted from the pole figure. These outputs are compared with the pole figures obtained during the experiments at three selected strain values.

The second order orientation tensor a (2) is a statistical representation of the texture anisotropy [START_REF] Woodcock | Specification of fabric shapes using an eigenvalue method[END_REF]. The c-axes distribution lies within an ellipsoid whose axes are the eigenvalues of a (2) with 1 ≥ a 1 ≥ a 2 ≥ a 3 > 0.

Conflicts of interest

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

Authors' contributions

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

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

Véronique Dansereau participated to the elaboration of the scientific question and codeveloped the model with TC. She participated to the writing of the manuscript.

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

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

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

All authors read and approved the final manuscript.

Appendix A. Demonstration of δ i identification

A.1. Constitutive equation CTI and single crystal behavior

The constitutive equation is:

The single crystal anisotropy can be described by three parameters derived for experimental data :

ψ 1 is the fluidity for parallel shearing to the basal plane :

β is the viscosity ratio between shear parallel to the basal plane η and the viscosity for shear within the basal plane :

γ is the viscosity ratio between compression (or traction) along the c axis η c and compression (or traction) in one direction g e r within the basal plane:

The same compressive (or extensive) stress parallel and normal to the c-axis S = g S 33 = g S r r results therefore in different strain rates:

g D r r = γ g D 33 (41)

A.2. δ i identification

To identify the parameters (δ 1 , δ 2 , δ 3 ) in the constitutive equation [START_REF] Treverrow | The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement[END_REF], we developed it and compared to the single crystal behavior. Using c = (0, 0, 1)