## Abstract

The flow of glaciers and polar ice sheets is controlled by the highly anisotropic rheology of ice crystals that have hexagonal symmetry (ice lh). To improve our knowledge of ice sheet dynamics, it is necessary to understand how dynamic recrystallization (DRX) controls ice microstructures and rheology at different boundary conditions that range from pure shear flattening at the top to simple shear near the base of the sheets. We present a series of two-dimensional numerical simulations that couple ice deformation with DRX of various intensities, paying special attention to the effect of boundary conditions. The simulations show how similar orientations of *c*-axis maxima with respect to the finite deformation direction develop regardless of the amount of DRX and applied boundary conditions. In pure shear this direction is parallel to the maximum compressional stress, while it rotates towards the shear direction in simple shear. This leads to strain hardening and increased activity of non-basal slip systems in pure shear and to strain softening in simple shear. Therefore, it is expected that ice is effectively weaker in the lower parts of the ice sheets than in the upper parts. Strain-rate localization occurs in all simulations, especially in simple shear cases. Recrystallization suppresses localization, which necessitates the activation of hard, non-basal slip systems.

This article is part of the themed issue ‘Microdynamics of ice’.

## 1. Introduction

Ice is one of the most common minerals found at the Earth's surface. Most of it is concentrated in polar ice sheets, formed originally from precipitation of snow. Owing to its ductile behaviour, ice flows towards the ice sheet margins driven by gravity. While this behaviour is commonly described by continuum-mechanical approaches, its underlying processes are of crystallographic nature, where deformation takes place by dislocation creep [1]. At terrestrial conditions, single ice crystals have a hexagonal symmetry (ice lh) and deform plastically by preferential glide of dislocations along the basal planes of the hexagonal lattice, where the slip of dislocations is almost isotropic. The stress required to activate dislocation glide along non-basal slip systems is about 60 times higher than that for the basal one [2], when subjected to similar strain rates.

Natural ice usually occurs in the form of polycrystalline aggregates, which develop a *c*-axis preferred orientation during deformation. This process is facilitated by recrystallization of the ice grains, as polar ice is always at high homologous temperatures (i.e. near to the melting point). The microstructural evolution of the mineral aggregate during deformation leads to progressively anisotropic macroscopic behaviour (e.g. [3]). The development of a fabric or a lattice preferred orientation (LPO) due to deformation influences the response of ice layers to imposed stresses [4], deforming an order of magnitude faster than an equivalent isotropic polycrystal in simple shear deformation. We define LPO here as the statistically preferred orientation of the crystalline lattices of a population of grains [5]. Activation of recrystallization processes during plastic deformation (dynamic recrystallization, DRX) can enhance ice flow, modifying the microstructure (defined as the collection of all microscopic deformation-related structures and the orientation stereology of polycrystals [5]) and LPO [6,7], or activating non-basal slip systems [8]. It is, therefore, critical to understand how DRX influences ice microstructures at different boundary conditions, because microstructural changes are related to climatic transitions in the palaeoclimate record stored in ice sheets [9].

Natural ice deformation fabrics and microstructures can be studied with different techniques. Firstly, the analysis of ice microstructures measured from ice core samples provides an insight into the nature of strongly deformed natural ice. There are abundant data from the upper 100 m of ice cores, where firn transitions into ice. However, only a limited number of datasets from deep ice cores exist. Very few campaigns have been able to acquire deep ice samples, because they are expensive and logistically demanding (see [5] for an overview). Furthermore, coring sites are preferably located in very special strain environments, mainly in zones of ice divides and domes. This is because deformation in these zones mostly takes place under vertical compression (e.g. pure shear or axial compression) as the ice is buried under annual snow layers. These locations are specifically chosen to get relatively undisturbed layers to extract the clearest possible palaeoclimate signals. Therefore, the deformation regimes that are subjected to high shear strain rates, which are found at the flanks or outlets of the ice sheets, are usually not prioritized for deep drilling campaigns, as this regime gives rise to folding of ice layers [9–12]. However, it is critical that we understand how ice behaves in all the parts of ice sheets and not only in ice divides and domes. Only few ice coring projects have focused to date on recovering samples with the aim of studying ice flow properties [13,14]. The EGRIP campaign will soon recover ice samples from the onset region of an ice stream (EGRIP website: http://eastgrip.org/).

Experiments are another widely used technique to learn about ice rheology during deformation and recrystallization. Samples of natural or artificial ice are deformed under controlled strain or stress boundary conditions and at a certain temperature in order to understand the processes that govern deformation and the resulting microstructures. These methods have significantly contributed to our understanding of ice rheology [15–17]. However, experiments are limited to certain geometries and usually do not achieve high finite strains. Another disadvantage here is that experimental samples are subjected to much higher strain rates than natural ice, especially in the case of slowly flowing polar ice sheets ( [13,18]). As the viscosity of ice depends on the applied stress and this relation is highly nonlinear, other techniques should be used in order to complement our understanding of ice rheology [19].

Numerical simulations provide a tool to fill in the gaps left by the study of natural and experimental ice deformation, because simulations are not limited by scale, stress or strain rate. Moreover, numerical models allow the systematic investigation of different aspects of DRX on ice flow, based on the underlying physical processes [20]. In this way, it is possible to link processes on the dislocation scale to the behaviour of ice at the scale of the polycrystal. Additionally, by evaluating the relation between stress and strain rate from the bulk behaviour at the mesoscopic (centimetre to metre) scale, we can also establish links with the macroscopic scale. Furthermore, modelling allows the comparison of ice rheology under different boundary conditions, such as pure- and simple-shear deformation. Thus, microstructural modelling can be a valuable addition to the analysis of ice microstructures in natural ice samples and laboratory experiments.

Despite the importance of DRX on the behaviour of natural ice, theoretical and computational studies have been restricted mostly to the simulation of (i) ice deformation at the scale of single crystals [21] or (ii) deformation of polycrystalline materials using finite-element methods, whereby isotropic rheological properties are used to model crystal plasticity [22,23].

The full-field formulation based on fast Fourier transforms (FFT) [24] has proved to be a successful alternative modelling tool to simulate plastic deformation of polycrystals [25,26], especially when coupled with algorithms for DRX processes [8,27]. Viscoplastic full-field numerical modelling (in short VPFFT modelling) of the evolution of a polycrystalline aggregate under controlled strain boundary conditions can be applied to investigate the relative importance of deformation mechanisms and recrystallization processes on the development of ice fabrics and flow (e.g. [20]). Here we use the VPFFT formulation inside the modelling platform ELLE [28,29], which has been widely used to model processes in rocks and has recently been adapted to experiments on polycrystals with the characteristics of natural ice. The aim of this paper is to present the state of the art of the VPFFT/ELLE modelling of ice, paying particular attention to the effect of vorticity boundary conditions (i.e. pure versus simple shear) and DRX.

## 2. Numerical procedure

### (a) ELLE numerical platform

The software platform ELLE [28,29] was used to perform the simulations presented in this contribution. The ELLE software is open source and can be downloaded from the online repository (http://www.elle.ws). ELLE processes reproduce the evolution of microstructures during deformation and metamorphism, such as grain growth [30–32], DRX [8,33,34], folding [12,27,35,36] and rotation of rigid objects in anisotropic materials [37,38].

In this contribution, we compare the microstructural evolution of ice lh during viscoplastic deformation and recrystallization. The numerical approach is based on the coupling of (i) a full-field viscoplastic code that uses a FFT approach (VPFFT [25,39]) to calculate the viscoplastic deformation of the polycrystalline aggregate and (ii) the microstructural modelling platform ELLE for the simulation of DRX.

Two main layers define the data structure of the models: (i) a contiguous set of polygons (figure 1*a*,*b*) that are used to define grains and that are themselves defined by straight boundary segments joined at boundary nodes (*bnodes*; figure 1*d*) and (ii) a high-resolution grid of unconnected nodes (*unodes*) that store physical properties within grains, such as the local lattice orientation, defined by three Euler angles, and the dislocation density (figure 1*e*). These *unodes* also serve as Fourier points to simulate crystallites for the viscoplastic deformation calculations [25]. The data structure of the VPFFT and ELLE codes is fully wrapping: grains touching one side of the model continue on the other side. Two different sets of initial microstructures are used to investigate three series of simulations (table 1): a single-phase ice polycrystalline aggregate with a cell size of *L* × *L* to be deformed in simple shear (figure 1*a*) and an *L* × 2*L* cell to be deformed in pure shear vertical shortening (figure 1*b*). A grid density of 256 × 256 *unodes* was used to define a unit cell, resulting in a total of 65 536 discrete *unodes*. Grains initially have a random lattice orientation (figure 1*a*–*c*), representing a bulk isotropic material.

In our approach, each small time increment (Δ*t*) is represented by a sequential loop of processes. For all the simulations presented, an experimental run (figure 1*c*) consists of iterative applications of an incremental 0.02 natural strain, defined as ln(*Lf*/*Li*), equivalent to 2% of shortening of vertical coaxial compression or an incremental shear strain of Δ*γ* = 0.04 in simple shear simulations. The deformation step is followed by a subloop of DRX, where grain boundary migration (GBM), recovery and polygonization processes are activated. Since the mobility parameter of recrystallization processes and time step are the same for all simulations, we model deformation at different strain rates by varying the ratio between the viscoplastic deformation step and the DRX step, activating the DRX subloop *N*_{DRX} = 0, 1, 10, 20 or 25 times per deformation time step (see values in table 1). When DRX is not active it is not possible to calculate the strain rate, because there is no time step defined in the model. When DRX is active, *N*_{DRX} is inversely proportional to the strain rate. Using the parameters defined further below, the *N*_{DRX} correspond to strain rates ranging from 3.17 × 10^{−11} down to 1.27 × 10^{−12} s^{−1}.

The first series of simulations (labelled by PSH *N*_{DRX}) models DRX, including GBM and recovery, of ice under pure-shear vertical shortening up to a natural strain of 1.2. This series includes a simulation (PSH 20_N) that incorporates polygonization. The second series (SSH *N*_{DRX}) simulates DRX of ice*,* including GBM and recovery, under dextral simple shear conditions up to a natural strain of 1.2. This series includes simulations (SSH 1_N and SSH 10_N) that incorporate polygonization. The simulation set-up is listed in table 1. Lattice orientation data, as well as misorientation and activity of slip systems fields, were plotted using the texture analysis software MTEX [40].

### (b) Fast Fourier transform: viscoplastic deformation by dislocation glide

VPFFT provides a strain rate and stress field that minimizes the average local work rate under compatibility and equilibrium constraints [25,39]. VPFFT is a full-field crystal plasticity formulation that explicitly resolves velocity and stress fields, assuming that deformation is accommodated by dislocation glide only. The constitutive relation between the strain rate and the deviatoric stress *σ*′(*x*) at a point *x* of the Fourier grid is given by
2.1
where the sum runs over all (*N _{s}*) slip systems

*s*in the crystal,

*m*

^{s}is the symmetric Schmid tensor,

*τ*

^{s}is the resistance to slip or critical resolved shear stress, is the shear strain rate, is a reference strain rate and

*n*is the rate sensitivity or stress exponent in Glen's law [41]. To simulate the mechanical properties of a polycrystalline aggregate of ice lh crystals, a mineral with hexagonal symmetry and deformation accommodated by glide of dislocations along basal and non-basal slip systems (i.e. prismatic and pyramidal) is assumed. The degree of anisotropy (

*A*) is defined by the ratio between the critical resolved shear stresses

*τ*

^{s}between the non-basal and the basal slip systems.

*A*is approximately 60 according to Duval

*et al.*[2], but the basal slip system is assumed to be 20 times weaker than the pyramidal and prismatic ones (

*A*= 20) in all the simulations presented here. This value of

*A*is a compromise to avoid too long a calculation time, while still maintaining an acceptable accuracy of the result [8]. The rate sensitivity exponent is

*n*= 3 for all the slip systems. The resulting micromechanical fields from VPFFT are used to estimate the geometrically necessary dislocation density and deformation-induced lattice rotation, required inputs to simulate intra-crystalline recovery and GBM.

### (c) ELLE: recrystallization processes

After each deformation step, recrystallization is simulated by means of three processes: polygonization, GBM and recovery. The polygonization routine creates new high-angle grain boundaries (HAGB) when a cluster of *unodes* within a grain have a misorientation angle with respect to the parent grain that is higher than a critical pre-defined angle. GBM is driven by the reduction of grain boundary energy and stored strain energy, reproducing the motion or displacement of HAGB [42]. Recovery reduces the intra-granular stored energy in a deformed crystal, simulating annihilation of dislocations and their rearrangement into low-angle subgrain boundaries (LAGB) [43].

#### (i) Polygonization

The polygonization routine modifies the *bnodes* layer through splitting a pre-existing grain by inserting new boundary nodes. This process commences with a routine that identifies subgrains, i.e. clusters of *unodes* with misorientation values above a pre-defined threshold (HAGB more than 15°; figure 2*b*). Once a critical subgrain is detected, the position of its boundary is obtained using a Voronoi tessellation of the cluster of unodes clipping the boundary to the pre-existing grain. A convex hull routine is used to define the smallest convex polygon enclosing the Voronoi points and clipped boundary nodes. The use of Voronoi decomposition allows us to prevent topological errors in the network of grains.

#### (ii) Grain boundary migration

GBM is modelled according to the front tracking approach explained in more detail in [8,33]. GBM is driven by a reduction of surface energy and stored strain energy. The driving stress (Δ*f*) is derived from the internal strain energy (Δ*H*) and the total grain boundary surface energy, which depends on the surface energy (*J*)and the local radius of curvature (*r*) of the grain boundary:
2.2

The velocity (** v**) of the grain boundary is calculated from the driving stress (Δ

*f*) and grain boundary mobility (

*M*) such that the work done equals the change in the local energy state: 2.3 where

**is a unit vector in the direction of movement. The mobility is highly dependent on temperature (**

*n**T*) and it is calculated from the intrinsic mobility (

*M*

_{0}= 0.023 m

^{2}kg

^{−1 }s

^{−1}[44]), the grain boundary diffusion activation energy (

*Q*=

*40 kJ mol*

^{−1}[45]) and the universal gas constant (

*R*): 2.4

The stored energy (Δ*H*) is calculated from the difference in dislocation density across a grain boundary (Δ*ρ*), the shear modulus (*G*, assumed to be isotropic) and the Burgers vector (** b**):
2.5

This routine picks single *bnodes* in a random order and uses four orthogonal small trial moves (parallel to the *x*- and *y*-axes, and of 1/100 of the average distance between *bnodes*) of that *bnode* to determine the direction (** n**) with the maximum driving stress (Δ

*f*) [33]. The

*bnode*is then moved over a small distance,

**·Δ**

*v**t*, in the direction of

**according to equation (2.3). In areas swept by the moving boundary, the stored strain energy (Δ**

*n**H*) is reduced to zero.

#### (iii) Recovery

The recovery module that simulates intra-granular reduction of the stored strain energy is applied after GBM. The recovery routine implemented in ELLE [43] simulates annihilation of dislocations and their rearrangement into LAGB. This process assumes that each *unode* represents a small crystallite or potential subgrain. The rotation rate of the lattice within the crystallite is proportional to the torque applied on the sides of the crystallite. This torque is calculated from the angular mismatch of the *unode* and its neighbours, which results in an effective interface energy [46]. The algorithm, described in detail in [8], is applied to an *unode* and its neighbours that belong to the same grain. The lattice orientation of the reference *unode* is rotated towards the value that results in the maximum reduction in interface energy calculated from trial rotations. This procedure is repeated for each *unode* in random order.

## 3. Results

### (a) Microstructure evolution

At the end of the simulations (natural strain of *ε* = 1.2) differences in grain size and grain shape are observed. They reflect the different strain rates applied and the concomitant differences in the activity of DRX (table 2). The case of ice deformation without DRX is represented by experiments PSH 0 and SSH 0, for pure and simple shear boundary conditions, respectively. In both cases, viscoplastic deformation leads to the development of a grain shape preferred orientation (SPO) by the irregular elongation of grains oriented parallel to the stretching direction (figure 3*a*,*e*). This direction is the horizontal *x-*axis in pure shear models, while in simple shear models it rotates clockwise from 45° with respect to *x* at the beginning of the simulations towards the shear direction that is parallel to the *x*-axis. Deformation localization in the simple shear models produces strongly elongated grains in high-strain bands, which form at low angle with respect to the shear direction (figure 3*a*). In pure shear models conjugate sets of such high-strain bands are observed, but oriented symmetrically with respect to the vertical maximum shortening direction (figure 3*e*). In both cases, the high-strain bands are obliquely oriented with respect to the grain SPO. These bands are additionally characterized by high values of misorientations (i.e. intra-granular heterogeneities calculated as the angular difference in lattice orientation between an *unode* and its neighbours), progressively developing high local misorientations and recognizable subgrains, as can be seen in the misorientation field maps (figure 3*e*).

Different degrees of DRX were simulated in experiments PSH and SSH 1, 10 and 25, where the number refers to *N*_{DRX}, i.e. the number of recrystallization steps per VPFFT steps in the loop. The activation of GBM and recovery produces a general grain size coarsening, with grains becoming more equant with decreasing strain rate (table 2 and figure 3*b*–*d*,*f*–*h*). As a result, the grain SPO qualitatively weakens, but preferred elongation of grains parallel to the stretching direction is still recognizable in models that have up to *N*_{DRX} = 25 steps per deformation increment (figure 3*b*–*d*,*f*–*h*). A reduction of the local misorientation by the arrangement of subgrain domains results in the development of low-angle boundaries and lower intra-granular heterogeneities than in models with *N*_{DRX} = 0. Subgrain size generally tends to increase with *N*_{DRX}. Contrary to what can be observed in models without recrystallization, bands with highly stretched grains do not develop in other simulations, except in models with *N*_{DRX} = 1. However, these bands are small compared with those in models that only include deformation. In the models *N*_{DRX} = 10 and *N*_{DRX} = 25, these bands are not recognizable and grain size distributions are relatively uniform. In general, grain boundaries show straight to smooth geometries. Triple junctions are typically characterized by angles close to 120° with increasing *N*_{DRX}, but large deviations from this angle can still be qualitatively observed in the *N*_{DRX} = 25 cases.

### (b) Lattice orientation evolution

Orientation maps and lattice orientation distributions at a natural strain of *ε* = 1.2 are displayed in figure 4. Differences in lattice orientation are recognizable from colour differences according to inverse pole figures relative to the *y* direction of the sample. All simulations present a final similar LPO with respect to the finite deformation directions, regardless of the applied boundary conditions (vorticity) and amount of DRX. The initial randomly oriented grains (initial microstructure in figure 4) evolved towards an LPO approximately defined by *c*-axes (i.e. {0001} axes in figure 4*a*–*d*) parallel to the maximum shortening direction, or perpendicular to the foliation, while the *a*-axes (i.e. {11–20} axes in figure 4) tend to be parallel to the stretching direction, or parallel to the foliation.

In the case of pure shear deformation, the *c*-axes maxima are oriented sub-parallel to the compression axis *σ*_{1} during all simulation time, displaying a broad single maximum for the PSH 0 case and becoming increasingly concentrated towards the *y*-axis with increasing DRX (see polar figures in figures 4*e*–*h* and 5*e*–*h*). In the case of simple shear deformation, *c*-axes are oriented obliquely to the maximum compression axis and rotate in a synthetic sense (i.e. in the same way as the imposed shear sense) towards the normal to the shear plane (or sample *y*-axis) with progressive deformation. It has to be noted that this condition was not achieved at the end of the simulations (figure 5*a*–*d*). Apart from the orientation of the *c*-axis maximum with respect to *σ*_{1}, differences between pure and simple shear are remarkably subtle (figure 5). The fabric generally tends to be slightly weaker in pure shear models and displays a more marked girdle component than in simple shear cases.

When DRX is active, the LPO develops faster and *c*-axes tend to be slightly inclined with respect to the *xy*-plane in both pure and simple shear simulations (figure 5). The effect of DRX is more pronounced on the distributions of *a*- and *b*-axes, since they lie on a broad great circle in models with low DRX, and this trend is divided into a number of distinct maxima for cases with higher DRX. This result is coherent with the enhancement of the LPO single maximum.

### (c) Activity of the slip systems

In the absence of DRX, the activity of the basal slip systems is 60–70% higher than that of the non-basal slip systems in both pure and simple shear simulations (figure 6*a*). The basal slip system activity is concentrated in the high strain-rate bands, oriented at low angles to the horizontal shear plane in simple shear (figure 7*a*), or in conjugate bands distributed symmetrically around the vertical *y*-direction in pure shear simulations (figure 7*c*). When DRX is active the activity of the non-basal slip systems is increased in both pure and simple shear simulations, especially in pure shear (figure 6*b*). The slip system activity is not distributed in bands because they are destroyed by DRX. Larger grains, produced by GBM, contain subgrains with different amounts of basal (tilt and twist boundaries) and non-basal activity (figure 7*b*,*d*).

### (d) Rheology

All simulations show an increment of 30% or more in differential stress with increasing strain (i.e. strain hardening; figure 8*a*). Differential stress is defined as the difference between the maximum and minimum principal stresses. Significant deviations of the stress evolution between pure and simple shear only start developing after *ε* ≈ 0.5, where differential stress stabilizes in simple shear models (light lines in figure 8*a*) and continues increasing in the pure shear cases (dark lines in figure 8*a*). A similar hardening behaviour, but followed by weakening at higher strain, was observed in simple shear deformation of polycrystalline aggregates of olivine [47,48] and quartz [49], and also in numerical simulations of two-phase materials [50]. Differences between laboratory experiments and numerical simulations are discussed further in §4b.

### (e) Strain localization

The observable strain localization in conjugate high-strain bands of strongly elongated grains (figure 3*a*,*e*) is quantified using a slight modification of the localization factor defined by Sornette *et al.* [51] and Gomez-Rivas [52]:
3.1
where *n*_{t} the total number of unodes and is the von Mises strain rate for every *unode n _{i}*, defined as
3.2

The strain localization factor (*f*) can range from 0 to 1, where 0 represents homogeneous deformation and 1 maximum localization.

All simulations in simple shear show an initial trend towards more homogeneous deformation, followed by a reversal when localization strengthens again (figure 9*a*). DRX suppresses localization, with the reversal occurring at a higher strain in cases with higher *N*_{DRX}. Pure shear simulations are characterized by a steady state up to a natural strain of *ε* = 0.3, followed by either a reduction of localization at high *N*_{DRX} or an increase in localization when *N*_{DRX} is low (figure 9*b*). In general, DRX suppresses localization of deformation, but only after about a natural strain of *ε* = 0.3, when a distinct LPO starts to develop.

To show the strain-rate intensity, the von Mises strain rate is calculated as a function of the symmetric strain-rate tensor (equation (3.2)) and plotted in figure 10. In simple shear simulations high strain-rate synthetic bands oriented at low angles to the horizontal shear plane (also called synthetic shear bands) develop (figure 10*a*). In pure shear simulations, high strain-rate conjugate bands form at angles of 66°–71° with respect to the vertical *y*-direction (figure 10*c*). In both simple and pure shear models, the number and intensity of shear bands are reduced by DRX (experiment 25 in figure 10*b*,*d*).

### (f) Ice simulations including polygonization

The influence of polygonization on the final microstructures was tested by analysing a series of simulations that include the ELLE polygonization subloop (i.e. experiments with suffix _N). These models were carried out in both pure and simple shear deformation and with different amount of recrystallization: 1, 10 and 20 DRX steps per deformation step, respectively (figure 11*b*,*d*,*f*). When polygonization is active the differential stress reaches a steady state earlier than in the case of the corresponding simulation without this routine (dashed lines in figure 8). Polygonization does not affect the final strong LPO regardless the amount of DRX. This applies to both simple and pure shear cases (figure 11*a*–*d* and *e*,*f*, respectively). The lack of influence of polygonization is probably due to the fact that new grains have lattice orientations close to those of their parent grains. However, the incorporation of the polygonization process, and therefore the development of proper HAGB out of tilt walls, reduces grain size compared to the cases in which this process is not included (table 2) [8]. Polygonization also changes the distribution of grain sizes, reducing the skewness and, in the case of the simple shear simulation SSH 10_N, develops a bimodal distribution (figure 12*b*,*c*).

## 4. Discussion

### (a) Effects of ice flow on the lattice preferred orientation and final microstructure

Natural deformation in glaciers and ice sheets mostly occurs under general strain conditions between the end members pure and simple shear. In general, there is an increase in vorticity with depth, with pure shear dominating at the top and simple shear at the base (e.g. [53]). It is therefore important to understand how ice behaves differently under these varying boundary conditions, especially considering the strong anisotropy (LPO) developed in ice sheets and glaciers [3,54]. Our numerical results show that the boundary conditions do not significantly change the evolution of the microstructure due to viscoplastic deformation, because elongated grains and an SPO parallel to the stretching direction and normal to the finite shortening direction develop in both pure and simple shear cases. Qualitatively, increasing DRX in simple shear simulations results in an SPO oriented at a high angle with respect to the shear plane. Simple shear deformation results in a rotation of the SPO towards the shear plane, according to the sense of shear (figures 3 and 4). GBM in isolation leads to an increase in grain size, as long as the grain boundary energy is a significant component of the driving force (equation (2.2)). It also reduces the grain elongation, hence masking strain localization. Although GBM is assumed to dominate in materials close to their melting point [55], recovery and polygonization probably play a significant role in ice sheets and glaciers. They counteract the GBM-induced grain growth, leading to a steady grain size where a balance is reached [32,56]. Adding polygonization to the model indeed leads to a reduction in grain size and grain elongation, although a steady state is not (yet) reached in our simulations (figures 11 and 12). The results are consistent with the observation that GBM and polygonization are active in large parts of ice sheets, as indicated by the absence of elongated grains in most ice core microstructures observed to date (e.g. [5,57]). Moreover, polygonization has recently been found to be also active in experimental ice creep tests close to the melting point [58].

In our simulations, a strong *c*-axes maximum initially develops parallel to the maximum compressive direction (*σ*_{1}) at a natural strain of *ε* = 0.3 in pure shear boundary conditions. When DRX is active this single maximum tends to shift towards being obliquely oriented with respect to *σ*_{1} (or *X* axis) in the *YZ* plane (figure 5). Pure shear and uniaxial compression experiments [59,60] of polycrystalline ice result in a girdle fabric of *c*-axes around the direction of *σ*_{1} (compression direction). The lack of *c*-axes oriented parallel to *σ*_{1} in our simulations is interpreted as due to the activation of recrystallization processes, coherently with temperatures near to the melting point in the experiments (between −2°C and −7°C).

In simple shear, the *c*-axes maximum initially develops parallel to *σ*_{1}, but rapidly tends to rotate with increasing strain towards being oblique or normal to the shear plane (figure 5). Similarly to the pure shear cases, activation of recrystallization processes shifts the maximum towards being oblique to the *XY* section. Ice deformation laboratory experiments revealed a similar pattern in simple shear, with the *c*-axes maximum oriented normal to the shear plane but also obliquely oriented with respect to *σ*_{1} (e.g. [61,62]).

Once an LPO forms, it does strengthen and rotate with increasing deformation according to the sense of shear in simple shear simulations (figures 4 and 5). The rotational evolution of ice fabrics in simple shear is typically characterized by the development of a strong single maximum approximately perpendicular to the shear direction [61]. With increasing strain, the obliquely oriented maximum disappears in experiments and only the maximum oriented perpendicular to the shear plane remains [61]. This observation differs from our simulations, in which the *c*-axes maximum is always obliquely oriented. This rotational evolution of ice fabrics has been recognized as vertical or near vertical single maximum LPOs in ice cores [5,63–66]. The vorticity of deformation is typically close to simple shear in the deep part of ice sheets, with the shear plane generally oriented parallel to the bedrock [57]. Another difference between laboratory experiments and numerical simulations is that the first only consider uniaxial deformation while simulations are performed under perfect two-dimensional plane strain (i.e. pure shear).

The statement that a single maximum LPO like that observed in nature under simple shear cannot be simulated unless recrystallization is taken into account [4] is not consistent with the results of our simulations, as the developing hard grains further rotate according to the sense of shear independently of the intensity of recrystallization *N*_{DRX} (figure 5, at a natural strain of *ε* = 1.2).

Our simulations show how such LPOs develop as a consequence of deformation, even when there is strong recrystallization. The implementation of the polygonization routine does not affect the developed LPO, because new grains that nucleate maintain the *c*-axis orientations close to that of their parental grains (figure 11*b*,*d*,*f*). An additional factor controlling the development of an LPO is the degree of strain localization, as illustrated by the difference in LPO developed in high and low strain-rate domains (e.g. figure 5*i*,*j*). The detailed evolution of the LPOs in different scenarios of strain localization needs to be systematically addressed in future studies.

### (b) Effects of vorticity of ice flow on rheology

All simulations present strain hardening behaviour during the first stages of deformation (figure 8), during which the degree of deformation localization remains relatively constant (figure 9*a*,*b*). At a natural strain of *ε* = 0.3 an LPO can be identified in both pure and simple shear simulations (figure 5), which result in the development of an intrinsic anisotropy that leads to an increment of strain localization in simulations performed at high strain rates (i.e. experiments SSH and PSH 0 and 1; figure 9). Differential stress reaches a steady state in simple shear simulations (figure 8*a*). This is associated with the rotation of the LPO and SPO to a more favourable position for slip, as basal planes and elongated grains progressively rotate towards the shear direction (figure 4*a*–*d*). However, the SPO and LPO develop in a less favourable orientation for easy glide in pure shear models, where they tend to rotate away from *σ*_{1} (figure 4*e*–*h*). In such cases, strain hardening keeps increasing until the end of the simulations (figure 8*a*). The development of an SPO perpendicular to the compression direction has been proposed as the cause of strain hardening in previous studies [50,67]. The formation of new grain boundaries by polygonization has a clear influence on the rheology, where the differential stress reaches a steady state earlier than in the equivalent simulation without polygonization (dashed lines figure 8*a*). This steady state is achieved when an equilibrium in grain size is reached [68], which depends on the relative amount of new grains nucleated by polygonization versus GBM [5].

A hardening-and-softening behaviour is observed in laboratory experiments under constant stress or constant strain rate in both pure and simple shear configurations [60]. In experiments, a stress peak is typically observed below a strain of 0.04. In our simulations, the transition between hardening and softening is observed at a higher strain, at approximately 0.30. There are some remarkable differences between laboratory experiments and numerical simulations. One of them is that in our models, we assume lower strain-rate conditions (10^{−10} to 10^{−12} s^{−1}) than those used in laboratory experiments, where higher strain rates are used (10^{−5} and 10^{−7} s^{−1}). This implies that our assumption of viscoplastic behaviour of ice polycrystals is perhaps not valid for the range of strain rates used in laboratory experiments, although it is more realistic for comparisons with natural cases. In these conditions, transient creep must be considered and an elastoviscous–plastic behaviour is more appropriate to describe the mechanical behaviour of the ice polycrystalline aggregate. For example, in [69], an elastoviscous–plastic model was used to simulate transient creep from the experimental data of [70]. However, this approach is limited to very low strains (less than 0.05), significantly lower than those expected in natural ice sheets. Another key difference is that a hardening-and-softening behaviour is always observed in laboratory experiments regardless the stress configuration. Contrarily, this is only observed in our models in simple shear conditions. We have to take into account that the strain rate and geometry of deformation are constant during the simulation time. Only strain hardening occurs in numerical pure shear models, at least up to the end of the simulation time (*ε* = 1.2). Our results indicate that tertiary creep is probably not steady state, and that hardening and softening can probably occur at large strain. This would be associated with an increase of the strength of the LPO and the viscous anisotropy of the material.

Simulations performed at low strain rates do not show high-strain localization bands (figure 10*b*–*d*), as deformation is more homogeneously distributed (figure 9, experiments SSH and PSH 10 and 25). GBM and polygonization result in microstructures dominated by larger and more equidimensional grains, preventing the formation of strong grain elongation and masking localization bands. In that situation and with a strong LPO, deformation cannot be accommodated by high-strain bands any more, leading to the activation of the non-basal slip systems in both pure and simple shear simulations (figure 6*b*). The favourable orientation of the pyramidal plane with respect to *σ*_{1} in pure shear models (figure 5*e*–*h*) enables dislocations to glide along this system of planes rather than along the basal ones (figure 6*b*, dashed lines). However, the pyramidal plane does not reach this favourable orientation with respect to *σ*_{1} in simple shear simulations, due to the rotation of the LPO with the sense of shear (figure 5*a*–*d*). This results in the basal slip system remaining more active than the non-basal ones throughout the simulations (figure 6*b*, solid lines). This evolution of changing slip system activity has to be further investigated, as its occurrence has been proposed to be a rate-limiting mechanism not only at high stresses [71], but also applicable to low stress conditions.

The activation of non-basal slip systems can be recognized from the rotation of the *a*-axes towards the maximum elongation axis, while *c*-axes rotate towards the direction of the maximum finite shortening (figure 4*c*,*d*,*g*,*h*). Recently, clustering of *a*-axes in a wide girdle at a high angle to the compression direction was observed in experiments with DRX analysed by electron backscatter diffraction (EBSD) [58]. Clustering of *a*-axes has been previously described in the deep part of the GRIP ice core, where the slip along non-basal systems is proposed as a possible cause of the *a*-axis concentration [72]. The activation of the non-basal slip systems is related to the development of subgrain boundaries, as shown in figure 7*b*,*d*. Subgrain boundaries composed of basal (tilt and twist boundaries) and non-basal dislocations (tilt boundaries) were described using EBSD datasets, where the non-basal dislocations play an important role in the formation of all subgrain boundaries [73,74]. The frequent abundance of different subgrain boundary types (figure 13), especially the ‘p-type’ subgrain boundaries (more than 50% of all subgrain boundaries in EDML [75]), which cannot be basal twist boundaries only but also non-basal tilt boundaries (more than 75% according to very first estimates, table 2 [74]), gives indications of the importance of non-basal dislocation activity. This has to be treated with care, as only recovered dislocations can be measured, forming subgrain boundaries, and not all dislocations involved in deformation. Clearly, more high-resolution full-crystal orientation measurements are needed to validate our models with respect to magnitude of recovery versus the activation of different slip systems during deformation.

## 5. Conclusion

We present a state of the art of a full-field numerical approach applied to ice. It includes simulations with and without recrystallization under pure and simple shear boundary conditions. Although the amount of recrystallization strongly affects the grain morphology, it has a much smaller effect on LPO development. In all cases, regardless of the amount of DRX and ice flow, a single *c*-axes maximum develops that is oriented approximately perpendicular to the maximum finite shortening direction and which in simple shear rotates towards the normal to the shear plane. This leads to distinctly different behaviour in pure and simple shear. In pure shear, the LPO and SPO are increasingly unfavourable for deformation, leading to hardening and to an increased activity of non-basal slip. The opposite happens in simple shear, where vorticity causes rotation of the LPO and SPO to a favourable orientation resulting in strain softening. The increase of recrystallization enhances the activity of the non-basal slip, due to the reduction of deformation localization. In pure shear conditions the pyramidal slip activity is thus even more enhanced and can become higher than the basal-slip activity.

Our results further show that subgrain boundaries can be developed by the activity of the non-basal slip systems. The implementation of the polygonization routine reduces grain size and SPO, but does not significantly change the final LPO, because newly nucleated grains approximately keep the *c*-axis orientations of their parental grains. However, it enables the establishment of a grain size equilibrium, and therefore the differential stress reaches a steady state.

## Authors' contributions

M.-G.L., P.D.B., A.G. and I.W. conceived the study. M.-G.L. wrote the manuscript together with A.G., P.D.B., E.G.R. and I.W. D.J. and R.A.L. critically discussed the content and edited the manuscript. A.G., R.A.L., M.-G.L., P.D.B., J.R., E.G.R. and F.S. contributed to the development of the numerical codes. All authors gave final approval for publication.

## Competing interests

The authors declare that there are no competing interests.

## Funding

M.-G.L. was funded by the programme on Recruitment of Excellent Researchers of the Eberhard Karls Universität Tübingen. F.S. was funded by DFG (SPP 1158) grant BO 1776/12-1. The Microdynamics of Ice (MicroDICE) research network, funded by the European Science Foundation, is acknowledged for funding research visits of M.-G.L.

## Acknowledgements

We thank all the members of the ELLE development group, in particular, Lynn Evans, Sandra Piazolo and Verity Borthwick for their contributions to the simulation code. We gratefully acknowledge A. Treverrow and an anonymous reviewer, whose constructive reviews greatly improved the manuscript, together with the editorial guidance of P. Sammonds. This work was carried out as part of the Helmholtz Junior Research group ‘The effect of deformation mechanisms for ice sheet dynamics’ (VH-NG-802).

## Footnotes

One contribution of 11 to a theme issue ‘Microdynamics of ice’.

- Accepted November 10, 2016.

- © 2016 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.