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

physically-based formulation for texture evolution during


Introduction
Dynamic recrystallization (DRX) is a mechanism by which local stress incompatibilities and local strain heterogeneities are relaxed [1][2][3].By local, we mean the scale of the grains that constitute crystalline materials such as ice, metals, and most rocks.Dynamic recrystallization occurs in plastically deformed materials at high temperatures.It is characterized by the nucleation of new grains and grain boundary migration and can lead to the complete re-organisation of the microstructure (grains shape and size) and drastic modification of the texture (crystallographic orientations of grains) [2,[4][5][6][7].These changes in the microstructure and texture may lead to modification of the rheological behavior, i.e., the response of the material to imposed stresses or strains.
Ice is a hexagonal crystalline material (Figure 1) that deforms on Earth at temperatures very close to its melting point.In the ductile regime, its deformation is characterized by a strong anisotropy, with dislocations gliding almost solely on the basal plane [8,9].Therefore, the orientation of an ice crystal relative to the direction of solicitation strongly controls its viscosity.Consequently, the distribution of crystallographic orientations in polycrystalline ice impacts the mechanical response in the form of texture-induced viscoplastic anisotropy [10].
Microstructural evidence for recrystallization is systematically observed along deep ice cores, where it is supposed to contribute variably to the measured texture, depending on the deformation conditions (temperature, impurities).To interpret texture evolution along ice cores, a model has been proposed by De la Chapelle et al. 1998 [11].This model states that DRX may either slow down the strain-induced texture evolution (the so-called "continuous" dynamic recrystallization regime) or drastically modify the texture, a final texture controlled by the stress state (the "discontinuous" dynamic recrystallization regime).Whether continuous or discontinuous DRX predominates depends on the temperature and strain conditions along the ice core [12,13].In deep sections of ice sheets, along ice streams, and in glacier areas where temperature and strain are high, DRX leads to highly anisotropic textures, particularly when simple shear is dominant (see e.g.[14][15][16][17]).These textures result in strong viscoplastic anisotropy, which can either enhance or slow down strain rates.For instance, the vertical cluster of c-axes usually observed along ice cores will accelerate horizontal shearing but slow down vertical compaction [18][19][20][21][22][23].To model the deformation of an ice sheet, it is therefore essential to consider the effect of DRX on texture evolution.
Dynamic recrystallization in ice has been extensively studied in the laboratory.These studies allowed us to characterize the basic mechanisms at play, nucleation and grain-boundary migration, and the impact of strain heterogeneities on their kinetics (see e.g.[3,[24][25][26][27][28][29][30][31]).Kamb (1972) [4] was a pioneer in performing several creep experiments in simple shear and compression to characterize the effect of DRX on both texture evolution and macroscopic mechanical behavior.He proposed that the strain rate acceleration, which is characteristic of the tertiary creep stage, could result from the development of DRX textures favorable to the activation of basal slip.This study is likely at the origin of the widely established belief that DRX favors "well-oriented" grains relative to the imposed stress.Since then, work by Grennerat et al. 2012 [32] has revealed that the strain field during deformation in polycrystalline ice is highly heterogeneous, with no simple relationship between the intensity of the local strain and the Schmid factor, a parameter that quantifies the orientation relationship between a crystal and the imposed macroscopic stress.More recently, direct measures of the evolution of the strain field with DRX by Chauve et al. 2017 [33] have illustrated the heterogeneity of the mechanical behavior at the grain scale and its impact on DRX.They showed that DRX tends to homogenize the strain field in relation to local microstructure changes.
Although this recent data show that DRX textures result from a complex interplay between the imposed stress or strain states, the local microstructure, and the resulting redistribution of stresses, it also corroborates the notion that the development of DRX textures is systematically associated with marked strain softening that signs the onset of tertiary creep [22,25,27].Further texture evolution during tertiary creep results in a quasi-constant strain rate for bulk strains higher than 10%; this "stationary" strain rate is independent of the initial texture and controlled by the imposed stress [34][35][36][37].This contrasts with the minimum strain rate observed at the onset of secondary creep, which depends on the initial texture and can be described by the classical Glen's law (a Norton-Hoff type law), which nonlinearly relates the strain rate to the applied stress (through an exponent of 3 and an Arhenius type of dependence to the temperature) [9,38].This law does not hold in tertiary creep, but it has been adapted by modifying the stress exponent to a value slightly higher than 3 and adding a strain-rate enhancement factor (the ratio between the strain rate measured at tertiary creep and the minimum strain rate at secondary creep).This enhancement factor, between 4 in compression and 8 in simple shear [22,37] accounts for the impact of the DRX texture on strain softening during tertiary creep.
Textures that develop in the ice during DRX are well characterized.Under uniaxial compression, the textures measured at tertiary creep are characterized by c-axes oriented within a girdle at an angle between 30 • and 50 • to the compression axis depending on the experimental temperature T and stress.These textures are stable, i.e., they evolve very slowly for strains higher than 10% [4,24,27,29,30,37,39].During simple shear, tertiary creep textures seem to evolve with increasing strain from a two-maxima pattern, which is symmetric relative to the principal finite stretching direction, toward a strongly clustered texture, characterized by c-axes aligned in a direction perpendicular to the shear plane [25,26,31,40,41].
Texture development in ice deforming by dislocation creep has been successfully simulated using approaches of variable complexity [18,19,42,43].However, these approaches fail to reproduce the tertiary creep experimental textures, and the rare models simulating the effect of dynamic recrystallization are rather unsatisfactory in terms of either texture patterns or kinetics of texture evolution [44][45][46][47][48].The discrepancies between models and observations result from the use of a simplified, mean-field approach that does not enable DRX mechanisms to be properly represented [44,45], or from limitations in the deformation modeling frame despite the high complexity of microstructure evolution [46].In [47,48], the simulated textures are very close to the observed ones, but the kinetics of evolution are not.This discrepancy may result from the use of a Sachs approximation for modeling the mechanical response of the polycrystal, which results in a very soft mechanical response and incoherence in the grain rotation calculation [49].
In this work, we propose a new formulation to predict the impact of DRX on the texture evolution of polycrystalline materials such as ice.This formulation couples a continuous transverse isotropic (CTI) law, which models the anisotropic response of the ice crystal, to a crystal orientation evolution equation that accounts for DRX.The contribution of DRX to texture evolution is modeled by a physically based orientation attractor that maximizes the resolved shear stress on the easiest slip system in the crystal (basal slip for ice).This set of equations, named CTI-RX, is solved in a full-field model using the finite-element method (FEM).
In Section 2, we present the CTI-RX model: its equations for the CTI law (Section 2.1) and the proposed formulation for the evolution of the c-axis (Section 2.2).The numerical scheme designed to solve these equations is described in Section 3. Section 4 presents the calibration and validation for cases without dynamic recrystallization.Section 5 evaluates the capability of the model to reproduce the evolution of the texture and the mechanical response by comparing the predictions with experimental data for tertiary creep of ice polycrystals subjected to uniaxial compression, uniaxial tension, and simple shear.

Continuous transverse isotropic behavior
Following [50][51][52][53][54], we account for the viscoplastic anisotropy of the ice single crystal by modeling its mechanical behavior using the continuous transverse isotropic (CTI) law.In this framework, the viscoplastic response of the ice single crystal depends only on the orientation of the c-axis, being isotopic in the basal plane (see Figure 1).This formulation is based on theoretical work by Boehler (1987) [55], which has shown, through geometrical arguments, that the relation between the deviatoric stress S and the deviatoric strain rate D in a transverse isotropic symmetry may be described by a potential φ S : where φ S is expressed on the B C T I base, which is built from a combination of the following four invariants: 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 nonlinear (φ (n) S ) cases.
e 3

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

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: Considering the following derivation rules, the constitutive relation is given by The three independent parameters δ i in Equation ( 5) can be fitted from the experimental data as detailed below.
In the present implementation of the model, the anisotropy of the ice single crystal (or grain) behavior is described using the following parameters, where 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: • β is the ratio between the fluidity for shear parallel to the basal plane and that for shear normal to the basal plane: • γ is the ratio of 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: 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 (9) Using the above constraints on single crystal behavior, identification gives (see Appendix A for details): The constitutive equation ( 5) can also be expressed in an inverse manner as follows: and in terms of the potential with : ) and η 1 = 1 ψ 1 (see [19] and appendix A for details).

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 : This hypothesis limits the number of independent parameters to 3 instead of 7 in the general case.A similar hypothesis was made in previous studies [19,48,56].Using the chain rule for deriving φ n D gives: k is further obtained by identification of the power exponent on Equation ( 14) in order to have Therefore, a non-linear CTI law can be written as: Leading to the constitutive equation: with effective viscosity 1−n 2n (18) and α i parameters: It is worth noting that β links g D 12 to g S 12 through Therefore, β depends on non-linearity, i.e., on the value of n.A similar enhancement is obtained for a linear behavior n = 1 and a non-linear behavior with n = 3 by setting β n=3 = β n=1 .

Evolution of crystal rotation in response to deformation and recrystallization
Deformation resulting from dislocation glide in the different slip plane results in a rotation of the crystal whose rate is described as in [57,58] by where are the spin-and strain-rate tensors.According to Equation (21), the rotation of the c-axis is due to the bulk spin, Wc, and the viscoplastic spin expressed by λ[Dc − (c T Dc)c].λ is a parameter that depends on the plasticity model chosen for the ice crystal.In this frame, we assume that deformation is only due to the contribution of dislocations gliding on the basal plane, λ is set to 1.A value of λ < 1 must match a plasticity model that include prismatic and pyramidal slip [57,58].This assumption is well funded for ice and has often been made in previous studies [45,[59][60][61][62].
To account for the impact of dynamic recrystallization on the evolution of the c-axis orientation, we add a new term to this equation based on the assumption that dynamic recrystallization induces a local re-orientation of the c-axis, either through nucleation of a grain with a new orientation or by migration of a grain boundary, and that this re-orientation leads to local softening.This continuous description of the impact of DRX on the evolution of the c-axis orientation is supported by experimental observations of DRX nucleation in ice.These observations show that nucleation by bulging and polygonization results in progressive misorientation of the recrystallized volumes relative to the parent grains, with misorientations occurring most often < 20 • [33].
With the addition of the dynamic recrystallization term, the c-axis rotation is described as follows: where Γ R X is a parameter that controls 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.Therefore, the evolution of the c-axis orientation is controlled by the balance between (i) the deformation-induced rotation and (ii) the rotation toward a recrystallization-driven attractor, c 0 .The meaning of Γ R X can be illustrated by considering the artificial case in which only dynamic recrystallization contributes to the c-axis rotation.In this case, Equation (22) becomes The norm of ∂c ∂t can be interpreted as the rotation rate due only to dynamic recrystallization.This rate increases as ∥c 0 − c∥.Considering a constant c 0 , the solution for the c-axis rotation is: 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 diatoric 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 toward, or/and nucleates in, an orientation that facilitates local viscoplastic deformation (dislocation glide).This assumption is supported by experimental observations that show that nucleation and grain boundary migration occur in areas of high strain incompatibility and produce orientations that tend to reduce this incompatibility [3,33,63].Considering that the ice single crystal deforms mainly by dislocation glide on the basal plane, the recrystallization-induced rotation should maximize the resolved shear stress (RSS) on the basal plane.RSS for a dislocation with a Burgers vector a i (Figure 1) gliding on the basal plane c in response to a diatoric stress S is given by RSS(a i , c, S) = S : µ, (25) 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 as the eigenvalues and v i (i = 1 − 3) as the corresponding eigenvectors of the diatoric stress tensor S, the maximum resolved shear stress is ∥s 1 − s 3 ∥ and the associated orientations are [64]: Degenerated cases exist if at least two eigenvalues are equal: • s 1 = s 2 , the possible attractors c 0 are the orientations at 45 • from v 3 .
• s 2 = s 3 , the possible attractors c 0 are the orientations at 45 • form v 1 .
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 recrystallization on the orientation of the c-axis.In the following section,we describe three end-member cases.
For uniaxial compression along the y-axis (Figure 2 (a)), the stable solution for the deformation term (Wc − λ[Dc − (c T Dc)c]) corresponds to c aligned with y (red dot, Figure 2 (a)) [19].For recrystallization ( 1 Γ R X (c 0 − c)), the stable solutions are the orientations at 45 • to the y-axis (blue line, Figure 2 (a)).When both deformation and recrystallization are active, the stable orientations of the c-axis are located on a girdle between these two end-members (black line, Figure 2 (a)).The position of this girdle depends on the value of Γ R X (black arrows, Figure 2 (a)).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 2 (b)) [19].When both deformation and recrystallization are active, the stable orientations (black circle, Figure 2 (b)) form a girdle between 45 • and 90 • to the y-axis depending on the value of Γ R X .
For simple shear parallel to the xOz plane (with ∂u x ∂y < 0, Figure 2 (c)), there are four stable solutions of the recrystallization term that are ±y and ±x (Figure 2 (c), blue dots), and two for the deformation term, ±y (Figure 2 (c), red dots).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 2 (c)).This latter stable orientation, whose exact position depends on the value of Γ R X , only exists for low Γ R X and vanishes for large  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).

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: (2) The non-linear CTI law: (3) The equation for the evolution of the c-axis: With 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 by creep (stress boundary conditions) under either uniaxial compression, uniaxial tension, or simple shear (Figure 3).A full-field approach using FEM was chosen to describe the field of orientations 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 [29][30][31].

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

Simple shear creep tests
In the simple shear tests, forcing is applied as a constant tangential stress on the top boundary (F x (top) = F ) of the sample.Other boundary conditions are: • u(bottom) = 0, • u y = u z = 0 for the lateral surfaces (left, right, front, back) to avoid solid rotation,

Characteristic Numbers and time
When adimensionalized with respect to the material viscosity η n in Pa s 1 n , the vertical extent of the cuboid L in meters (Figure 3) and the applied macroscopic stress Σ (in Pa), the system of equations becomes: 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 is the ratio between the deformation characteristic time τ = ( ) 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
Because the evolution of the c-axis orientation depends on the strain rate ( D( u)) 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 this problem, we linearize the equation describing the evolution of the c-axis using a second-order backward differentiation formula (BDF2).The time discretization of the CTI-RX model results in the following system of equations: Finite elements and variational methods were used to solve the time-discretized problem on a Lagrangian grid within the C++ environment Rheolef [65].The domain is spatially discretize using hexahedral finite elements (unless otherwise mentioned), and an initial c-axis orientation is prescribed for each element.The FE approximation for the tensorial mechanical parameter M and the scalar mechanical parameter c is computed with polynomials of the order 0 (P 0 ).The FE approximation for the tensorial strain rate D, spin rate W and stress S fields is computed with polynomials of order 1 (P 1 ).The FE approximation for the displacement field is computed with polynomials of the order 2 (P 2 ).
The Rheolef problem_mixed solver is used to solve the displacement field u and the pressure field p under the incompressibility hypothesis for the linearized version of Equation (35), 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 (36).In practice, it is recommended for numerical efficiency to flatten the two nested fix loops.
?? 1 presents the numerical algorithm of two nested fix points implemented in R 3 iCe.This algorithm is a semi-implicit solver for the coupled Equations ( 31) and (32).The simulated orientation field c is thus consistent with the displacement field u and therefore with the diatoric strain rate D and the deviatoric stress S.This implementation of CTI-RX model is distributed as the R 3 iCe code.
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.

CTI parameter values for the ice single crystal
As stated in Section 2, the values for the CTI parameters, γ, η n , n and β (Equations ( 17) and ( 19)) are defined on the basis of experimental data.γ corresponds to the ratio of viscosities measured in compression and tension for an ice single crystal that is well oriented for basal slip.Laboratory experiments have shown that this ratio is very close to one.Therefore, in the present numerical simulations, we set γ = 1, as in previous studies [50,51].
The following Duval et al. (1983) [9] and previous modeling work based on the CTI law [20,48,51], we set n = 3.While experiments in ice single crystals that are well oriented for basal slip are best fitted by n = 2, experiments on crystals that are poorly oriented for basal slip result in n = 3 [9].Moreover, a value of 3 must be imposed in the CTI formulation to recover n = 3 for the polycrystal response.
β is linked to the ratio between the basal and non-basal shearing viscosities.The strong viscoplastic anisotropy of ice is mainly carried by this parameter.In creep experiments at 1 MPa [9], the ratio between basal and non-basal shear viscosities is on the order of 10 4 .This ratio is related to β through Equation ( 20), β n+1 2 = 10 4 .For n = 3, β = 10 −2 is a value that is consistent with those considered in previous studies using the non-linear CTI formulation [20,48,51].
The viscosity η n can be adjusted by comparison with either the experimental response of a single crystal or that of isotropic polycrystals.Solving the CTI law for the uniaxial compression creep of a single crystal oriented at 45 • from the compression axis and of a polycrystal with random orientations picked from a uniform texture results in a macroscopic strain rate of D si m sc = 1.7 × 10 −2 ( ) MPa s 1/3 for the polycrystal.The value of 7.5 MPa s 1/3 is taken for η 3 in order to reproduce the experiment duration as observed during uniaxial compression tests that are taken here as a reference (see Section 5.2).The parameters of the CTI law for the ice single crystal are summarized in Table 2.

Validation of the CTI-based implementation for predicting 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 for improving the CTI formulation with a representation of DRX.The mechanical response of polycrystalline ice is known to be characterized by strong strain and stress heterogeneities, which are very likely precursors of dynamic recrystallization [32].We first compare the R 3 iCe predictions to those of the CraFT-EVP model [32], which has been shown to adequately reproduce the heterogeneous strain fields observed in experiments [32,66].We then investigate the reliability of a simplified mesh configuration for simulating the evolution of the texture and its impact on the macroscopic mechanical response in a numerically efficient manner.

Stress and strain field predictionscompared with CraFT-EVP
CraFT-EVP has been used to provide full-field predictions of the mechanical response of ice during primary and secondary creep [32].It is based on an elasto-viscoplastic formulation proposed by Suquet et al. 2012 [67] and solved using CraFT software.The predictions of CraFT-EVP for the mechanical response and strain fields were validated by direct comparisons with mechanical tests on ice polycrystals that included strain field measurements by Digital Image Correlation [32,67].The macroscopic strain fields predicted by CraFT-EVP agreed well with the experimental observations at 1% strain (secondary creep, before any texture evolution owing to DRX).With no texture evolution, CraFT-EVP and R 3 iCe both simulate a stationary secondary creep.
Reproducing accurate stress fields in R 3 iCe is key because the recrystallization attractor c 0 depends on the deviatoric stress S at the element scale (Section 2.3).Note that the Elasto-ViscoPlastic (EVP) law implemented in CraFT-EVP [67] is based on a full description of crystal plasticity, in which all possible slip systems contribute to the deformation as a function of their critical resolved shear stress and crystal orientation, whereas the CTI law in R 3 iCe is a simplified parameterization of viscoplastic anisotropy, which has been calibrated directly from experimental data (Section 3.4).The objective of the comparison performed below is to verify that the stress field predicted by R 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 the two predictions can be considered independent.
Figure 4 (a-b) shows the polycrystal microstructure generated by a grain growth algorithm using Neper [68] used for both CraFT-EVP and R 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 was 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 perfectly compare between the CraFT-EVP and R 3 iCe predictions although the overall patterns do not show a strong mismatch.The locations of areas in compression (blue) and extension (red) are similarly prescribed (Figure 5 (a)).Significant differences between the responses of the models CraFT-EVP and R 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 surface results in a constant force applied on each element that can lead to heterogeneous displacement on this surface.As a result, some "weaker" grain deform faster, as represented by the dark blue areas on Figure 5 (a).In CraFT-EVP, periodic boundary conditions are imposed because equations are solved with a fast Fourier transform-based method.The simulated stress fields (Figure 5 (b)) compare satisfactorily in terms of amplitude and global pattern, although the precise stress localization shows more scattering.Some of this mismatch may also be explained by the differences in boundary conditions between the two models.
Similarly, the distributions of predicted Von Mises equivalent stress and strain by the two models show some discrepancies (Figure 6).These discrepancies are stronger between the two configurations in which grains are described by multiple elements (Figure 6) and where the impact of the boundary conditions mentioned above is stronger.Distributions for all components of both tensors ε 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 .The 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 7 (a)), and R 3 iCe with a configuration with several elements per grain (Figure 7 (b)) and one orientation per element (Figure 7 (c)).All three configurations predict similar v 1 and v 3 orientation distributions and therefore similar c 0 orientation distributions are expected.
In summary, if CraFT-EVP simulations are taken as a reference, the strain and stress fields at secondary creep, prior to DRX activation, are qualitatively and statistically well predicted by the CTI constitutive law implemented in R 3 iCe.In particular, the principal stress tensor directions used in R 3 iCe to simulate the dynamic recrystallization impact on orientation evolution are similarly predicted by the two models.

Stress and strain field prediction in the simplified mesh configuration
Resolving CTI-RX with sub-grain scale meshes is numerically very costly.In the following, we test a microstructure representation composed of one orientation (or grain) per mesh element (Figure 4 (c)) in order to keep a statistically reasonable number of orientations on a mesh size (∼ 2700 elements).
The distributions of the Von Mises equivalent stress and strain in these one-orientation-permesh-element simulations are compared with those of the simulations with meshes at the subgrain 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, which is here considered as a reference.The fact that the CTI-RX fields predictions obtained with a sub-grain scale mesh depart slightly from the Craft-EVP predictions is attributed, as mentioned before, to the discrepancy in the applied boundary conditions (periodic versus Neumann applied on the top surface).Please note that this discrepancy is reduced by a larger number of grains (orientations) in contact with the top surface in the simplified microstructure.
Figure 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.There is a superb 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 to have a large enough number of orientations while maintaining reasonable computation times.We assume that this choice has a limited impact on the simulated stress fields and therefore on the evaluation of the attractor c 0 at the sample scale.

Texture evolution: comparison with experimental data
This section tests 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 because of the fast strain rates and high temperature required.Deformation-only textures may only be sampled in deep ice cores from the central parts of Greenland or Antarctica [11].These textures were well reproduced by uniaxial compression simulations performed using the mean-field ViscoPlastic Self Consistent model [13,43].R 3 iCe prediction for texture evolution without recrystallization are therefore compared with those of VPSC simulations extracted from [13].In contrast, there is a large experimental dataset in which the texture evolved in response to both deformation and DRX.We use this dataset to constrain the predictions of R 3 iCe for the effect of DRX on texture evolution and mechanical behavior under uniaxial compression, uniaxial tension, and simple shear.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 [30].These conditions enabled us 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 per element are performed, starting with initial orientations randomly selected in a uniform texture.The parameters for all simulations are given in Table 3.
For all configurations, the simulated texture evolution is represented by (1) the evolution of the eigenvalues a i of the second-order orientation tensor a (2) (a (2)    The second-order orientation tensor a (2) is a statistical representation of texture anisotropy [69].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.
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 the absence of recrystallization
To test the capability of R 3 iCe to predict texture evolution, we simulated uniaxial compression under a 1 MPa deviatoric stress using one orientation per mesh element with initial orientations randomly selected from a uniform texture.
Figure 8 (a-c) presents the evolution with strain (i.e.time) of the eigenvalues (a i , i = 1, 3) of the second-order orientation tensor a (2) , the strain rate, and the texture in the Z plane (perpendicular to the compression direction).Figure 8 (d) shows textures predicted at 22%, 40% and 60% total strain.The strain-rate curve discontinuities (Figure 8 (b)) are technical artifacts imposed by the operating rules of the GRICAD clusters.The shared access restricts the simulation duration to 48 hours, which means that the code must be relaunched every 48 hours.To reduce the calculation time, the convergence criterion is set tougher for the first step and relaxed for the following steps.The steps observed at each restart reveal the convergence error during the simulation, which appears to have a limited impact.
A strong cluster (single maximum of c-axes) texture develops as strain increases.This texture evolution induces geometrical hardening of the polycrystal, which is expressed by a decrease in the strain rate by a factor ∼ 1.5 between 20% and 60% total strain, when the texture is highly clustered (a 1 ≈ 0.71).
Similar cluster textures are observed along deep cold ice cores, for instance, along the Talos Dome ice core [13] where a clustered texture with a 1 ≈ 0.7 was measured at about 650 m depth, where an accumulated strain of 60% is estimated.The simulated texture evolution with strain is also in good agreement (although slightly weaker) with the one predicted by Montagnat et al. [13] using the VPSC model when departing from an isotropic texture (see red dots in Figure 8 (a)).

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.This is compared with experimental measurements from Montagnat et al. 2015 [29].These latter were obtained from initially isotropic granular ice compressed at −5 • C under a constant load with initial stresses of 0.7 and 0.8 MPa, and up to various strains (from 2 to 17.8% bulk shortening).
Comparison of the c-axis pole figures shows good qualitative agreement between the simulated and experimental textures (Figure 9 (d-e)).Both are girdle textures that become more concentrated and show an evolution of the cone angle to the compression axis from 45 • to 35 • up to ∼ 20% macroscopic strain and remain nearly constant thereafter (Figure 9 (d)).There is also good agreement between the predicted evolution of the texture eigenvalues and those measured ).Red dots are predictions from the VPSC model taken from [13]; b) evolution with strain of the macroscopic strain rate; c) evolution of the intensity of the texture with the colatitude θ over a section defined on the Z plane (dashed blue lines in d)).Data at 22%, 40% and 60% macrostrain are shown.The white dashed lines define the angle at ±45 • from the compression axis Y ; d) textures presented as pole figures for 22% (i ), 40% (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).at 7%, 12% and 17.8% (Figure 9 (a)).This good agreement between the 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 shorter time (∼ 130 hours), but the tertiary creep (maximum strain-rate) regime is reached between 10 to 20% strain, similar to the experimental observations (see also [27,37]).
The macroscopic response predicted by R 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 a cone angle around ∼ 45 • to the compression axis.This softening is followed by a slight hardening as the girdle texture slowly evolves toward a cone angle at ∼ 35 • to the compression axis.
The bulk strain associated with the maximum strain rate and the simulated enhancement factor cannot be directly compared to that measured in Montagnat et al. 2015 [29] because these experiments were performed under constant load and the area of the sample section evolved during the tests, whereas the model simulations were performed in a constant stress configuration.The simulated strain rate may be corrected by assuming a constant sample volume during the experiment, which leads to εcor r si m = εsim × (1 + ε) 3 .This corrected strain rate is plotted in blue in Figure 9 (b).This brings the maximum strain rate to ∼ 10% and the enhancement factor to 3.0.These values are closer to the macroscopic response measured by Montagnat et al. [29], which is characterized by an enhancement factor between 3.3 and 3.8.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 [29]; b) simulated macroscopic strain rate (black) and corrected simulated macroscopic strain rate (blue) to match the constant force boundary conditions instead of the constant stress one used in the simulation.EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude θ over a section defined on the Z plane (dashed blue lines in (d)).Data at 7%, 12% and 17.8% macrostrain are shown.The white dashed lines define the angle at ±45 • from the compression axis Y ; d) simulated textures (pole figures) for 7% (i ), 12% (i i ), and 17.8% (i i i ) macroscopic strains; e) measured textures for the same macroscopic strains [29].All figures presenting texture data are color-coded using the same color bar (scales in multiples of a uniform distribution).

Uniaxial tension
Experimental creep tests performed under 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 to the definition given in Jacka et al. 1984 [70], to 13.1% axial strain.At this strain level, they did not attain the quasi-constant strain rate typically of tertiary creep.The texture obtained, measured using the classical manual Rigsby stage technique with one orientation per grain, is presented in Figure 10.The texture is characterized by a small circle girdle with a mean half angle of 50.4 • and a standard deviation of 15.3 • .Although tertiary creep was not reached in this experiment, an enhancement factor of ∼ 3, similar to that observed in compression under the same conditions, was measured.
The texture evolution and macroscopic response predicted by a simulation under uniaxial tension creep are shown in Figure 11.The texture evolves with increasing strain toward a girdle texture with a cone angle slightly around 45 • to the extension direction.This texture appears to stabilize at ∼ 5% axial strain.The distribution of the c-axis angle from the tension axis is given in Figure 12 for the simulated texture at 13% strain.It is similar to the one measured by Jacka et al. 1984 [27] (Figure 10).An enhancement factor of approximately 10 is simulated at 10% axial strain.

Simple shear
The simulated texture evolution in simple shear is compared with observations made by Journaux et al. [31] for torsion experiments on isotropic granular ice performed at −7 • C, under a constant torque corresponding to a maximum shear stress between 0.4 and 0.6 MPa.
The experimental textures are characterized by two submaxima, one almost perpendicular to the shear plane, called M 1, and the second one, M 2, initially at a low angle to the shear plane.Similar evolution was also observed in simple shear in [26,41,71].As the shear strain increases, 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 toward 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, is also well reproduced (Figure 13 (a)).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 strains are shown.The white dashed lines define the angle at ±45 • from the compression axis Y ; d) simulated textures (represented as pole figures) for 5% (i ), 13% (i i ), and 15% (i i i ) axial macrostrain.All figures presenting texture data are color-coded using the same color bar (scales in multiples of a uniform distribution).
In both experiments and simulations, the strain rate increased rapidly up to γ = 0.3 while the double-maxima texture develops (Figure 13 (d)).It continues to increase, but at a progressively slowing rate as the texture evolves toward a single maximum normal to the shear plane, up to γ ∼ 3 where is reaches a quasi-steady state.As explained in Section 5.2, the discontinuities in the strain-rate curve (Figure 13 (b)) are technical artifacts imposed by the operating rules of GRICAD clusters.
In phase with Treverrow et al. 2012 [37], we define an enhancement factor for simple shear as the ratio between the strain rate during secondary creep (here it corresponds to the beginning of the run) and that obtained when the strongest simulated texture is reached.Using this definition, we obtain an enhancement factor (E F ) of 6.7 that matches Treverrow et al. 2012 [37] data for a similar applied shear stress.

Testing the relative effects of deformation and recrystallization on texture evolution under high stresses
The work of Qi et al. [30] provides some of the rare existing highly resolved textures obtained for temperatures and types of samples similar to those of the experiments simulated above, but under imposed displacement rates that induce higher stress conditions (see Table 3).The stress values reported in Table 3 correspond to the flow stress, i.e., the nearly constant stress reached after peak stress, at about 20% strain.To reach deviatoric stresses higher than 1 MPa without failure, a confining pressure of 10 MPa was applied.
To test the sensitivity of the dynamic recrystallization formulation to the level of stress, we performed simulations under creep conditions using the stress values from Qi et al. 2017 [30] experiments.As CTI-RX equations are relating the deviatoric stress to the deviatoric strain and are solved under the Stokes hypothesis (incompressibility), it is not necessary to numerically apply the confining pressure.Although boundary conditions differ between the simulations (macroscopic constant stress) and experiments (macroscopic constant displacement rate), we assumed the quasi-steady state behavior sampled by the flow stress measurements in the experiments to be equivalent to the quasi-steady state part of the tertiary creep in the creep simulations [72].The stable textures expected at both stages can therefore be compared.
A comparison between the simulated and observed textures at different compressive flow stresses is shown in Figure 14.The experimental observations are well reproduced.Both measured and simulated textures show a transition from a girdle texture to a weak single maximum with increasing flow stress for a macroscopic strain of ∼ 20%.To correctly fit the experimental data, a value of Γ 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 [29] (Table 3).This implies a higher contribution of recrystallization compared with 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 correctly reproduces the textures observed along deep ice cores where dynamic recrystallization is assumed to have a low impact.3 and comparison with experimental data from [31].a) Evolution with strain of the eigenvalues of the second-order orientation tensor (a 1 > a 2 > a 3 ).The blue dots correspond to the experimental values from Journaux et al. 2019 [31]; b) simulated macroscopic strain rate.EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude θ 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 [31].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 is based on the simulation of the evolution of the texture (orientations of c-axes) based on a new formulation that (i) introduces an attractor controlled by the local stress field to reproduce the impact of DRX mechanisms, (ii) is coupled with the continuous transverse isotropic (CTI) law to account for strain and stress field heterogeneities at the crystal scale arising from the crystal viscoplastic anisotropy, and (iii) solves these coupled equations using a semi-implicit numerical scheme in a full-field implementation using FEM, ensuring consistency between texture evolution and mechanical behavior.The first two ingredients are derived from experimentally based knowledge of dynamic recrystallization processes.
Indeed, recent experimental observations have (i) clearly documented the role of the crystal viscoplastic anisotropy in generating local stresses and strains markedly different from the macroscopic ones [3,32,66] and (ii) revealed that the formation of new orientations (nucleation) during dynamic recrystallization occurs through bulging (heterogeneous grain boundary migra- tion) and polygonization (formation of new grains boundaries by organization of dislocations within a grain) [3,33,63,66,73].These two processes are controlled by the local strain and stress fields and result in recrystallized orientations that are strongly related to, but slightly different from, the parent grain orientations.The last important experimental observation is that the development of recrystallization textures leads to a weakening of the polycrystal mechanical response.All these observations justify the choice of a continuous formulation for the c-axis orientation evolution (Equation ( 22)), with a targeted recrystallized orientation that is calculated to maximize the local basal resolved shear stress (Equation ( 26)).
Comparison with tertiary creep experiments performed in compression [29,37,39], in simple shear [31,37,41] and in tension [27] supports the fact that the R 3 iCe model can reproduce texture evolution in conditions similar to those of laboratory experiments (T = −5 to −10 • C, S ∼ 1 MPa).For this, a single parameter is adjusted: the recrystallization parameter Γ R X , which weights the relative contribution of deformation and recrystallization to the rotation of the c-axes.As expected, under similar stress conditions, a lower temperature requires a higher Γ 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 also correctly simulated, compared with experimental results performed under similar conditions [22,27,29,36,37].Because the present model does not consider any softening mechanism associated with dislocation interactions or grain boundary migration, this comparison seems to confirm that texture-induced viscoplastic anisotropy explains most of the mechanical softening measured during tertiary creep.Therefore, softening associated with dislocation annihilation and grain boundary migration during dynamic recrystallization appears to play a secondary role.
Simulations at higher uniaxial compressive stresses (from 1 to 4 MPa) have been performed to compare with the experiments of Qi et al. 2017 [30].Contrary to the previously mentioned experiments, these experiments were performed under a range of imposed strain-rate conditions, which resulted in different flow stresses.A confining pressure of 10 MPa was applied to prevent failure at the relatively high strain rates considered.The observation of more clustered textures at higher flow stress implies that higher strain rates result in a larger impact of deformation relative to DRX on texture evolution.R 3 iCe run under imposed stress without confining pressure can correctly represent texture clustering under stress.Indeed, experimental textures measured at 20% strain are satisfactorily reproduced by R 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 counterbalance the impact of deformation on the c-axes rotation.We interpret this apparent contradiction as due to the activation of brittle-ductile processes in the experiments.This assumption is based on observations by Kalifa et al. [74], which performed triaxial tests under the same range of strain rate (∼ 10 −5 − 10 4 s −1 ) than Qi et al. 2017 [30], but under variable confining pressure from 0 to 10 MPa.By careful observations made immediately after the peak stress, they estimated the density of cracks in the microstructure and revealed that micro-cracks were visible up to 10 MPa confining pressure.Therefore, microcracking is also expected to have occurred in Qi et al. 2017 [30] experiments.Strong stress concentration occurs at crack tips, and Chauve et al. 2017 [73] have recently shown that, for ice, this stress concentration can be released through DRX.The R 3 iCe formulation does not account for microcracking or the effect of confining pressure.Therefore, the enhancement of recrystallization due to the plastic energy available around the crack tips needs to be reproduced by increasing the recrystallization rate (Γ R X ).
Recent models of texture evolution during dynamic recrystallization have been proposed by Rathmann et al. 2021 [54] and Richard et al. 2021,2022 [47,75], aiming to account for 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 with that in the CTI-RX formulation.The deformability term can be obtained from Equation ( 22) by assuming homogeneous stresses (Sachs model).The diffusive term, which is necessary to reproduce the texture dispersion effect owing to recrystallization, is spontaneously obtained in full-field approaches such as R 3 iCe, due to the heterogeneous stress field.Although these models accurately predict steady-state textures for different strain geometries, the predictions for transient texture evolution are less accurate than that from R 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 [47,54,75] formulations is their numerical efficiency, which makes them suitable for large-scale modeling (glaciers, ice sheets), provided that transient behaviors are not of concern.In such a context, the R 3 iCe model could be a good candidate for providing constraints to develop parameterizations for the temporal evolution, destined to inexpensive formulations.
The CTI-RX model is not limited to the loading conditions prescribed here or to polycrystalline ice.Its formulation is generic and can be extended to other polycrystalline materials, provided that their single-crystal viscoplastic behavior can be adequately described using the continuous transverse isotropic law defined in Equation (17).For instance, the versatility of the CTI-RX model could enable the analysis of the mechanical behavior of materials such as e.g.magnesium and quartz, which also show strong viscoplastic anisotropy with hexagonal symmetry.

Conclusions
In this paper, we present a new formulation, CTI-RX, to model the effect of dynamic recrystallization on the texture evolution of polycrystalline materials, provided that their viscoplastic behavior can be described using a continuous transverse isotropic (CTI) law.This formulation is validated in the context of polycrystalline ice.
The integration of an orientation attractor, denoted as c 0 along with an anisotropic flow law in full-field resolution allows for the accurate replication of texture evolution during dynamic recrystallization in tertiary creep under compression, tension, and simple shear.For the latter, an accurate reproduction of the two-maxima texture and its evolution with strain is obtained.This is made possible by defining c 0 such as to maximize the local resolved shear stress, thus providing a physically based formulation for recrystallization-induced c-axes rotation.This formulation is as simple as possible, given our knowledge of the mechanisms.
The R 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 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.
Accurate reproduction of the textures in the R 3 iCe model leads to a good prediction of the mechanical softening associated with dynamic recrystallization.This confirms that textureinduced viscoplastic anisotropy may explain most of the mechanical softening observed during tertiary creep and suggests that recovery and grain boundary migration play a secondary role in this softening.
Under conditions of high stress, where confining pressure is required in the laboratory to prevent failure, we propose that local fracturing at the crystal scale likely enhances dynamic recrystallization processes and accelerates texture evolution with strain.This is consistent with the lower Γ R X required to fit the texture measured during the confined compression experiments [30] relative to that used to fit the low stress unconfined compression creep experiments [29], despite the higher experimental temperature of the latter (−5 • C vs. −10 • C, cf.Table 3).These results provide an additional illustration of the ability of the R 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 for direct implementation in large-scale flow modeling frameworks.The model is nonetheless highly valuable for constraining the parameterization required to properly account for texture evolution and its impact on the mechanical response in complex or changing boundary conditions such as those prevailing at the bottom of deep ice cores, ice streams, and glaciers.
(2) δ 3 is given by applying stress g D • This gives: • and using δ 1 and δ 3 as defined above:

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

:
Double contracted product A : B A i j B i j ∥.∥ Norm ∥a∥ a 2

Figure 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 .
Figure 3. Schematic representation of the rectangular cubic simulation domain, with surface and point boundaries labeled.

Figure 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 a hexahedral mesh of 2740 elements.

Figure 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 .
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 the results obtained for simulations where grains are defined by multiple mesh elements (grains microstructure) (Figure 4 (a-b)).The full line shows a R 3 iCe simulation where a single orientation is attributed per element (Figure 4 (c)).
with N being the number of orientations), (2) pole figures representing the c-axes distribution at various strains, and (3) the evolution with strain of the texture in a section extracted from the pole figure.These outputs are compared with the pole figures obtained during the experiments at three selected strain values.

Figure 7 .
Figure 7. Distributions of the principal stress tensor directions v 1 and v 3 for uniaxial compression on (a) CraFT-EVP simulation (b) R 3 iCe simulation with grains (c) R 3 iCe simulation with one orientation per mesh element.

Figure 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[13]; b) evolution with strain of the macroscopic strain rate; c) evolution of the intensity of the texture with the colatitude θ over a section defined on the Z plane (dashed blue lines in d)).Data at 22%, 40% and 60% macrostrain are shown.The white dashed lines define the angle at ±45 • from the compression axis Y ; d) textures presented as pole figures for 22% (i ), 40% (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 .
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[29]; b) simulated macroscopic strain rate (black) and corrected simulated macroscopic strain rate (blue) to match the constant force boundary conditions instead of the constant stress one used in the simulation.EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude θ over a section defined on the Z plane (dashed blue lines in (d)).Data at 7%, 12% and 17.8% macrostrain are shown.The white dashed lines define the angle at ±45 • from the compression axis Y ; d) simulated textures (pole figures) for 7% (i ), 12% (i i ), and 17.8% (i i i ) macroscopic strains; e) measured textures for the same macroscopic strains[29].All figures presenting texture data are color-coded using the same color bar (scales in multiples of a uniform distribution).

Figure 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 center of the pole figure shows the direction of the tension axis.From Jacka et al. 1984 [27].

Figure 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 strains are shown.The white dashed lines define the angle at ±45 • from the compression axis Y ; d) simulated textures (represented as pole figures) for 5% (i ), 13% (i i ), and 15% (i i i ) axial macrostrain.All figures presenting texture data are color-coded using the same color bar (scales in multiples of a uniform distribution).

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

Figure 13 .
Figure 13.Simulations for a 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[31].a) Evolution with strain of the eigenvalues of the second-order orientation tensor (a 1 > a 2 > a 3 ).The blue dots correspond to the experimental values from Journaux et al. 2019[31]; b) simulated macroscopic strain rate.EF stands for Enhancement Factor; c) evolution with strain of the intensity of the texture with the colatitude θ 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[31].All figures presenting texture data are color-coded using the same colorbar (scales in multiples of a uniform distribution).

Figure 14 .
Figure 14.Simulation of the 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).

-•-
) = 0, (MS + SM)D = S (A10)The constitutive equation gives g D 13 = (2δ 1 + δ 3 ) δ 2 identification is less trivial and presented in more detail.•First, uniaxial stress σ is applied along the c-axis.The stress tensor is g σ = ( using the constitutive equation, we obtain for g D 33 : Then a uniaxial stress σ is applied normal to the c-axis, for instance in the direction 1.The stress tensor is using the constitutive equation, we obtain for g D 11 :

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

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

Table 2 .
Parameters of the CTI law for the ice single crystal.

Table 3 .
Values of the numerical and physical parameters used in the model R 3 iCe.
Finally from the grain behavior, we can write that γ g D 33 = g D 11 :