Can changes in ice-sheet flow be inferred from crystallographic preferred orientations?

25 Creep due to ice flow is generally thought to be the main cause for the formation of crystallographic preferred orientations (CPOs) in polycrystalline anisotropic ice. However, linking the development of CPOs to the ice flow history requires a proper understanding of the ice aggregate's microstructural response to flow transitions. In this contribution the influence of ice deformation history on the CPO development is investigated by means of 30 full-field numerical simulations at the microscale. We simulate the CPO evolution of polycrystalline ice under combinations of two consecutive deformation events up to high strain, using the code VPFFT/ELLE. A volume of ice is first deformed under co-axial https://doi.org/10.5194/tc-2021-224 Preprint. Discussion started: 17 September 2021 c © Author(s) 2021. CC BY 4.0 License.

by providing a study of the multi-stage deformation of ice samples and the resulting CPOs for different settings. We present a series of numerical simulations of polycrystalline ice that is affected by two consecutive deformation events. We choose combinations of subsequent deformation kinematics that may be expected to occur in ice sheets. We analyse the CPO 105 developed during these deformation events in order to determine the kinematic conditions for the preservation, modification or destruction of CPO's in ice. 110 We analyse different examples of flow changes that represent relevant and/or common deformation regimes in ice sheets. Ice flows in all directions from the accumulation zone towards the edges of the ice mass driven by the gravity force. Due to precipitation and accumulation, ice is gradually compressed ranging from vertical axial compression in domes, to plane strain vertical compression in ridges (zone I in Fig. 1a). Evidence for extensional 115 flow transverse to the ridge has been found in ridges (Fig. 1b), as in the NorthGRIP (Wang et al., 2002;Faria et al., 2014) and EDML ridges (Weikusat et al., 2017).  (Wang et al., 2002). Flow turns progressively more non-coaxial with depth and distance from the dome, where strongly non-coaxial deformation dominates close to the bedrock (indicated by zone II). The dominant kinematic regime at the ice stream is extension along the flow direction, where simple shear deformation with a vertical shear plane can occur in shear margins of the stream (indicated 125 by zone III). At the ice shelf, the predominant regime assumed in this contribution is simple shear with vertical shear plane (zone IV) (Lutz et al., 2020). https://doi.org/10.5194/tc-2021-224 Preprint. Discussion started: 17 September 2021 c Author(s) 2021. CC BY 4.0 License. At depth and away from the centre of the ice sheet (zone II in Fig. 1a), ice experiences a vertical gradient in velocity, resulting in simple shear deformation parallel to the bedrock 130 (Hudleston, 2015). Inside a glacier, ice stream or in a flank flow (Voigt, 2017), flow acceleration may dominate, resulting in extension along the flow direction (zone III in Fig.   1a). In a shear margin of glaciers and ice streams, as well as in some ice shelves, ice experiences a gradient in velocity in the lateral direction perpendicular to that velocity, resulting in simple shear deformation with a vertical shear plane (zone III and IVa in Fig.1a-135 b) (Young et al., 2002;LeDoux et al., 2017;Lutz et al., 2020). Notice that, in the presented configuration the xyz reference frame is defined with y vertical, x horizontal and parallel to flow and z normal to xy.

Numerical simulation of ice fabrics and postprocessing
Ice polycrystalline viscoplastic deformation was simulated using the Fast Fourier Transform algorithm (VPFFT; Lebensohn and Rollett, 2020), within the numerical open-source platform 145 ELLE (http://www.elle.ws; Bons et al., 2008). ELLE has been successfully used for the simulation of a variety of studies of rock microstructure evolution during deformation and metamorphism (see Piazolo et al., 2019). The coupling of the full-field VPFTT code and ELLE allows simulating deformation of a polycrystalline aggregate by dislocation glide up to high strains (>10) as often found in nature in ice and rocks. The reader is referred to 150 Lebensohn and Rollett (2020) for a detailed description of the VPFFT approach, and to Griera et al. (2013) and Llorens et al. (2016a) for a detailed description of the coupling between VPFFT and ELLE. The local mechanical response of a nonlinear heterogeneous material can be calculated as a convolution integral of the Green functions associated with a linear homogenous medium and a polarisation field. The VPFFT formulation is used to 155 transform the polarisation field that contains all the information on the heterogeneity and non-linearity of the material's behaviour into Fourier space. By the conversion of the real space convolution integrals to simple products in the Fourier space, the mechanical fields are calculated, and the convolution product is transformed back to real space. In this full-field approach, the semi two-dimensional data structure is discretised by a regular mesh of 160 256x256 unconnected nodes (unodes) or Fourier points (Fig.2). The strain rate and stress field under compatibility and equilibrium constraints related by the constitutive equation (1) https://doi.org/10.5194/tc-2021-224 Preprint. Discussion started: 17 September 2021 c Author(s) 2021. CC BY 4.0 License.
is obtained by an iteratively solving for the flow law for every unode (u): Where ̇ is the strain rate, ′ is the deviatoric stress, is the symmetric Schmid tensor, ̇ the shear strain rate, is the critical resolved shear stress defined for the slip system (Ns), 0 is a reference strain rate and n is the rate sensitivity exponent. In all three slip systems, the shear strain rate is assumed to be related to the deviatoric stress 170 by a stress exponent of n=3, in accordance with the "Glen's law" (Haefeli,1961). After convergence, the VPFFT code calculates the associated lattice rotations from the velocity gradient and stress fields. The actual stress exponent in ice 1h may actually be closer to n=4 (Bons et al., 2018). The first tests show that raising n to 4 has little effect on the CPO, and in order to be coherent with previous ice CPO simulations (Llorens et al., 2016a, 2016b, 175 Steinbach et al., 2016Qi et al., 2018;Piazolo et al., 2019) all simulations presented here are carried out with the more commonly used n=3. Uniaxial shortening of a single crystal by glide along the non-basal planes only requires about 60 times higher stress than when the crystal can deform entirely by basal-plane glide (Duval et al., 1983). We therefore set 60 times lower for basal glide than for prismatic and pyramidal glide. 180 Previous numerical studies have shown that about the same single-maximum CPO develops at low temperature and a high strain rate conditions, regardless of the set rate of dynamic recrystallisation (Llorens et al., 2016a;2016b). As this contribution aims to study the CPO response to a change in the deformation regime, the numerical procedures that simulate dynamic recrystallisation are deliberately not incorporated in this study. Every unode 185 represents a crystallite where the crystal orientation, strain rate, dislocation density and local stress is stored (Fig.2b). All unodes within a grain have the same initial crystal orientation ( Fig.2a). The crystal orientation is defined by the three Euler angles. We use the three Euler  Fourier points (also termed unconnected nodes or unodes), and (c) a boundary node (bnodes) layer that define grains as polygons. 205 Deformation was applied to the microstructure by incremental steps of a shear strain of γ=0.04 or, alternatively, 2% of shortening or extension (Table 1), both equivalent to 0.02 of natural strain, defined as ln(Lf/Li) (with Lf and Li being the final and initial length of a line in the direction of maximum extension, respectively). Each examples considered in this contribution consists of two deformation regimes that are applied in succession (see Table 2). 210 We considered four different model series to simulate flow transitions according to combinations of deformation regimes described in Figure 1. Table 1. Velocity gradient tensors V considered for the different regimes (x 10 -8 s -1 ) and corresponding zones in and outflow in all horizontal directions (zone I in Fig. 1). As ice flows away from domes and gets buried deeply it enters the zone dominated by simple shear parallel to the bedrock (zone II in Fig. 1). To simulate this transition, we carried out a series of simulation with first vertical uniaxial compression parallel to y (V1), followed by dextral simple shear on horizontal plane xz (V2) ( Table 1 and 2). 225 On ridges away from ice domes, ice flows away and stretches in the direction perpendicular to the ridge (Wang et al., 2002;Faria et al., 2014) (zone I in Fig. 1). Here we consider the end-member case that this leads to uniaxial extension in the flow direction. Similar to Series A, the ice is then assumed to be buried and to enter zone II (Fig. 1) where flow is dominated 235 by bedrock-parallel simple shear, again in the direction perpendicular to the ridge (zone II in Fig. 1). For this case, we considered simulations in which a polycrystalline ice aggregate is first deformed by V3, horizontal uniaxial extension parallel to x, followed by V2, dextral simple shear on horizontal plane xz (Table 1 and 2).

Series C: Ice flowing from an ice dome to an ice flank or stream
In this configuration we assume that ice is first gradually flattened by vertical uniaxial compression parallel to y in the ice dome (V1; zone I in Fig. 1) and subsequently enters a zone in which the flow accelerates in, for example, an ice stream, leading to uniaxial 245 extension in the flow direction (V3; zone III in Fig. 1). Series C simulations include a set of simulations where V3 deformation starts at different strains in the first stage with V1.

Series D: Ice flowing from an ice-stream or glacier to an ice shelf or shear margin
250 Dominantly simple shear deformation with a vertical shear plane takes place in the margins of ice streams (Hudleston, 2015) and can be found in some ice shelves (Young et al., 2002;LeDoux et al., 2017;Lutz et al., 2020) (zone IV in Fig.1). Ice affected by this shearing may have experienced different deformation types. Here we consider the case of uniaxial extension in the flow direction (zone III in Fig.1), which would represent the ice above the 255 deepest layers, that thus did not experience of bedrock-parallel simple shear before.
A series of simulations were run to simulate the scenarios involved in this scheme: first uniaxial extension in the x direction (V3; zone III in Fig. 1), and subsequently dextral simple shear on vertical plane xy (V4; zone IV in Fig. 1).

285
At the beginning of the second regime the CPO slightly rotates antithetically (i.e., opposite to the imposed shear sense) (see ε2=1 in Fig. 3b), while the CPO intensity (Fig. 3d) and girdle component decrease (Fig. 3c). After that the CPO intensity (point maximum) gradually intensifies again and the c-axes align at about 10º to the normal to the shear plane (see ε2=3 in Fig. 3b and Fig. 3c). After the two deformation events the resulting CPO resembles the one 290 resulting from simple shear deformation only, in both CPO intensity and symmetry (green and blue markers in Fig. 3d and Fig. 3c). This implies that the previous uniaxial flattening regime is not recognizable after only a moderate amount (ε1>1) of subsequent simple shear, as may be expected as the change in CPO only involves a small rotation and intensification of the c-axes point maximum.

Series B: Ice flowing from the centre of the ridge to deep lateral zones
In series B, the initially random distributed c-axes develop a CPO characterised by a strong vertical girdle, in coherence with the horizontal extension applied at the first deformation 300 regime (see ε1=0.92 in Fig. 4a). From the final step onwards, the aggregate is deformed under horizontal top-to-the-right simple shear up to a natural strain of ε2=4 (Fig. 4b). C-axes reorient following the new flow, destroying the vertical girdle fabric (Fig. 4c), and forming a broad single point maximum almost normal to the shear plane (see ε2=4 in Fig. 4b). Although the final CPO symmetry is coherent with simple shear deformation, its shape after a strain of 305 ε2=4 still differs from that of the previous case (series A) or that of simple shear only (see the last step for the second regime in Figs. 3a and 4a). At the highest modelled strain (ε2=4) the CPO intensity is clearly lower than that developed under solely simple shear deformation ( Fig. 4d). Changing a girdle to point maximum requires a significantly higher strain than the formation of a point maximum only (green versus blue graph in Fig.4d).  320 The microstructure in series C first develops a strong CPO (Fig. 5b) with a point maximum in the y-direction of the axial compression applied during the first deformation regime (see first row in Fig. 5a). In this series the second deformation regime was applied at different steps of natural strain of the first regime (at ε1=0.2, 0.4, 0.8 and 1.2, respectively). The intensity of the initial deformation regime has a notable influence on the CPO developed during the second 325 regime. In the absence of first-regime deformation (ε1=0.0), where the microstructure only experiences uniaxial extension in the x-direction, a vertical girdle fabric develops (see first column in Fig. 5a and green marker Fig. 5b). Regardless of the step at which the second regime is applied, the CPO tends to evolve towards a girdle fabric (see evolution of all series in Fig. 5a-b), which reduces slightly the CPO intensity (Fig. 5c). If the second deformation 330 regime is applied when the microstructure has only been slightly affected by the first regime (ε1 < 0.4), the c-axes tend to rotate forming a final girdle fabric but with a recognizable point maximum component, resulting in an "hourglass shape" (see second and third columns in Fig. 5a). However, if a strong point maximum CPO has been developed during the first deformation regime (ε1 > 0.8), the second regime cannot modify the inherited c-axis 335 orientation at a natural strain of ε2 = 1.2 (see fifth column in Fig. 5a), and the final CPO continues being dominated by a strong point maximum (Fig. 5b). The effect of the second regime is observable in the a-axis distribution, because regardless of the intensity of the first regime, the extension along the x direction produces three a-axis maxima (see last row in Fig.   5a). The result suggests that converting a point maximum into a girdle takes even more strain   350 In series D a vertical girdle fabric develops during the first deformation regime (ε1 = 0.92 in Fig. 6a), which then rapidly evolves towards a point maximum due to simple shear along a vertical shear plane imposed during the second flow regime (Fig. 6b). However, after significant strain there is still a girdle component.

365
Experiments and simulations with an initial random distribution of crystallographic orientations predict a quick alignment of the crystallographic axes according to the imposed deformation conditions (Azuma and Higashi, 1984;Qi et al., 2017, Llorens et al., 2016aLlorens et al., 2016b). The types and orientations of CPO's are best described in terms of the principal axes of the finite-strain ellipsoid (FSE), also known as finite stretching axes (FSA's) 370 that are typically labelled X, Y and Z from longest to shortest (Fig.7) (Passchier, 1990). These axes are parallel to the instantaneous stretching axes (ISA's) in the case of coaxial deformation, but not in the case of simple shear. C-axes tend to align themselves parallel to the shortening FSA's. In the case of uniaxial compression (V1 in Table 1; X=Y>Z) this leads to the development of a strong point maximum. A girdle forms in uniaxial extension (V3 in 375 Table 1; X>Y=Z). In simple shear (V2 and V4) the orientation of the Z-axis rotates from initially 45° to 0° to the shear plane at infinite strain, while the X-axis remains in the shear plane. A point-maximum CPO is expected and observed parallel to Z in the plane between the instantaneous shortening axis and the normal to the shear plane. a-axes {11-20} align themselves perpendicular to the c-axes in the plane of the two least compressive FSA's (i.e. 380 the X and Y-axes of the FSE) (see first column in Fig. 7).   (Azuma et al., 1999), Talos Dome (Montagnat et al., 2012), and Greenland ice cores GISP2 (Gow et al., 1997) and GRIP (Thorsteinsson et al., 1997). They also match 395 laboratory experiments carried out at low temperatures and high strain rates (Craw et al., 2018;Fan et al., 2020). The vertical girdle fabric predicted in our simulations under uniaxial extension (first-regime deformation in series B and D) is also observed in Antarctic and Greenland ice cores, such as the Vostok ice core (Lipenkov et al., 198;Voigt, 2017), EPICA DML (Weikusat et al., 2017), NEEM (Montagnat et al., 2014a) or in the Styx Glacier in 400 Northern Victoria Land (Kim et al., 2020). As described in the NorthGRIP Greenland ice core, the vertical girdle is interpreted as evidence for extension transverse to the ridge, where the vertical girdle plane is oriented perpendicular to the axis of horizontal extension (Wang et al., 2002;Faria et al., 2014;Weikusat et al., 2017). Greenland, where an approximately vertical single maximum CPO is found (e.g., Faria et al., 2014;Montagnat et al., 2014a;Weikusat et al., 2017).
Simulations predict that the land-based ice CPO is destroyed at the ice shelf, assuming the example of shearing on the vertical plane as the dominant deformation regime, as observed in the Antarctic Amery (Young et al., 2002) or Western Ross ice shelfs (LeDoux et al., 2017;420 Lutz et al., 2020). In this case, although the inherited vertical girdle CPO is not completely destroyed at the end of the simulation (ε =2), it is progressively overprinted by the shearing on the vertical plane, (series D in Fig. 7).
When an inherited single vertical maximum is affected by extension along the flow, similar to ice flowing in an ice stream that previously has been accumulated in an ice dome, the CPO 425 is also progressively overprinted, but the required strain is clearly much higher than that for the other cases (ε >1.2) (series C Fig. 7). Simulations predict that ice moving into a shear margin or potentially developing a shear margin during the initiation of an ice stream, will basically loose its original CPO as shearing on the vertical plane overprints the inherited fabric. 430 As observed in the simulations, the effectivity of the second flow regime on the reorientation of the CPO depends not only on the strength of the inherited CPO, but also on the relationship between the crystal orientation and the Finite Strain Ellipsoid (FSE) of the new imposed regime. When the inherited a-axis preferred orientation and the elongation axis of the new FSE are parallel (see series C and D in Fig. 7) c-axis orientations are more slowly (glide direction) along the a-axes.
The flow behaviour from the pole figure observed in Series C at ε1 = 1.2 and ε2 = 1.2 (see lower right corner in Fig. 5a) could easily be misinterpreted as axial compression in the y-440 direction (Fig. 3b), while the true current flow is axial extension in the x-direction.
The effectivity of the second flow regime is also observable in the relative activity of slip systems (Fig.8). In all cases, the new imposed regime produces a remarkable increment of basal slip, while the pyramidal slip system activity is reduced (see series A, B and D in  All these results demonstrate that an inherited CPO can change entirely within strains of ε <1.2 with a range of transition fabrics (Fig.9). Therefore, these simulations do not support the assumption that an inherited CPO from a previous regime is not destroyed by the current flow (Smith et al., 2017). This process depends on both the intensity of the inherited CPO, and its experiments of natural samples with a pre-existing CPO, where the application of a stress field in a non-favourable orientation with respect to the inherited CPO destroys it (Jun andJacka, 1998, Craw et al., 2018).
One particular case is an inherited strong CPO, where a-axes are optimally oriented with 460 respect to the second flow regime, because the FSE elongation axis from the first and second regime are in the same plane (see FSA in series C and D in Fig.8). In this situation the second regime does not produce a rapid re-orientation of the a-axes with respect to the elongation axis, the basal activity is not enhanced and therefore the CPO is only slightly modified (series C and D in Fig.9). In this case the required finite strain for overprinting the initial CPO is 465 clearly much higher than that for the other cases. 2. Depending on the sequence of deformation regimes, an inherited CPO can be completely overprinted by the later deformation event. An inherited CPO can change entirely within natural strain lower than 1.2, with a range of transitional fabrics. This process depends both 485 on the intensity of the inherited CPO and its orientation with respect to the new stress field.
3. More specifically, when the inherited a-axis preferred orientation and the elongation axis of the finite strain ellipse of the flow are parallel, lattice orientation only needs to be reoriented slightly. The required finite strain for overprinting the initial CPO in this configuration is much higher (ε >1.2) than that for the other cases. This is the case of ice 490 flowing from an ice dome to an ice stream. This situation could lead to a misinterpretation of the second flow regime from observed c-axis preferred orientation.
4. According to these results, ice flow interpretations from observed CPOs can be reliable, but must be carried out with caution in areas with complex (multi-stage) deformation histories. 495