Skip to main content Accessibility help


  • Access
  • Open access
  • Cited by 13



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Dynamic recrystallisation of ice aggregates during co-axial viscoplastic deformation: a numerical approach
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Dynamic recrystallisation of ice aggregates during co-axial viscoplastic deformation: a numerical approach
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Dynamic recrystallisation of ice aggregates during co-axial viscoplastic deformation: a numerical approach
        Available formats
Export citation


Results of numerical simulations of co-axial deformation of pure ice up to high-strain, combining full-field modelling with recrystallisation are presented. Grain size and lattice preferred orientation analysis and comparisons between simulations at different strain-rates show how recrystallisation has a major effect on the microstructure, developing larger and equi-dimensional grains, but a relatively minor effect on the development of a preferred orientation of c-axes. Although c-axis distributions do not vary much, recrystallisation appears to have a distinct effect on the relative activities of slip systems, activating the pyramidal slip system and affecting the distribution of a-axes. The simulations reveal that the survival probability of individual grains is strongly related to the initial grain size, but only weakly dependent on hard or soft orientations with respect to the flow field. Dynamic recrystallisation reduces initial hardening, which is followed by a steady state characteristic of pure-shear deformation.


Ice is a common mineral on the Earth's surface. Thick ice sheets cover Greenland and Antarctica. How much ice is stored in these ice sheets depends on its mechanical properties, as the weight of the ice masses causes a steady flow of ice towards the surrounding oceans. The mechanical properties of ice, therefore, directly impact on human society through its role in controlling sea level (Kennicutt and others, 2014). At natural strain-rates of <10−6 s−1 (e.g. Budd and Jacka, 1989), ice is a ductile material, assumed to deform dominantly by dislocation creep (e.g. Schulson and Duval, 2009), similar to other rock-forming minerals, such as quartz and olivine at much higher absolute, but similar homologues temperatures. In contrast to these minerals, ice on Earth only occurs very close to its melting temperature.

Ice on Earth is exclusively Ih, with a hexagonal structure. The strong viscoplastic anisotropy of single ice crystals constitutes the first complexity in describing and modelling the behaviour of the flow of ice (Duval and Castelnau, 1995). The critical resolved shear stress (CRSS) to activate dislocation slip on the basal plane is ~60–100 times lower than on prismatic and pyramidal planes (Duval and others, 1983; Wilson and Zhang, 1996; Montagnat and others, 2013). Strain is therefore mostly achieved by slip on the basal plane, with relatively minor contributions of the other slip planes (Ashby and Duval, 1985).

The bulk behaviour of large ice masses results from the behaviour of the ensemble of individual ice grains. This behaviour is thus strongly influenced by the viscoplastic anisotropy of these grains and the orientation of their lattices, usually expressed in the orientation of their c-axes. During flow, ice quickly develops a lattice preferred orientation (LPO) by lattice rotation (Azuma and Higashi, 1985; Alley, 1988, 1992; Jacka and Li, 2000). The ice microstructure further evolves through time by recrystallisation, that is any re-orientation of the lattice caused by grain boundary migration (GBM) and/or formation of new grain boundaries (GBs; see Glossary in Faria and others, 2014a, b, c), which is driven by GB surface energy and strain energy. The mechanical properties of ice masses are therefore not constant, but are the result of complex, time-dependent coupling between deformation and microstructural evolution. Understanding this coupling is of paramount importance for studies on the flow behaviour of polar ice sheets.

1.1. Flow of polar ice sheets

Experiments on ice show it has a strongly strain-dependent behaviour (Jacka, 1984; Azuma, 1994; Treverrow and others, 2012; Budd and others, 2013). Initial, transient or primary creep is dominated by dislocation glide on the basal plane and is associated with an increase in strength (Montagnat and others, 2006). Secondary creep is reached at a few percentage strain, when recrystallisation sets in, the maximum strength is reached and thereby the material begins to soften again (Glen, 1955). An LPO develops during tertiary creep where a steady state may be reached after ~10% strain (Budd and Jacka, 1989; Treverrow and others, 2012). Although secondary creep is highly transient, it is the stage on which the commonly used flow law is based. ‘Glen's law’ (Glen, 1958) describes the relation between deviatoric stress (σ) and strain-rate ( $\dot \varepsilon $ ) by a power law with a stress exponent, n:

(1) $$\dot \varepsilon = {A_{\left( {T{\rm,} \mu s} \right)}}{\sigma ^n}$$

This law is based on the maximum strength that is reached during secondary creep. In most studies, it is assumed that n = 3, although values ranging from 1 to 5 have also been reported, for example, based on experiments (Goldsby and Kohlstedt, 2001) and observations on ice sheets (Cuffey and Kavanaugh, 2011; Gillet-Chaulet and others, 2011). The pre-exponential factor depends on temperature (T), but also on the microstructure (μs), in particular the LPO. The effect of LPO is usually incorporated by using an LPO-dependent scalar ‘enhancement factor’ (see Faria and others, 2014b). This enhancement factor changes the absolute strain-rate at a given stress, but not the sensitivity of strain-rate to a change in stress, defined by n.

1.2. Microstructures of deforming ice

Laboratory experiments and deep drill cores through glaciers and ice sheets are the most important sources of information on ice microstructures (e.g. Weikusat and others, 2009; Piazolo and others, 2012; Wilson and others, 2013; Faria and others, 2014a, b). The three dominant processes that affect the shape and size of grains are sub-grain formation or polygonisation, GBM driven by strain energy and grain growth by surface energy (Urai and others, 1986). All three processes operate concurrently, albeit in different proportions, mainly depending on strain-rate and temperature (Faria and others, 2014b). Grain boundary surface energy-driven grain growth in isolation (‘normal grain growth’; NGG) leads to an increase in grain size by reduction of GB surface area (Azuma and others, 2012). Polygonisation, in geology also termed ‘rotational recrystallisation’ (RRX) (Urai and others, 1986), is a recovery process whereby the lattices in regions within grains rotate relative to each other by the formation of tilt walls. Progressive rotation causes tilt walls to become true GBs, which causes a reduction of the grain size. A dynamic balance between grain growth and polygonisation is invoked by some authors (Mathiesen and others, 2004; Durand and others, 2008; Roessiger and others, 2011, 2012; Azuma and others, 2012). The third process, also called ‘strain-induced boundary migration’ (SIBM) (Urai and others, 1986; Drury and Urai, 1990; Tullis and others, 1990; Stipp and Tullis, 2003) is driven by the energy stored within the lattice in the presence of dislocations. Variations in stored energy (or strain energy) cause GBs to bulge into regions with highest stored energy, leading to the formation of irregular or lobate GBs.

Based on grain-size profiles, it is proposed that the three recrystallisation processes dominate at different depths (De la Chapelle and others, 1998): NGG in the upper few hundreds of meters, followed by RRX and finally SIBM in the lowest few hundreds of meters in ice sheets, where a sudden increase in grain size is observed in drill cores (e.g. Faria and others, 2014a). However, several authors (Kipfstuhl and others, 2006, 2009; Weikusat and others, 2009; Faria and others, 2014c) have challenged the tri-partite model for ice sheets and showed that all three processes operate at all levels, albeit in different proportions, decoupling the recrystallisation phenomenon from simple depth considerations and describing grain size evolution by strain-rate and temperature conditions (Faria and others, 2014c).

Dislocation creep is assumed to be largely insensitive to grain size and shape (Poirier, 1985; Schulson and Duval, 2009). However, recrystallisation is expected to strongly influence flow properties of ice, if it influences the LPO (Duval and Castelnau, 1995; Duval and others, 2000; Kocks and others, 2000). For example, if crystals with certain orientations are more prone than others to be consumed by migrational recrystallisation owing to their higher dislocation density, migrational recrystallisation would directly influence the LPO, and hence the strength. The study of the interaction between deformation and recrystallisation is thus clearly relevant to improve our understanding and quantification of the flow of ice, from a small volume to that of polar ice sheets.

1.3. Numerical modelling of microstructural evolution in ice

Much insight in the mechanical and microstructural behaviour of ice has been gained from experiments (Azuma and Higashi, 1985; Wilson, 1986; Budd and Jacka, 1989; Jacka and Li, 2000; Treverrow and others, 2012; Piazolo and others, 2012). However, one major limitation of experiments is that deformation rates are inevitably higher than those typical in nature by several orders of magnitude, especially in case of slowly flowing polar ice sheets (ε ≈ 1010−10−13 s−1; Budd and Jacka, 1989; Montagnat and Duval, 2004). Numerical modelling provides a solution to this, as it is limited by computational power and the sophistication of the algorithms, but not by scale, stress or strain-rate.

Theoretical and computational studies have been performed for many years in order to understand the deformation and microstructural evolution (Sachs, 1928; Taylor, 1938). Crystal plasticity based finite element (FE) implementations have been used in the past to model plastic deformation of polycrystalline materials (Becker, 1991; Mika and Dawson, 1998; Raabe and others, 2001; Delannay and others, 2003; Diard and others, 2005). However, the large number of degrees of freedom needed in FE calculations, which restricts the size of the simulated microstructure and the use of averaged properties, are important limitations. To overcome these and to be able to predict the micromechanical behaviour, a full-field formulation based on fast Fourier transforms (FFT) (Moulinec and Suquet, 1994) is used to simulate plastic deformation in heterogeneous materials, such as composites (Moulinec and Suquet, 1998; Michel and others, 2000; Idiart and others, 2006) or polycrystals (Lebensohn, 2001; Lebensohn and others, 2004, 2005, 2008). Polycrystal full-field models differ from mean-field models (Sachs, 1928; Taylor, 1938; Lebensohn and Tomé, 1993) in that the latter explicitly or implicitly assume homogeneous deformation at single crystal level.

Previous full-field simulations applied to ice were mostly limited to low-finite strain (<4% of shortening; Montagnat and others, 2011) with the exception of a brief example of high-strain FFT modelling by Montagnat and others (2013). Most studies dealing with FFT-modelling of ice only considered viscoplastic deformation and only a few include recrystallisation (Durand, 2004; Montagnat and others, 2013). Considering that ice is near to its melting point at natural conditions, a numerical study on ice recrystallisation processes coetaneous with deformation up to high-strain, the subject of this paper, is long overdue.

In this contribution, we present models combining FFT-based simulations to determine local plastic deformation with recrystallisation up to high-strain. The model will be illustrated with pure-shear, plane-strain deformation of pure, polycrystalline Ih ice. We will show how recrystallisation has a major effect on the microstructure, but, perhaps counter-intuitively, only a minor effect on LPO development. Although LPO's do not vary much, recrystallisation appears to have a distinct effect on the relative activities of slip systems. A limitation of the present approach is that only GBM and sub grain formation were considered in the simulations. However, formation of new high-angle GBs due to subgrain rotation or nucleation in highly deformed regions is not yet implemented.


In this study, we numerically model the microstructural evolution of an aggregate of pure ice grains during deformation and recrystallisation in 2-D pure shear. The approach is based on the coupling between modelling deformation of the viscoplastic polycrystal with a full-field theory formulation (numerically implemented in the FFT code by Lebensohn (2001) and Lebensohn and others (2008)) and the microstructural modelling platform, ELLE (Jessell and others, 2001a; Bons and others, 2008) designed to simulate recrystallisation processes. The FFT approach provides the stress and velocity fields, from which deformation induced lattice rotation and an estimate of geometrically necessary dislocation densities are calculated. These data, as well as a map of the grain boundary network are used to simulate intracrystalline recovery and GBM driven by the reduction of surface energy and stored strain energy generated by dislocations.

The data structure of the models consists of two basic layers (Fig. 1a and b): (1) a contiguous set of polygons (termed flynns) that are themselves defined by boundary nodes (bnodes) and connected by straight boundary segments and (2) a set of unconnected nodes (unodes) that provide a higher resolution of physical properties within grains. The correspondence between both codes is done using the unodes layer that allows a direct mapping between the regular computational grid required by the FFT code and the material information required by ELLE.

Fig. 1. Basic data structure and set up of the simulations. The ELLE data structure has two different layers: (a) boundary nodes (double or triple nodes) that define polygons (grains), (b) a regular mesh of unconnected nodes to store dislocation density and lattice orientations used for the FFT calculation. Initially 2 × 1 model with 3260 grains (c) is deformed in pure shear with 2% vertical shortening every step (d).

In the next sections, we briefly provide the basics of the governing and constitutive equations and numerical implementation of both approaches.

2.1. Viscoplastic deformation – FFT

The FFT-based formulation provides an exact solution of the mechanical problem by finding a strain-rate field and a stress field, associated with a kinematically admissible velocity field, which minimises the average local work rate under the compatibility and equilibrium constraints (see Lebensohn, 2001; Lebensohn and others, 2008 for a detailed description of the method, and Montagnat and others, 2013 for its application to ice). The method is based on the fact that the local mechanical response of a periodic heterogeneous medium can be calculated as a convolution integral between the Green function of a linear reference homogeneous medium and the actual heterogeneity field. The use of the FFT algorithm reduces the convolution integrals in Cartesian space to simple products in Fourier space and allows us to obtain the mechanical fields by transforming that product back to real space. As the heterogeneity field depends on the unknown mechanical field, an augmented Lagrangian-based iterative scheme is used to determine stress and strain-rate fields that fulfil the compatibility and equilibrium constraints (Michel and others, 2001). In general, the FFT approach has a better numerical performance than simulations based on the FE Method (Prakash and Lebensohn, 2009; Roters and others, 2011), but it is restricted to periodic boundary conditions and requires discretisation into a regular grid.

The constitutive behaviour of a polycrystal aggregate is defined using a nonlinear viscous rate-dependent crystal plasticity approach where deformation on the crystal scale is assumed to be accommodated by dislocation glide only. The constitutive equation between the strain-rate and the deviatoric stress σ'(x) at position x is defined using a nonlinear viscous response given by

(2) $$\eqalign{{{\dot \varepsilon} _{ij}}\left( x \right) & = \mathop \sum \limits_{s{\rm = 1}}^{{N_{\rm s}}} m_{ij}^{\rm s} \left( x \right){{\dot \gamma} ^{\rm s}}\left( x \right) \cr & = {{\dot \gamma} _0}\mathop \sum \limits_{s = 1}^{{N_{\rm s}}} m_{ij}^{\rm s} \left( x \right){\left \vert {\displaystyle{{{m^{\rm s}}\left( x \right):\sigma ^{\prime}\left( x \right)} \over {{\tau ^{\rm s}}\left( x \right)}}} \right \vert ^n}sgn\left( {{m^{\rm s}}\left( x \right):\sigma ^{\prime}\left( x \right)} \right)} $$

where the sum runs over all N s slip systems, m s is the symmetric Schmid tensor defined by the dyadic product of the vector normal to the glide plane and the Burgers vector of slip system s, and ${\dot \gamma ^{\rm s}}$ and τ s are, respectively, the shear strain-rate and the CRSS of slip system s, ${\dot \gamma _0}$ is a reference strain-rate and n the stress exponent (inverse of the strain-rate sensitivity). The CRSS associated with each material point can potentially change according to a hardening or softening law (e.g. Lebensohn and others, 2008). This was not considered in this study, where all CRSSs  are kept constant during simulations.

After each time increment Δt, the microstructure is updated using an explicit scheme. The new position of a point X of the Fourier grid is determined using the velocity fluctuation term ${\tilde v_{\rm i}}\left( {\bf x} \right)$ (Eqn (25) in Lebensohn, 2001) arising from the heterogeneity field as

(3) $${X_i}\left( {\bf x} \right) = x_i^0 + \left( {{{\dot E}_{ij}}x_j^0 + {{\tilde v}_i}\left( {\bf x} \right)} \right) \times \Delta t,$$

where ${\dot E_{ij}}$ is the macroscopic strain-rate. The local crystallographic orientations are updated according to the local lattice rotation,

(4) $${\omega _{ij}}\left( {\bf x} \right) = \left( {{\dot \Omega} \left( {\bf x} \right) + \widetilde{{\dot \omega}} \left( {\bf x} \right) - \dot \omega _{ij}^{\rm p} \left( {\bf x} \right)} \right) \times \Delta t,$$

where ${\dot \Omega} \left( {\bf x} \right)$ is the average rotation rate of the macroscopic velocity gradient tensor, $\widetilde{{\dot \omega}} \left( {\bf x} \right)$ is the local fluctuation in rotation-rate obtained by taking the anti-symmetric part of the fluctuation term ${\tilde v_i}\left( {\bf x} \right)$ and $\dot \omega _{ij}^{\rm p} \left( {\bf x} \right)$ is the plastic rotation-rate of the crystal lattice, calculated as $\sum\nolimits_{s{\rm = 1}}^{{N_{\rm s}}} {\alpha _{ij}^{\rm s} \left( {\bf x} \right){{\dot \gamma} ^{\rm s}}\left( {\bf x} \right)} $ , where $\alpha _{ij}^{\rm s} $ is the anti-symmetric Schmid tensor.

An estimate of the dislocation-density field can be found using the strain gradient plasticity theory (e.g., Fleck and Hutchinson, 1997; Gudmundson, 2004; Wen and others, 2005), based on Ashby's concept of geometrically necessary dislocations (GND) (Ashby, 1970). Following this approach, GND are required to maintain strain compatibility across the material. The presence of GND (ρ) produces a strain energy storage associated with lattice distortion that, according to Ashby (1970), is related to the absolute value of the effective plastic strain gradient (e.g. Brinckmann and others, 2006; Estrin and Kim, 2007)

(5) $$\rho = \displaystyle{{{\eta ^{\rm p}}} \over b},$$

where b is the magnitude of the Burgers vector and η p is the effective plastic strain gradient derived by Gao and others, (1999) as

(6) $${\eta ^{\rm p}} = \sqrt {\displaystyle{1 \over 4}\eta _{ijk}^{\rm p} \eta _{ijk}^{\rm p}}, $$


(7) $$\eta _{ijk}^{\rm p} = {\varepsilon _{ik{\rm,} j}} + {\varepsilon _{\,jk{\rm,} i}} - {\varepsilon _{ij{\rm,} k{\rm \;}}}, $$

where ε ik,j indicates the partial derivatives of the plastic strain tensor with respect to the reference coordinates. Although other, higher order approaches to calculate the dislocation densities resulting from the activity of each slip systems are possible (see e.g. Ma and others, 2006; Roters and others, 2010), or the Nye's dislocation density tensor can be used (e.g., Nye, 1953; Kröner, 1962; Ashby, 1970), we simplify our approach by assuming a constant Burgers vector for all systems, and therefore, dislocations arising from basal and nonbasal slip systems are not differentiated. Furthermore, statistically stored dislocations, dislocation tangles, etc. are not considered as their type, activity and density cannot be determined from the model calculations, but would require additional assumptions. Dynamic dislocation effects, such as avalanches (Weiss and Marsan, 2003) are also not considered for the same reason.

2.2. Recrystallisation – ELLE

The numerical platform, ELLE was used in this study to simulate the recrystallisation processes. The modular programming of ELLE provides a highly versatile data transfer and allows coupling between various processes that affect a microstructure (Jessell and others, 2001a, b; Bons and others, 2008; Piazolo and others, 2010). ELLE has been used to simulate a variety of microstructural processes, including dynamic recrystallisation (DRX) (Piazolo and others, 2002, Montagnat and others, 2013), grain growth (Jessell and others, 2001b; Jessell and others, 2003; Roessiger and others, 2011), strain localisation (Jessell and others, 2005; Griera and others, 2011), porphyroclast rotation (Griera and others, 2013), deformation of two phase materials (Jessell and others, 2009) and folding (Llorens and others, 2013a, b).

In this study, GBM and recovery processes are used in order to simulate the microstructural evolution by recrystallisation. GBM covers the motion or displacement of high-angle boundaries (HAGB), whereby lattice distortions in swept regions are restored. Recovery is the decrease in intracrystalline heterogeneities by means of local rotation of the lattice without motion of HAGBs.

2.3. GBM

The model aims to simulate GBM that is driven by reduction of the GB energy and stored strain energy resulting from lattice distortions (i.e. dislocation density). For a single-phase, polycrystalline aggregate, the velocity of a grain boundary is expressed as

(8) $$\vec v = M\Delta f\vec n,$$

where M corresponds to the boundary mobility, Δf to the driving pressure (force divided by  GB area it acts on) and $\vec n$ to the unit vector normal to the boundary in the direction of movement. The GB mobility is highly dependent on temperature and can be expressed using an Arrhenius or exponential equation

(9) $$M = {M_0}ex{\,p^{( - Q/RT)}},$$

where M 0 is the intrinsic mobility, Q is the boundary-diffusion activation energy, R is the universal gas constant (e.g. Gottstein and Mecking, 1985). The driving pressure Δf is defined as

(10) $$\Delta f = \Delta H - \displaystyle{{2{\gamma _{\rm s}}} \over r},$$

where ΔH is the difference in stored strain energy density across the grain boundary, γ s to the boundary energy and r to the local radius of curvature of the grain boundary. The stored energy ΔH is calculated using the difference in dislocation density Δρ across a grain boundary and the energy of a dislocation per unit length E ρ (Cottrell, 1953; Karato, 2008):

(11) $$\Delta H = \mathop \sum \limits_{s = 1}^{{N_{\rm s}}} {\left( {\Delta \rho {E_\rho}} \right)_{\rm s}} = \mathop \sum \limits_{s = 1}^{{N_{\rm s}}} {\left( {\rho \displaystyle{{G{b^2}} \over {4\pi P}}\log \displaystyle{{{R_{\rm e}}} \over {{b_0}}}} \right)_{\rm s}}\sim \Delta \rho G{b^2},$$

where G is the shear modulus, P is a constant that depends on dislocation type and Poisson's ratio, b 0 is the radius of the dislocation core and R e is the upper limit for integration of the distance from the dislocation line. As the dislocation density for each slip system is unknown, we make the assumption that the elastic energy is the same for all slip systems and Eqn (11) can be approximated to $\Delta H \sim \Delta \rho G{b^2}$ (e.g. Cottrell, 1953; Karato, 2008).

GBM is simulated using a front-tracking approach based on the algorithm by Becker and others (2008) and recently used by Roessiger and others (2012) to simulate grain growth of ice with air bubbles. In this algorithm, the motion of boundaries is driven by minimisation of the Gibbs free energy, according to Eqn 10. For each time increment, each bnode is moved over a small distance assuming that the velocity does not change during the small time increment. In order to predict the direction and magnitude of the bnode displacement, four orthogonal trial positions at very small distance from the current node position are chosen. The change of total energy (E j ) for each trial position is due to (1) the change in total boundary length, and therefore the change of total boundary energy and (2) the decrease of dislocation density due to the removal of all dislocations in the swept area,

(12) $${E_j} = {A_j}\overline {{\rho _j}} {E_\rho} - \gamma \mathop \sum \limits_{i{\rm = 1}}^N {l_i}\;, $$

where A j , $\; \overline {{\rho _j}} $ and l i correspond to the swept area, to the average dislocation density of the swept area and the length of the boundary segments due to the j-th trial position, respectively. N is the number of segments (two or three) that connect at the node under consideration. From the four E j values, the spatial gradient in free energy can be calculated using a central finite difference approach. The direction of bnode displacement is thus determined by the direction of maximum energy reduction. Then, the velocity of a boundary is calculated using Eqn (8) assuming that the rate of energy dissipation is equal to the sum of the rates of work done by each individual boundary segment linked to the moving node. A more detailed description of this node movement algorithm can be found in Becker and others (2008) and Bons and others (2008).

During each time step, bnodes are sequentially moved in randomised order. The movement is calculated and applied immediately for each individual bnode. After each displacement, the topology of the grain network is examined and modified as needed. When a boundary sweeps across a unode, the local information stored in this material point is updated with the following rules: (1) its dislocation density value is set to zero and (2) its lattice orientation is reassigned to that of the nearest unode in the grain it now belongs to (Fig. 2).

Fig. 2. (a) Example of the simulation of GBM using two neighbouring grains with different dislocation densities, stored in the unconnected nodes (squares). (b) The difference in dislocation density drives the grain boundary into the grain with the highest dislocation density. The dislocation density in swept nodes is set to zero. (c) Surface energy strives to reduce grain boundary length, which is achieved by moving the boundary in the direction r towards the centre of curvature. Actual movement in one calculation step are much smaller, <1% of the distance between boundary nodes.

2.4. Recovery

In this study, the term recovery refers to intragranular reduction of the stored energy in a deformed crystal by annihilation of distributed dislocations and their rearrangement into lower-energy configurations in the form of regular arrays or low-angle subgrain boundaries, i.e. polygonisation (Sellars, 1978; Urai and others, 1986; McQueen and Evangelista, 1988; Borthwick and others, 2013; Faria and others, 2014b). To simulate this numerically, the approach by Borthwick and others (2013) is implemented in the ELLE modelling platform. The numerical approach was validated using the results from intracrystalline evolution during annealing experiments of deformed salt single crystals (Borthwick and Piazolo, 2010; Borthwick and others, 2013). Briefly, the model assumes that the rotation rate of a crystallite/subgrain is proportional to the torque (Q') generated by the change of surface energy associated with the reduction in misorientation (e.g. Li, 1962; Erb and Gleiter, 1979; Randle, 1993). Analogous to GBs, where the velocity is proportional to a driving force, a linear relationship between the angular velocity ω and a driving torque Q' is assumed:

(13) $$\omega = M^{\prime}Q^{\prime},$$

where M' is the rotational mobility (Moldovan and others, 2002). Considering a 2-D microstructure, the torque acting on a subgrain delimited by sb subgrain boundaries is given by

(14) $$Q^{\prime} = \mathop \sum \limits_{{\rm sb}} {l_{{\rm sb}}}d{\gamma _{{\rm sb}}}/d{\theta _{{\rm sb}}},$$

where l sb denotes the boundary length with grain boundary energy γ sb and misorientation angle θ sb across the boundary between the reference subgrain and a neighbouring subgrain sb. For low-angle boundaries, the boundary energy as a function of misorientation angle is calculated with the Read–Shockley equation (Read and Shockley, 1950). For this misorientation angle, the equivalent rotations by crystal symmetry are considered in order to obtain the minimum misorientation θ between two crystallite/subgrains.

Note that Eqn (14) assumes that the torque is calculated with respect to an axis through the centre of mass of the subgrain. If grain rotation is assumed as a viscous process and it is accommodated by either cooperative motion/rearrangement of dislocations in the boundary (Li, 1962) and/or by boundary and lattice diffusion (Moldovan and others, 2001, 2002), M' can be expressed as:

(15) $$M^{\prime} = {D_{{\rm gb}}}\displaystyle{{\delta {\rm \Omega}} \over {KT}}\displaystyle{1 \over {{d^{\rm p}}}},$$

where D gb is the grain boundary self-diffusion, δ the diffusion width of (sub-)grain boundary, Ω is the atomic volume, K the Boltzmann's constant, T is the temperature and d is the (sub-)grain-size. The exponent p varies according to the assumed accommodation process. Results from atomistic simulations and phase-field models by Upmanyu and others (2002) indicated a value of exponent p = 3. The assumption of crystallite geometry as columnar with unit column height includes a geometric factor that raises p to 5. For a temperature of 243 K, the rotational mobility M' is 5.5 × 10−7 Pa−1 s−1 m−2 (list of values in Table 1).

Table 1. Explanation of symbols used in text. Values used in the current models are given in parentheses

Recovery is implemented in ELLE using the unconnected nodes layer and assuming that each unode is a potential subgrain. As previously indicated, the goal of this routine is to reduce the intracrystalline heterogeneities by means of local rotation of the crystal lattice at each unode resulting in a decrease of the local misorientation, development of areas of homogenous lattice orientations and the evolution of subgrain boundaries (Borthwick and others, 2013). Briefly, the process follows the same philosophy as for GBM. However, instead of using trial positions, trial rotations of the lattice are now used. The algorithm starts by choosing a random unode and finding the neighbouring unodes that belong to the same grain. Then, the average misorientation $\widetilde{{{\theta _{\rm i}}}}$ and boundary energy $\widetilde{{{\gamma _{\rm i}}}}$ between the reference unode i and the neighbour unodes j are calculated.

(16) $$\widetilde{{{\theta _i}}} = \mathop \sum \limits_{\,j{\rm = 1}}^N {\theta _{ij}},\; \widetilde{{{\gamma _i}}} = \mathop \sum \limits_{\,j{\rm = 1}}^N \displaystyle{{\gamma ({\theta _{ij}})} \over N},$$

where θ ij is the minimum misorientation angle after taking the symmetry operators into account.

A number of potential rotation axes are initially selected according to the potential specific active slip-system during deformation. During the routine, the predefined crystallographic axes are selected sequentially and a very small angle of rotation is carried out. Then, the new $\widetilde{{{\gamma _i}}}$ is calculated and compared with the initial energy: the crystal orientation of the reference unode is rotated towards the value that results in the maximum reduction in energy. If all trials result in an increase in the energy of the neighbourhood, the crystal orientation at the datum point is left unchanged. A difference to the original approach by Borthwick and others (2013), where a constant rotation rate is applied independently of the torque component, the rotation rate is calculated here by using the torque as a driving force according to Eqn (13). To simplify the routine, the boundary length of all datapoints is considered to be equal. This procedure is repeated for each unode in a random order.

Finally, dislocation density in the model updated as a change of misorientation angle implies a modification of them following a proportional relationship between both variables.

2.5. Program flow and coupling the viscoplastic FFT code with ELLE

A resolution of 256 × 256 Fourier points (unodes) is used to map lattice orientations, resulting in a unit cell defined by 65,536 discrete nodes. Each unode represents a small area or crystallite with a certain lattice orientation and dislocation density. The lattice orientation is defined by three Euler angles. The FFT module uses these unodes, not the GB network. The data structure of ELLE is fully periodic wrapping: grains touching one side of the model continue on the other size. This feature not only reduces boundary effects, but is a requirement for the FFT code. The shape of grains changes by changing the position of bnodes, either according to the velocity field (deformation) or by GBM. Distances between nodes are kept between 5.5 × 10−3 and 2.5 × 10−3, the unit distance, by removing double nodes when their neighbours are too close or adding double nodes when two nodes are too far apart. When two triple nodes approach each other within a distance of 2.5 × 10−3, a neighbour switch is carried out. Small two- or three-sided grains are removed when they are defined by triple nodes only and their distance becomes <2.5 × 10−3, the unit distance.

In ELLE, each process is simulated with a separate module that acts on the data structure (Jessell and others, 2001a). Each process is activated sequentially in a loop that represents a small time increment (∆t). The optimal ∆t depends on the process. To avoid using a small ∆t, required by one process, for a computationally expensive process that allows a larger ∆t, processes may be called on a different number of times within each loop. An experimental run (Fig. 3) consists of iterative applications of small increments of 2% of shortening in vertical coaxial compression by viscoplastic deformation (FFT), followed by GBM and recovery. To avoid large steps during simulation of recrystallisation, GBM and recovery are simulated using a sub-loop where both processes are combined using small time increments (∆t/5). This methodology avoids the influence of order in running the recrystallisation processes and increases the stability of the numerical solution. Resolution tests were performed to ensure that the temporal and spatial resolutions are adequate.

Fig. 3. Schematic program flow. The initial microstructure is subjected to a closed loop of two alternating processes or modules: viscoplastic deformation (FFT) and DRX. The DRX step itself consists of a number of alternating GBM, and recovery calculations.

Each simulation of coupled deformation and recrystallisation begins with a viscoplastic deformation step. After numerical convergence of the viscoplastic FFT, the new positions of Fourier points are updated using Eqn (17). As Fourier points and unodes in ELLE are equivalent, the position and material properties (i.e. Euler angles, dislocation density) of unodes are directly updated. During transfer information between programs, we assumed that the micromechanical fields are constant during the incremental step. Based on the velocity field, the new locations of bnodes are calculated using a linear inverse weighted interpolation,

(17) $$x_i^{{\rm t + \Delta t}} = x_i^{\rm t} + \Delta t\mathop \sum \limits_{i{\rm = 1}}^{\rm n} \displaystyle{{\mathop \sum \nolimits_{i{\rm = 1}}^{\rm n} d\left( {{x_i} - {x_j}} \right)} \over {d\left( {{x_i} - {x_j}} \right)}}{v_j},$$

where only the velocity v j from unodes located at a distance d(x i x j ) below a threshold is taken into account. After an increment of deformation, the unodes no longer define a regular grid and a new regular mesh of Fourier points is required in order to run the next step. A particle-in-cell approach is used to remap all properties from the current configuration to a new regular Fourier grid. In order to avoid unrealistic lattice orientations, Euler angles of new Fourier nodes are not interpolated and they are defined using the lattice orientation of the nearest unode that belongs to the same grain. For the case of viscoplastic materials, as simulated here, the constitutive relations are defined by means of quantities defined in the current configuration only and the approach allows us to run numerical simulations up to high-strains. This approach was verified by Griera and others (2011, 2013) using Eshelby's problem (Eshelby, 1957, 1959) of deformation of a hard inclusion in a softer matrix.

2.6. Numerical experiment setup

The 20 × 10 cm2 initial microstructure with 3264 grains, each with an initially homogeneous, random lattice orientation (Fig. 1c), was deformed up to 70% of shortening in plane-strain pure shear. Dislocation glide of ice single-crystal was defined by slip on the basal plane {0001}{11–20}, prismatic plane {1–100}{11–20} and the pyramidal systems {11–22}{11–2–3}. This approach is similar to that of other authors (Montagnat and others, 2013). In the simulations, the ratio of CRSS for non-basal versus basal slip systems was set to A = 20. This value, <A = 60–100 known for Ih ice, was chosen as a compromise between calculation time (dependent on A) and accuracy of the simulations. Tests showed that the effect of increasing A diminishes strongly for A > 10 (Lebensohn and others, 2007). A constant stress exponent of n = 3 was set for all slip systems. Grain boundary properties were assumed isotropic, i.e. independent of the lattice orientation of the adjacent grains.

The FFT calculation is numerically the most expensive of all. Each deformation step was thus kept at a constant strain increment during calculation. Different strain-rates were achieved by varying the number of steps, each equal to 20 years, in the recrystallisation loop, giving strain-rates from 1.27 × 10−12 to 3.17 × 10−11 s−1 for 1–25 DRX steps per numerical time step, respectively (Table 2). One simulation (Experiment 0) with no DRX is included for reference.

Table 2. Numerical experiment set-up

2.7. Postprocessing

Crystallographic preferred orientation data are presented as pole figures. Furthermore, LPO data were plotted in inverse pole figures, where the sample x-, y- and z-directions are projected with respect to the crystallographic directions. These data were plotted using the texture analysis software MTEX ( (Bachmann and others, 2010; Mainprice and others, 2011) from the orientation distribution function.

The c-axis orientation of the crystal-axis orientation c k can be determined by two angles: the colatitude or tilt angle θ k, and the longitude or azimuth φ k, both in a range [0–2π], considering a reference system RS where the third axis is perpendicular to the xy-plane. c k can expressed as (Durand and others, 2006):

(18) $${c^{\rm k}} = (\cos {\varphi _{\rm k}}\sin {\theta _{\rm k}},{\rm} \sin {\varphi _{\rm k}}\sin {\theta _{\rm k}},{\rm} \cos {\theta _{\rm k}}),$$

The second-order orientation tensor a (2) was used to characterise the c-axis orientation distribution. Following Woodcock (1977), a (2) is defined as:

(19) $${a^{\left( 2 \right)\;}} = \; \mathop \sum \limits_{{\rm k = 1}}^{{N_{\rm g}}} {\,f_{\rm k}}{c^{\rm k}} \otimes {c^{\rm k}},$$

where f k is the volume of a grain k from a 2-D section, that is proportional to its measured cross-sectional area A k and the areas containing non-intersecting crystals N g (Gagliardini and others, 2004). For these simulations, we assume that this volume is proportional to the total number of unodes N p over which the c k is obtained:

(20) $${a^{\left( 2 \right)\;}} = \; \displaystyle{1 \over {{N_{\rm p}}}}\mathop \sum \limits_{{\rm k = 1}}^{{N_{\rm p}}} {c^{\rm k}} \otimes {c^{\rm k}},$$

assuming that a (2) is symmetric, there is a principal RS R sym in which the second-order orientation tensor is diagonal. Then, we can define a i (2) (i = 1, 2, 3) and e i (i = 1, 2, 3) as the expression of the three eigenvalues and their corresponding eigenvectors, associated with the three base vectors of the R sym. The eigenvectors give the axis direction of the ellipsoid that best fits the density distribution of the grain orientations.

In this study, ice crystal preferred orientation symmetry is expressed as the proportion of point (or single maximum), girdle and random components of the [0001] crystallographic axis (Zaffarana and others, 2014). These fabric-type indices were calculated from the three eigenvalues (a 1 (2), a 2 (2), a 3 (2)) as P = a 1 (2)a 2 (2), G = 2a 2 (2) – 2a 3 (2) and R = 3a 3 (2) (Vollmer, 1990).


3.1. Experiment 0: viscoplastic deformation without recrystallisation

The initial microstructure underwent viscoplastic deformation without DRX. The resulting microstructure has irregular grains, elongated parallel to the stretching direction (i.e. sample y-direction). The deformation is localised into conjugate high-strain bands of strongly elongated grains, with an aspect ratio of approximately 15 : 1, surrounded by low-strain areas or ‘microlithons’ (Passchier and Trouw, 2005), where grains are moderately elongated (aspect ratio in the order of 5 : 1) (Fig. 5a). Due to the high differences in CRSS in basal versus non-basal slip systems in ice, the ability to deform a grain depends on the angle ϕ 0 between the c-axis and the compressive stress axis, expressed by the Schmid factor S g (Azuma and Higashi, 1985):

(21) $${S_{\rm g}} = \cos {\phi _{0\;}} \sin {\phi _{0}} .$$

S g is at a maximum when the c-axis is oriented at ϕ 0 = 45° (soft) and at a minimum when ϕ 0 = 0° or 90° (hard) with respect to the compressive axis. The high-strain bands affected soft grains with low resistance to the flow rate (S g). In the absence of recrystallisation all grains are preserved throughout the experiment. Pinching of extremely stretched grains leads to a small increase in the number of grains and concomitant reduction in mean grain size (Table 2). Deformation produced intragranular kinkbands and misorientation, calculated as Kernel Average Misorientation (KAM), was increased in GBs and triple junctions (Fig. 4a and b).

Fig. 4. Kernel Average Misorientation (KAM) map for (a) initial, (b) after viscoplastic deformation and (c) after 25 steps of DRX, showing intragranular heterogeneities for grains 1 and 2. Right column shows the polar figures at the three steps of the simulation, shown as lower hemisphere stereographic projections. The KAM is defined as the average misorientation angle of a given unode with all of its neighbour's unodes.

Fig. 5. Grain network and c-axis orientation (a–d) and boundary misorientation (e–h) at 70% of shortening for simulations with (a, e) no recrystallisation, (b, f) 1 step, (c, g) 10 steps and (d, h) 25 steps of recrystallisation per deformation step. The original size is double that of the images shown. For better visibility images for Experiments 0 and 1 have been enlarged (2x). Crystal orientations are shown as shortening direction y in crystal reference frame (IPF).

From an initially random c-axis distribution (Fig. 6a), LPO evolved into a cone distribution with an opening angle of 90° at 30% shortening (Fig. 6b). C-axes rotated towards the compression axis. The LPO further evolved into a strong single maximum pattern (Fig. 6c, cone opening angle about 60°). In general, the LPO evolves from a random to a point-maximum distribution, with a small girdle component (Fig. 7). Crystallographic orientations that deviate significantly from the point-maximum are mostly found in low-strain regions (Fig. 5a).

Fig. 6. Pole figures of lattice preferred orientation (LPO) for all the simulations performed in this study, showing (a) the initial random distribution, (b) at 30% of shortening, (c) at 70% of shortening. (d) Inverse polar figures of the x-, y- and z-axis of the sample with respect to crystal orientations. Colour bar for polar and inverse figure indicates the multiples of uniform distribution.

Fig. 7. Ice LPO symmetry expressed as the proportion of point (or single maximum), girdle and random components for the c-axes <0001>. Different grey values represent simulations with different amounts of DRX: 0, 1, 10 and 25 GBM steps per deformation step.

3.2. Experiments 1, 10, 25: viscoplastic deformation with recrystallisation

Three different ratios between DRX (GBM and recovery processes) and deformation were modelled (1, 10 or 25 steps of DRX per deformation step), to observe the influence of recrystallisation on the structure evolution. The final grain boundary geometry was interlobate with irregular lobate GBs (Fig. 5f and g) as expected for high temperature deformation (Passchier and Trouw, 2005). The curved shapes of boundaries indicate the direction of GBM such that regions with high dislocation densities were swept away. Despite the amount of recrystallisation, grains remain elongated and parallel to the stretching direction. The main differences caused by DRX in comparison with Experiment 0 are: (1) larger grains, with lower intracrystalline heterogeneities; (2) a more equidimensional microstructure with smooth boundaries; (3) without localisation bands of strongly elongated grains (Fig. 5b–d).

The mean grain area increased, however not proportionally to the number of GBM steps (Table 2). The final mean grain area is 5.3 and 13.4 times larger in Experiment 10 and Experiment 25 compared with Experiment 1. For a linear increase in grain area with time as expected for NGG one could expect factors of 10 and 25. At 70% strain, 41% of the initial grains disappeared in Experiment 1, and only 11 and 5% of grains survived in Experiment 10 and 25, respectively (Table 2). The recovery process produced high-misorientation areas (>15°), resulting in the development of subgrains in Experiments 10 and 25 (Fig. 5g and h). Subgrain formation is controlled by both viscoplastic deformation and intracrystalline recovery. It is a function of plastic rotation rate gradient, while growth and evolution is controlled by misorientation minimisation simulated by the recovery process. With active DRX, the internal boundary misorientation in general first increased, followed by a reduction when grains start to disappear (grey shading in Fig. 5f and h). Misorientation remained high at GBs and triple junctions after DRX (Fig. 4c).

Although activation of recrystallisation causes a remarkable change in grain boundary microstructure, it does not significantly modify the development of a single maximum distribution of c-axes compared with Experiment 0 (Fig. 5b–d). The single-maximum c-axes LPO at low angles to the shortening direction (parallel to the compression axis) develops at medium strains (Fig. 6b), leading to a strong single maximum at 70% strains (Fig. 6c). With increasing DRX the c-axes maxima develop a small angle with respect to the shortening direction. Figure 7 shows that c-axes LPO evolve quicker than in the case of no DRX and first show a point-to-girdle distribution, followed by a decreasing girdle component. The a-axes distribution shows a stronger effect of DRX with the development of two sets of 3 distinct point maxima (3 axes <11–20> and 3 prism plane normal axes {10–10}), distributed around a girdle in the yz-plane (Fig. 6d).


We have applied a FFT formulation to model viscoplastic deformation by dislocation glide, together with GBM and recovery processes, in order to simulate the development of the microstructure and LPO in ice. Results of four experiments that represent deformation at different strain-rates were compared. Our approach serves to predict the viscoplastic behaviour of an isotropic polycrystal on the microscopic scale.

4.1. Microstructure evolution

In Experiment 0 (only deformation) there is no grain growth, as GBM (static and dynamic) is not allowed and, therefore, all grains survive (Table 2). A slight increase in total number of grains is observed at the end of the experiment (Table 2) due to splitting by topological changes of the grain network during simulations. Grain growth is clearly observable in Experiments 1, 10 and 25, with the final grain size increasing with number of DRX steps (Fig. 8a). In all experiments the grain-growth rate initially follows that for NGG (i.e. dominated by boundary energy) (k = 0.04 mm2 a−1), followed by a stage of faster growth (Fig. 8c). This effect is stronger at a high-strain-rate (Experiment 1, k = 0.08 mm2 a−1) than at the lowest train-rate (Experiment 25, k = 0.064 mm2 a−1). The difference in the strain stored energy produces higher bulk effective mobilities by increasing the strain-rate. The increase in grain-growth rate is caused by GBM driven by strain energy, which is the strongest at high-strain-rates where recovery has less time to recover internal strains. It should, however, be borne in mind that grain size reduction by polygonisation is not fully implemented in the model.

Fig. 8. Average grain area evolution for all experiments with respect to the amount of shortening (a) and (b) aspect ratio evolution during deformation showing in grey lines the vertical shortening. The average grain area increases with time, being the effective mobility, higher in Experiment 1 than in Experiment 10 and 25 compared with the mobility for a NGG experiment (c). The “only NGG curve” (NGG experiment) and “only deformation” curve (Experiment 0) are considered the two extreme cases. Experiment 0 curve in (a) is almost parallel to the shortening axis.

The evolution of grain morphology depends on the ratio between deformation and DRX. Deformation without recrystallisation produces elongated grains, while DRX leads to larger and equidimensional grains. In Experiment 0 the average grain area remains constant during deformation, while aspect ratios of grains range from 5 : 1 in low-strain regions (or microlithons) to 15 : 1 in high-strain localisation bands (Fig. 5e). These low values of grain's aspect ratio highly contrast with the theoretical aspect ratio of 11 at 70% of shortening, which thus cannot be due to the highly intrinsic anisotropy of the ice but due to extensive DRX and associated GBM.

Figure 8b shows the influence of the amount of DRX in the relationship between grain size and grain elongation. When DRX is activated, grain elongation is strongly reduced, even in Experiment 1 where the grain aspect ratio only reaches a value of <3 at 70% shortening (Figs 5b and 8b). Shear localisation is difficult to discern in the microstructure in Experiment 1 and completely masked by DRX and associated GBM in the lower strain-rate experiments (Fig. 5c and d).

In pure ice, GBM inhibits the development of strong grain elongations at natural strain-rates. Very elongated grains are only found in polar ice cores with a very high content of impurities, which reduces GB mobility (Fitzpatrick and others, 2014). In the Greenland ice core project (GRIP) (central Greenland) and Vostok (east Antartica) ice core, the aspect ratio of grains observed has an average value of about 1.35 (Thorsteinsson and others, 1995) and 1.5 (Lipenkov and others, 1989), respectively. Slightly higher values (average 1.8 to maximum 2.2; unpublished data) were observed in the EPICA Dronning Maud Land (EDML) core (Antarctica). This suggests that DRX in polar ice is at least as important as in our Experiment 25.

When intense DRX (experiments 10 and 25) is activated, intragranular kink bands, developed by folding of the basal plane due to the highly anisotropy of ice crystal, disappear with increasing strain and grains tend to reduce the internal misorientation (Fig. 4). Therefore, DRX reduces the internal misorientation and sweeps away the kink bands produced by deformation. This observation contrasts with experimental observations and previous predictions by numerical simulations where kink bands are preserved (Lebensohn and others, 2008; Montagnat and others, 2011; Piazolo and others, 2015), mainly due to the limited shortening and high-strain-rate of the experiments, or because DRX processes were not numerically simulated. In contrast, kink band structures are not observed in polar ice cores. A plausible simple explanation of this is that intragranular recovery relaxes these structures and/or evolves to develop subgrain boundaries.

Heterogeneous distribution of stored strain energy and differences in grain curvature are the driving forces to explain the dynamics of growth/shrinkage of grains during the simulations. A summary of surviving grains is shown in Table 2. The number of grains disappearing during deformation increases with increasing DRX, as a result of grain growth. In general, grains oriented with high Schmid factors (Eqn (21)) (i.e. soft orientations) are expected to deform strongly and increase their internal stored strain energy (Thorsteinsson, 2002). Therefore, these soft grains are expected to be overgrown by grains with low Schmid factors (i.e. hard orientations). In our experiments, the survival of grains seems to be only weakly dependent on their initial lattice orientation relative to the stress field (S g): only grains oriented with their c-axes parallel to the xy-plane have a slightly enhanced chance of disappearing compared with grains with other orientations (Fig. 9a–c). The initial Schmid factor thus appears to play a minor role in the survival chances of grains, possibly due to the heterogeneous strain distribution within the polycrystal, which is consistent with findings by Thorsteinsson (2002). Initial grain size, however, is found to be a primary factor in determining grain-survival chances. As GB energy is a dominant driving force for GBM, grain survival probability is highly correlated with the initial grain size (Fig. 9d). This effect is observable in growth curves (Fig. 8c), where all the experiments follow the NGG curve at the beginning of the simulations.

Fig. 9. Initial grain boundary network for Experiment 25, showing (a) Schmid factor (S g) for grains that survive until the end of the experiment. Initial pole figures of (b) all grains, and (c) surviving grains. (d) Initial average grain area distribution and survival rate at the end of Experiment 25. Larger grains have a higher probability of survival, while small grains disappear.

In terms of the numerical model, a larger grain size implies larger regions of homogenous lattice orientation, and therefore, regions where incompatibilities of strain field are expected to be minimised. Therefore, deformation can be distributed more homogenously without developing large strain gradients. In the other extreme, regions of small grains and/or high internal misorientation are associated with high-strain gradients and, therefore, high values of stored strain energy. According to these observations, no correlation between Schmid factor and dislocation density or misorientation values is observed. Although our models are limited because nucleation of new grains, and therefore grain size reduction, was not simulated, our observation can be used to interpret situations where grain coarsening is dominant (e.g. during initial grain coarsening in the upper levels of ice sheets).

4.2. LPO evolution

In all simulations the initially random c-axis orientation distribution of grains is destroyed, and c-axes rotate towards the maximum shortening direction (Figs 5 and 6). This evolution has been observed before in the laboratory experiments by Jacka and Maccagnan (1984), where ice samples deformed at −3°C on compression develop circle girdle fabrics around the compression axis (Castelnau and others, 1996). When DRX is activated the c-axes start to develop an angle with the xy-plane (Fig. 6c and d). The influence of the amount of recrystallisation on the LPO evolution can be observed in Figure 7, where the initial random c-axis orientation distribution is destroyed first, followed by the formation of a combined girdle- and single-maximum fabric and finally tending towards a single maximum pattern. The girdle fabric is most pronounced at a high number of DRX steps, where the random orientations are also destroyed earlier than in experiments with little or no DRX.

In experiments with a large number of recrystallisation steps (Experiments 10 and 25), the a-axes simultaneously rotate, defining a girdle parallel to the xz-plane and maximum parallel to the elongation axis.

The FFT algorithm allows tracking of the relative activity of individual slip systems. Although the LPO development in the different experiments is qualitatively similar (Fig. 7), the amount of DRX causes significant differences in the relative activities of the slip systems (Fig. 10). Basal slip activity decreases with strain, compensated by non-basal slip (Fig. 10). DRX amplifies this behaviour up to a maximum shortening of 70%. The highest amounts of DRX (Experiment 25) exhibit the extreme case: pyramidal slip becomes the most important slip system at ~55% shortening. At higher strains this behaviour may be reverted at differing strain thresholds, indicated by a changing inflection point of the curves in Figure 10.

Fig. 10. Relative activity of basal and non-basal slip systems during deformation for all the simulations, calculated from Eqn (2).

The decrease of basal slip activity and its replacement by pyramidal slip has been described before in simulations using a visco-plastic self-consistent approach in uniaxial compression (Castelnau and others, 1996). However, the rise of pyramidal slip activity with more DRX does not significantly affect the final c-axis distribution, which is similar for all experiments, but instead produces an anisotropic distribution in a- and prism plane normal axes (<11–20> and <10–10>, Fig. 6). Such an anisotropic distribution of a-axes is difficult to observe in ice cores, as full crystal orientation measurements of statistically significant numbers of grains are not standard. In spite of these difficulties, it has been described in one ice core (Miyamoto and others, 2005): the deep part of the GRIP ice core (with a strong single-maximum region of c-axis). A significant concentration of a-axis is observed at certain depths. The anisotropy of slip along different crystallographic directions on the basal plane is proposed as a possible explanation (Miyamoto and others, 2005). Alternatively, the a-axes distribution could indicate a simple shear component to an overall compressional strain regime: a preferred glide direction in the basal plane will then cause alignment of a-axes along the simple shear direction perpendicular to the c-axes. However, as our simulations were performed in pure shear the first explanation proposed by Miyamoto and others (2005) may hold: if a slip system on a prismatic or/and pyramidal plane is active, a-axes may be able to concentrate gradually with depth (Miyamoto and others, 2005).

4.3. Stress evolution

Pure shear simulations are characterised by an initial strain hardening stage followed by a steady state at large strain (Fig. 11d). This hardening is associated with the alignment of c-axes close to the compression axis (Figs 11 and 6). The plateau is reached at lower equivalent stress (Duval and others, 1983) for experiments with more DRX, and is still not reached at 70% shortening in Experiment 1. With more DRX steps the strain hardening stage is contracted, so that simulations reach a steady state at 60% of shortening in Experiment 10 and 50% of shortening in Experiment 25. In Figure 11a–c the evolution of the stress tensors during deformation shows how σ xx and σ yy (extension and compression direction) evolve differently. At the beginning, the out-of-plane deviatoric stress (σ zz) is ~0. As the material hardens, this σ zz increases in absolute magnitude to maintain plane-strain deformation. If we assume that the macroscopic stress exponent can be calculated for the same amount of shortening for simulations at different strain-rates, in order to compare materials of similar states (e.g. similar anisotropy), we can observe an increase of the stress exponent during deformation (Fig. 11e). The exponent increases from the initial n = 3, imposed by the exponent used in the constitutive equation for the individual slip systems (Eqn (2)) to ~3.5 at the end of the experiments. An explanation for this increase of the bulk stress exponent is the DRX-dependent activation of the pyramidal slip system and ensuing strain hardening. It should be noted that the strain hardening is not caused by an imposed hardening law for slip on individual systems, but emerges owing to the development of an LPO, and possibly a shape-preferred orientation (SPO) as well. An SPO with grains elongated perpendicular to the shortening direction has been shown to produce strain hardening (Treagus and Lan, 2004; Takeda and Griera, 2006).

Fig. 11. Deviatoric stress components evolution during deformation for Experiments (a) 1, (b) 10 and (c) 25. Differential stress versus bulk shortening for Experiments 1, 10 and 25 in (d). Differential stresses in (d) were normalised using the experimental values by Duval and others (1983), assuming an average grain size of 1 mm, a temperature of −30°C and the correspondent strain-rate for each experiment. (e) Logarithmic relationship between strain-rate and deviatoric stresses for all simulations performed with DRX in (b). The macroscopic stress exponent n resulting from the first step (i.e. 0.02%), 30 and 70% of shortening are indicated.

4.4. Limitations of the model and further processes

In our simulations nucleation, in the sense of formation of new GBs is not implemented yet, but its influence on the microstructure evolution of our simulations has to be considered for interpretation of natural ice (RRX and SIBM-N) (Faria and others, 2014b). Potential nucleation regions will correspond to areas of heterogeneous intracrystalline deformation and high misorientation. For Experiments 10 and 25, the size of the nucleated grains by RXX (Poirier, 1985; Alley, 1992) will be small (as observed by misorientations in Fig. 5g and h) due to the rapidly increasing average grain size. Stored strain energy is not expected to be a strong enough driving force for the growth of the nuclei, due to low intracrystalline heterogeneities, which will produce new grains that are quickly overgrown by the large neighbouring grains. At high rates of GBM, the effect of spontaneous nucleation can therefore be expected to be minor.

We can assume that the possibility of reaching a grain-size steady state depends on the relative amounts of formation of new grains by splitting and nucleation versus GBM (RXX/SIBM), as explained by Faria and others, (2014b) using observations described for example, by Alley and others (1995). If creation of new grains is controlled by polygonisation (i.e. progressive rotation of subgrains up to development of high misorientation boundaries), we should not expect a strong influence on LPO development, as new grains have lattice orientations that are close to those of their parent grains. However, the numerical implementation of the development of proper high-angle GBs out of tilt walls would reduce both the grain size and elongation compared with the results presented here.

Laboratory creep deformation experiments show a hardening process with a peak strength at a few percentage of strain, followed by a steady state at 10–20%, in either uniaxial compression or simple shear experiments (e.g. Treverrow and others, 2012). Our results differ from this evolution, in that steady state creep is not yet reached at 70% of shortening. Even though our simulations have been performed at constant strain-rate and cannot be directly compared with models with three creep regimes, we can hypothesise that this difference could be caused by the nucleation effect on the microstructure, where a steady state is associated with an equilibrium grain size (Jacka and Li, 1994). When this steady state is reached depends on the equilibrium grain size, relative to the starting grain size and the relative roles of nucleation and polygonisation in controlling the steady state. At the high-strain-rates in experiments and, therefore, relatively slow GBM, nucleation probably plays a larger role than it would do in nature at much lower strain-rates. The absence of nucleation in the current simulations may thus limit their applicability to experimental strain-rates, but not necessarily to natural ones.

Other deformation mechanisms apart from dislocation glide can potentially be active, such as diffusion creep or grain boundary sliding for ice (Duval and others, 2000; Goldsby and Kohlstedt, 2001). As ice in ice sheets generally deforms at rather low-strain-rates, diffusion as an accommodating process cannot be ruled out completely when the grain size is small. However, microstructural evidence for this is scarce, whereas there is abundant evidence for dislocation activity (Faria and others, 2014a). Although the grain size at the beginning of our experiments is small, grain growth makes the assumption of dislocation glide as the dominant deformation mechanism increasingly valid, at least to a first-order.

Spatial resolution of the unodes grid has an influence on the distribution of the predicted strain-rate and stress fields, and therefore has a potential impact on subgrain development inside grains. An increase of the resolution provides more degrees of freedom to adjust local gradients and an improved resolution of intracrystalline heterogeneities. Figure 12 shows a comparison of the predicted misorientation maps for models with 256 × 256 and 512 × 512 unodes. In general, a higher-resolution grid produces sharper and better-defined subgrain boundaries than lower resolution models where heterogeneities are more distributed. Although smaller subgrain bands are visible in the higher-resolution model, both models show qualitatively very similar grain and subgrain sizes. Most subgrains are located near triple junctions and high-angle GBs, where problems maintaining the strain compatibility arise. This observation indicates that subgrain sizes and shapes inside grains are primarily controlled by the local crystallographic orientation and size of the grain and the local evolution of the neighbouring grains, but also clearly by the recrystallisation processes as observed in experiments with more DRX (see Fig. 5). The potential effects that the coupling between viscoplastic deformation and recrystallisation processes may have on length scale and evolution of subgrains need to be addressed in future studies

Fig. 12. Comparison of the predicted misorientation maps for models with (a) 256 and (b) 512 resolution unode grids. A higher resolution produces better defined subgrain boundaries than lower resolution, where heterogeneities are more distributed.

4.5. Implications on interpretation of natural polar ice microstructures

Experiment 0 is an end-member simulation, in the sense that ice at high homologous temperatures (~0.7 to 1) on Earth will not deform without a significant contribution from recrystallisation. This experiment illustrates the effect of formation mechanisms of  LPOs on polycrystals of highly anisotropic ice. Dislocation glide, preferably along basal planes, leads to continuous rotation of crystals towards the bulk compression axis. LPOs are typically the most acknowledged microstructural feature in observations in ice cores and experiments (e.g. Rigsby, 1951; Thorsteinsson and others, 1995; Azuma and others, 1999; Wang and others, 2003; DiPrinzio and others, 2005; Treverrow and others, 2012; Fitzpatrick and others, 2014; Montagnat and others, 2014 more refs in Faria and others, 2014a). The strains reached in the simulations (up 70% shortening) nominally correspond to ice sheet depths down to ~1000 m (e.g. Huybrechts and others, 2007). It is roughly at this depth that a change is observed from weak c-axes LPO to strong LPOs (Faria and others, 2014c). Strong LPO development, thus sets later in ice sheets than in our simulations. It remains to be tested, whether this is because nucleation and polygonisation are not yet fully implemented in our model. The simulations do, however, indicate that GBM is not the cause for this, as LPO-development is very similar with (Experiments 1, 10 and 25) and without GBM (Experiment 0). This is in accordance with the observations in ice cores that LPO is continuously developing with depth while other RX measures (grain shape, subgrain boundary density; Weikusat and others, 2009; Fitzpatrick and others, 2014) do not change systematically with depth.

Recrystallisation is the cause of dynamic grain growth (DGG) in natural ice, usually observed in clean interglacial ice cores (e.g. Faria and others, 2014c). Experiment 0 (Fig. 8a and b) microstructure differs significantly from natural polar ice microstructures, which emphasises the necessity to include recrystallisation processes in consideration of deformation by dislocation glide. Particularly interesting is the question concerning the dominance of deformation versus recrystallisation, which was investigated in Experiments 1, 10 and 25. Grain growth occurs when DRX is activated. Although DGG is not proportional to the amount of GBM, grain size and grain elongation depend on the amount of DRX (Fig. 8b).

The GBM module includes static and dynamic driving forces, which determine the motion of boundaries according to the local energy field from grain boundary curvature and the difference in dislocation density (stored strain energy) across the grain boundary. Previous authors have mentioned that a difference in strain between the considered grains and bulk values is needed to model the LPO evolution in ice (Van der Veen and Whillans, 1994); however the strictly local energy consideration of our simulations is necessary for a realistic GBM, as grain boundaries do not ‘see’ far into the neighbouring grain, as concluded from data on ice core microstructures (Weikusat and others, 2009) and creep deformation samples (Hamann and others, 2007). Local driving forces for GBM lead to increasing grain growth with increasing amounts of RX, which confirms the concept that DGG accounts for ice grain size evolution (Faria and others, 2014b). All microstructural studies on polar ice cores so far report certain amounts of grain growth with depth, thus time and deformation – at least as long as we observe similar types of ice (e.g. ‘clean’ interglacial ice) (e.g. Faria and others, 2014a and references therein). However, nucleation and formation of new GBs are not yet implemented into the suite of simulations presented here and will change the grain size evolution again. Our model is so far able to simulate grain size evolution in the grain growth field of the recrystallisation diagram of Faria and others, (2014b) (below the attractor surface Dss).


Observations of microstructure and LPO evolution during numerical simulations at different strain-rates have been presented, providing new information about the influence of recrystallisation during deformation. DRX can significantly increase grain size, avoiding elongation of grains and masking the deformation and strain localisation. DRX does not modify the c-axes LPO evolution, which in all cases is oriented parallel to the shortening direction. However, in experiments at low-strain-rates, the a-axes rotate towards the maximum elongation axis simultaneously. The anisotropic distribution of a-axes is related to the increment of activity of non-basal slip systems. Although GBM is driven by stored strain and boundary energy reduction, the grains’ survival probability is highly related to the initial grain size, and only slightly related to the crystal orientation with respect to the flow rate. Grain growth in experiments at high-strain-rate is faster than in experiments performed at low-strain-rate, due to the stored strain energy.


We thank Daniela Jansen for discussions and ideas. We thank all the members of the ELLE development group, in particular Sandra Piazolo and Verity Borthwick for their contributions to the simulation code. The study was funded by HGF (VH-NG-802). We thank the editor and two anonymous reviewers for their critical, but very helpful comments and corrections.


Alley, RB (1988) Fabrics in polar ice sheets: development and prediction. Science, 240(4851), 493495 (doi: 10.1126/science.240.4851.493)
Alley, RB (1992) Flow-law hypotheses for ice-sheet modeling. J. Glaciol., 38(129), 245256
Alley, RB, Gow, AJ, Meese, DA (1995) Mapping c-axis fabrics to study physical processes in ice. J. Glaciol., 41(137), 197203
Azuma, N (1994) A flow law for anisotropic ice and its application to ice sheets. Earth Planet. Sci. Lett., 128(3), 601614 (doi: 10.1016/0012-821X(94)90173-2)
Azuma, N and Higashi, A (1985) Formation processes of ice fabric pattern in ice sheets. Ann. Glaciol., 6, 130134
Azuma, N and 6 others (1999) Textures and fabrics in the Dome F (Antarctica) ice core. Ann. Glaciol., 29, 163168 (doi: 10.3189/172756499781821148)
Azuma, N, Miyakoshi, T, Yokoyama, S and Takata, M (2012) Impeding effect of air bubbles on normal grain growth of ice. J. Struct. Geol., 42, 184193 (doi: 10.1016/ j.jsg.2012.05.005)
Ashby, MF (1970) The deformation of plastically non-homogeneous materials. Philos. Mag., 21(170), 399424
Ashby, MF and Duval, P (1985) The creep of polycrystalline ice. Cold Regions Sci. Technol., 11(3), 285300 (doi: 10.1016/ 0165-232X(85)90052-7)
Bachmann, F, Hielscher, R and Schaeben, H (2010) Texture analysis with MTEX – free and open source software toolbox. Solid State Phenom., 160, 6368 (doi: 10.4028/
Becker, R (1991) Analysis of texture evolution in channel die compression. 1. Effects of grain interaction. Acta Metall. Mater., 39(6), 12111230 (doi: 10.1016/0956-7151(91)90209-J)
Becker, JK, Bons, PD and Jessell, MW (2008) A new front-tracking method to model anisotropic grain and phase boundary motion in rocks. Comput. Geosci., 34(3), 201212 (doi: 10.1016/j.cageo.2007.03.013)
Bons, PD, Koehn, D and Jessell, MW (2008) Lecture notes in earth sciences. In Bons, P, Koehn, D and Jessell, M eds. Microdynamic simulation. Springer, Berlin, number 106. 405
Borthwick, VE and Piazolo, S (2010) Complex temperature dependent behaviour revealed by in-situ heating experiments on single crystals of deformed halite: New ways to recognize and evaluate annealing in geological materials. J. Struct. Geol., 32, 982996 (doi: 10.1016/j.jsg.2010.06.006)
Borthwick, VE, Piazolo, S, Evans, L, Griera, A and Bons, PD (2013) What happens to deformed rocks after deformation? A refined model for recovery based on numerical simulations. Geol. Soc. London Spec. Publ., 394(1), 215534 (doi: 10.1144/SP394.11)
Brinckmann, S, Siegmund, T and Huang, Y (2006) A dislocation density based strain gradient model. Int. J. Plast., 22(9), 17841797 (doi: 10.1016/ j.ijplas.2006.01.005)
Budd, W and Jacka, T (1989) A review of ice rheology for ice sheet modelling. Cold Regions Sci. Technol., 16(2), 107144 (doi: 10.1016/0165-232X(89)90014-1)
Budd, WF, Warner, RC, Jacka, TH, Li, J and Treverrow, A (2013) Ice flow relations for stress and strain-rate components from combined shear and compression laboratory experiments. J. Glaciol., 59(214), 374392 (doi: 10.3189/ 2013JoG12J106)
Castelnau, O, Duval, P, Lebensohn, RA and Canova, GR (1996) Viscoplastic modeling of texture development in polycrystalline ice with a self-consistent approach: Comparison with bound estimates. J. Geophys. Res.: Solid Earth (1978–2012) , 101(B6), 1385113868 (doi: 10.1029/96JB00412)
Cottrell, AH (1953) Theory of dislocations. Prog. Met. Phys., 4, 205264 (doi: 10.1016/ 0502-8205(49)90004-0)
Cuffey, KM and Kavanaugh, JL (2011) How nonlinear is the creep deformation of polar ice? A new field assessment. Geology, 39(11), 10271030 (doi: 10.1130/G32259.1)
De La Chapelle, S, Castelnau, O, Lipenkov, V and Duval, P (1998) Dynamic recrystallization and texture development in ice as revealed by the study of deep ice cores in Antarctica and Greenland. J. Geophys. Res.: Solid Earth (1978–2012) , 103(B3), 50915105 (doi: 10.1029/97JB02621)
Delannay, L, Logé, RE, Chastel, Y, Signorelli, JW and Van Houtte, P (2003) Measurement of in-grain orientation gradients by EBSD and comparison with finite element results. Adv. Eng. Mater., 5(8), 597600 (doi: 10.1002/adem.200300376)
DiPrinzio, CL and 5 others (2005) Fabric and texture at Siple Dome, Antarctica. J. Glaciol., 51(173), 281290 (doi: 10.3189/172756505781829359)
Diard, M and 6 others (2006) The heterogeneous nature of slip in ice single crystals deformed under torsion. Philos. Mag., 86(27), 42594270
Diard, O, Leclercq, S, Rousselier, G and Cailletaud, G (2005) Evaluation of finite element based analysis of 3D multicrystalline aggregates plasticity. Application to crystal plasticity model identification and the study of stress and strain fields near grain boundaries. Int. J. Plast., 21(4), 691722 (doi: 10.1016/j.ijplas.2004.05.017)
Drury, MR and Urai, JL (1990) Deformation-related recrystallization processes. Tectonophysics, 172(3), 235253 (doi: 10.1016/0040-1951(90)90033-5)
Durand, G (2004) Microstructure, recristallisation et deformation des glaces polaires de la carotte EPICA. (Thèse de doctorat, Dôme Concordia, Antartique)
Durand, G and 5 others (2006) Ice microstructure and fabric: an up to date approach to measure textures. J. Glaciol., 52(179), 619630 (doi: 10.3189/172756506781828377)
Durand, G, Persson, A, Samyn, D, and Svensson, A (2008) Relation between neighbouring grains in the upper part of the NorthGRIP ice core—Implications for rotation recrystallization. Earth Planet. Sci. Lett., 265(3), 666671 (doi: 10.1016/j.epsl.2007.11.002)
Duval, P and Castelnau, O (1995) Dynamic recrystallization of ice in polar ice sheets. J. Phys. IV, 5(C3), C3197 (doi: 10.1051/jp4:1995317)
Duval, P, Ashby, M and Anderman, I (1983) Rate controlling processes in the creep of polycrystalline ice. J. Phys. Chem., 87(21), 40664074 (doi: 10.1021/ j100244a014)
Duval, P, Arnaud, L, Brissaud, O, Montagnat, M, and de La Chapelle, S (2000) Deformation and recrystallization processes of ice from polar ice sheets. Ann. Glaciol., 30(1), 8387 (doi: 10.3189/172756400781820688)
Erb, U, and Gleiter, H (1979) The effect of temperature on the energy and structure of grain boundaries. Scripta Metall., 13.1, 6164 (doi: 10.1016/0036-9748(79)90390-9)
Eshelby, JD (1957) The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. London, A241, 376396 (doi: 10.1098/rspa.1957.0133)
Eshelby, JD (1959) The elastic field outside an ellipsoidal inclusion. Proc. R. Soc. London, 252(1271), 561569 (doi: 10.1098/rspa.1959.0173)
Estrin, Y and Kim, HS (2007) Modelling microstructure evolution toward ultrafine crystallinity produced by severe plastic deformation. J. Mater. Sci., 42(5), 15121516 (doi: 10.1007/s10853-006-1260-8)
Faria, SH, Weikusat, I and Azuma, N (2014a) The microstructure of polar ice. Part I: highlights from ice core research. J. Struct. Geol., 61, 2149 (doi: 10.1016/j.jsg.2013.09.010)
Faria, SH, Weikusat, I and Azuma, N (2014b) The microstructure of polar ice. Part II: state of the art. J. Struct. Geol, 61, 2149 (doi: 10.1016/j.jsg.2013.11.003)
Faria, SH, Kipfstuhl, S and Lambrecht, A (2014c) The EPICA-DML deep ice core. Springer, Heidelberg
Fitzpatrick, JJ and 11 others (2014) Physical properties of the WAIS divide ice core. J. Glaciol., 60(224), 11811198 (doi: 10.3189/2014JoG14J100)
Fleck, NA and Hutchinson, JW (1997) Strain gradient plasticity. Adv. Appl. Mech., 33, 295361 (doi: 10.1016/S0065-2156(08)70388-0)
Gagliardini, O, Durand, G and Wang, Y (2004) Grain area as a statistical weight for polycrystal constituents. J. Glaciol., 50(168), 8795 (doi: 10.3189/172756504781830349)
Gao, H, Huang, Y, Nix, WD and Hutchinson, JW (1999) Mechanism-based strain gradient plasticity – I. theory. J. Mech. Phys. Solids, 47(6), 12391263 (doi: 10.1016/ S0022-5096(98)00103-3)
Gillet-Chaulet, F, Hindmarsh, RC, Corr, HF, King, EC and Jenkins, A (2011) In-situ quantification of ice rheology and direct measurement of the Raymond Effect at Summit, Greenland using a phase-sensitive radar. Geophys. Res. Lett., 38(24), L24503 (doi: 10.1029/2011GL049843)
Glen, JW (1955) The creep of polycrystalline ice. Proc. R. Soc. A: Math. Phys. Eng. Sci., 228, 519538 (doi: 10.1098/rspa.1955.0066)
Glen, JW (1958) The flow law of ice: a discussion of the assumptions made in glacier theory, their experimental foundations and consequences. IASH Publ., 47, 171183
Goldsby, DL and Kohlstedt, DL (2001) Superplastic deformation of ice: experimental observations. J. Geophys. Res.: Solid Earth (1978–2012) , 106(B6), 1101711030 (doi: 10.1029/2000JB900336)
Gottstein, G and Mecking, H (1985) Recrystallization. Preferred orientation in deformed metals and rocks: an introduction to modern texture analysis, ed. HR Wenk 183218 (doi: 10.1016/B978-0-12-744020-0.50014-6)
Griera, A and 5 others (2011) Strain localization and porphyroclast rotation. Geology, 39, 275278 (doi: 10.1130/G31549.1)
Griera, A and 6 others (2013) Numerical modelling of porphyroclast and porphyroblast rotation in anisotropic rocks. Tectonophysics, 587, 429 (doi: 10.1130/G31549.1)
Gudmundson, P (2004) A unified treatment of strain gradient plasticity. J. Mech. Phys. Solids, 52(6), 13791406 (doi: 10.1016/j.jmps.2003.11.002)
Hamann, I, Weikusat, C, Azuma, N and Kipfstuhl, S (2007) Evolution of ice crystal microstructure during creep experiments. J. Glaciol., 53(182), 479489 (doi: 10.3189/002214307783258341)
Huybrechts, P, Rybak, O, Pattyn, F, Ruth, U and Steinhage, D (2007) Ice thinning, upstream advection, and non-climatic biases for the upper 89% of the EDML ice core from a nested model of the Antarctic ice sheet. Clim. Past Discuss., 3(3), 693727 (doi: 10.5194/cpd-3-693-2007)
Idiart, MI, Moulinec, H, Ponte Castaneda, P and Suquet, P (2006) Macroscopic behavior and field fluctuations in viscoplastic composites: second-order estimates versus full-field simulations. J. Mech. Phys. Solids, 54(5), 10291063 (doi: 10.1016/j.jmps.2005.11.004)
Jacka, TH (1984) Laboratory studies on relationships between ice crystal size and flow rate. Cold Regions Sci. Technol., 10(1), 3142 (doi: 10.1016/0165-232X(84)90031-4)
Jacka, TH and Maccagnan, M (1984) Ice crystallographic and strain rate changes with strain in compression and extension. Cold Regions Sci. Technol., 8(3), 269286 (doi: 10.1016/0165-232X(84)90058-2)
Jacka, TH and Li, J (1994) The steady-state crystal size of deforming ice. Ann. Glaciol., 20, 1328 (doi: 10.3189/172756494794587230)
Jacka, TH and Li, J (2000) Flow rates and crystal orientation fabrics in compression of polycrystalline ice at low temperature and stresses. In Hondoh T ed. Physics of Ice Core Records, Hokkaido University Press, Sapporo. 83113
Jessell, MW, Bons, PD, Evans, L, Barr, T and Stüwe, K (2001a) Elle: the numerical simulation of metamorphic and deformation microstructures. Comput. Geosci., 27, 1730 (doi: 10.1016/S0098-3004(00)00061-3)
Jessell, MW, Evans, L, Barr, T and Stüwe, K (2001b) Modeling of anisotropic grain growth in minerals. Tectonic modeling: a volume in honor of Hans Ramberg, eds A Koyi and NS Mancktelow. 193, 39
Jessell, MW, Kostenko, O and Jamtveit, B (2003) The preservation potential of microstructures during static grain growth. J. Metamorphic Geol., 21(5), 481491 (doi: 10.1046/j.1525-1314.2003.00455.x)
Jessell, MW, Siebert, E, Bons, PD, Evans, L and Piazolo, S (2005) A new type of numerical experiment on the spatial and temporal patterns of localization of deformation in a material with a coupling of grain size and rheology. Earth Planet. Sci. Lett., 239(3), 309326 (doi: 10.1016/j.epsl.2005.03.030)
Jessell, MW, Bons, PD, Griera, A, Evans, LA and Wilson, CJL (2009) A tale of two viscosities. J. Struct. Geol., 31, 719736 (doi: 10.1016/j.jsg.2009.04.010)
Karato, SI (2008) Deformation of earth materials: an introduction to the rheology of solid earth. Cambridge University Press, Cambridge (doi: 10.1017/CBO9780511804892)
Kennicutt, II and 6 others (2014) Six priorities for Antarctic science. Nature, 512(7512), 2325 (doi: 10.1038/512023a)
Ketcham, WM and Hobbs, PV (1969) An experimental determination of the surface energies of ice. Philos. Mag., 19(162), 11611173 (doi: 10.1080/14786436908228641)
Kipfstuhl, S and 6 others (2006) Microstructure mapping: a new method for imaging deformation-induced microstructural features of ice on the grain scale. J. Glaciol., 52(178), 398406 (doi: 10.3189/172756506781828647)
Kipfstuhl, S and 6 others (2009) Evidence of dynamic recrystallization in polar firn. J. Geophys. Res.: Solid Earth (1978–2012) , 114(B5) (doi: 10.1029/2008JB005583)
Kocks, UF, Tomé, CN and Wenk, HR (2000) Texture and anisotropy: preferred orientations in polycrystals and their effect on materials properties. Cambridge University Press, Cambridge
Kröner, E (1962) Dislocations and continuum mechanics. Appl. Mech. Rev., 15(8), 599606
Lebensohn, RA (2001) N-site modelling of a 3D viscoplastic polycrystal using fast Fourier transform. Acta Mater., 49(14), 27232737
Lebensohn, RA and Tomé, CN (1993) A self-consistent approach for the simulation of plastic deformation and texture development of polycrystals: application to Zirconium alloys. Acta Metall. Mater., 41(9), 26112624 (doi: 10.1016/0956-7151(93)90130-K)
Lebensohn, RA, Liu, Y and Castañeda, PP (2004) Macroscopic properties and field fluctuations in model power-law polycrystals: full-field solutions versus self-consistent estimates. Proc. R. Soc. London. Ser. A: Math. Phys. Eng. Sci., 460(2045), 13811405 (doi: 10.1098/rspa.2003.1212)
Lebensohn, RA, Castelnau, O, Brenner, R and Gilormini, P (2005) Study of the antiplane deformation of linear 2-D polycrystals with different microstructures. Int. J. Solids Struct., 42(20), 54415459 (doi: 10.1016/j.ijsolstr.2005.02.051)
Lebensohn, RA, Tomé, CN and Ponte Castañeda, P (2007) Self-consistent modeling of the mechanical behavior of viscoplastic polycrystals incorporating intragranular field fluctuations. Philos. Mag., 87(28), 42874322 (doi: 10.1080/14786430701432619)
Lebensohn, RA, Brenner, R, Castelnau, O, Rollett, AD (2008) Orientation image-based micromechanical modelling of subgrain texture evolution in polycrystalline copper. Acta Mater., 56(15), 39143926 (doi: 10.1016/j.actamat.2008.04.016)
Li, JCM (1962) Possibility of subgrain rotation during recrystallization. J. Appl. Phys., 33(10), 29582965 (doi: 10.1063/1.1728543)
Lipenkov, VY, Barkov, NI, Dural, P and Pimenta, P (1989) Crystalline texture of the 2083 m ice core at Vostok Station, Antarctica. J. Glaciol., 35, 392398
Llorens, M-G, Bons, PD, Griera, A, Gomez-Rivas, E and Evans, LA (2013a) Single layer folding in simple shear. J. Struct. Geol., 50, 209220 (doi: 10.1016/j.jsg.2012.04.002)
Llorens, M-G, Bons, PD, Griera, A and Gomez-Rivas, E (2013b) When do folds unfold during progressive shear? Geology, 41, 563566 (doi: 10.1130/G33973.1)
Ma, A, Roters, F and Raabe, D (2006) A dislocation density based constitutive model for crystal plasticity FEM including geometrically necessary dislocations. Acta Mater., 54.8, 21692179 (doi: 10.1016/j.actamat.2006.01.005)
Mathiesen, J and 6 others (2004) Dynamics of crystal formation in the Greenland NorthGRIP ice core. J. Glaciol., 50(170), 325328 (doi: 10.3189/172756504781829873)
McQueen, HJ and Evangelista, E (1988) Substructures in aluminium from dynamic and static recovery. Czech. J. Phys. B, 38(4), 359372 (doi: 10.1007/BF01605405)
Mainprice, D, Hielscher, R and Schaeben, H (2011) Calculating anisotropic physical properties from texture data using the MTEX open source package. Geol. Soc. London Spec. Publ., 360(1), 175192 (doi: 10.1144/SP360.10)
Michel, JC, Moulinec, H and Suquet, P (2000) A computational method based on augmented Lagrangians and fast Fourier transforms for composites with high contrast. Comput. Model. Eng. Sci., 1(2), 7988
Michel, JC, Moulinec, H and Suquet, P (2001) A computational scheme for linear and non-linear composites with arbitrary phase contrast. Int. J. Numer. Methods Eng., 52(1–2), 139160
Miyamoto, A and 5 others (2005) Ice fabric evolution process understood from anisotropic distribution of a-axis orientation on the GRIP (Greenland) ice core. Ann. Glaciol., 42(1), 4752 (doi: 10.3189/172756405781812501)
Mika, DP and Dawson, PR (1998) Effects of grain interaction on deformation in polycrystals. Mater. Sci. Eng. A, 257(1), 6276 (doi: 10.1016/S0921-5093(98)00824-7)
Mohamed, G and Bacroix, B (2000) Role of stored energy in static recrystallization of cold rolled copper single and multicrystals. Acta Mater., 48(13), 32953302 (doi: 10.1016/S1359-6454(00)00155-5)
Moldovan, D, Wolf, D and Phillpot, SR (2001) Theory of diffusion-accommodated grain rotation in columnar polycrystalline microstructures. Acta Mater., 49(17), 35213532 (doi: 10.1016/S1359-6454(01)00240-3)
Moldovan, D, Wolf, D, Phillpot, SR and Haslam, AJ (2002) Role of grain rotation during grain growth in a columnar microstructure by mesoscale simulation. Acta Mater., 50(13), 33973414 (doi: 10.1016/S1359-6454(02)00153-2)
Montagnat, M and Duval, P (2004) The viscoplastic behaviour of ice in polar ice sheets: experimental results and modelling. Comptes Rendus Physique, 5(7), 699708 (doi: 10.1016/j.crhy.2004.06.002)
Montagnat, M, Weiss, J, Chevy, J, Duval, P, Brunjail, H, Bastie, P and Gil Sevillano, J (2006) The heterogeneous nature of slip in ice single crystals deformed under torsion. Philos. Mag., 86(27), 42594270
Montagnat, M, Blackford, JR, Piazolo, S, Arnaud, L and Lebensohn, RA (2011) Measurements and full-field predictions of deformation heterogeneities in ice. Earth Planet. Sci. Lett., 305(1–2), 153160 (doi: 10.1016/j.epsl.2011.02.050)
Montagnat, M and 11 others (2013) Multiscale modeling of ice deformation behavior. J. Struct. Geol., 61, 78108 (doi: 10.1016/j.jsg.2013.05.002)
Montagnat, M and 9 others (2014) Fabric along the NEEM ice core, Greenland, and its comparison with GRIP and NGRIP ice cores. Cryosphere, 8(4), 11291138 (doi: 10.5194/tc-8-1129-2014)
Moulinec, H and Suquet, P (1994) A fast numerical method for computing the linear and nonlinear mechanical properties of composites. Comptes rendus de l'Académie des sciences. Série II, Mécanique, Physique, Chimie, Astronomie, 318(11), 14171423
Moulinec, H and Suquet, P (1998) Numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput. Methods Appl. Mech. Eng., 157(1), 6994 (doi: 10.1016/S0045-7825(97)00218-1)
Nasello, OB, Di Prinzio, CL and Guzmán, PG (2005) Temperature dependence of “pure” ice grain boundary mobility. Acta Mater., 53(18), 48634869 (doi: 10.1016/j.actamat.2005.06.022)
Nye, J (1953) Some geometrical relations in dislocated crystals. Acta Mater., 1(2), 153162 (doi: 10.1016/0001-6160(53)90054-6)
Passchier, CW and Trouw, RAJ (2005) Deformation mechanisms. Microtectonics, Springer, Berlin. 2566 (doi: 10.1007/978-3-662-08734-3_3)
Piazolo, S, Bons, PD, Jessell, MW, Evans, L and Passchier, CW (2002) Dominance of microstructural processes and their effect on microstructural development: insights from numerical modelling of dynamic recrystallization. Geol. Soc. London Spec. Public., 200(1), 149170 (doi: 10.1144/GSL.SP.2001.200.01.10)
Piazolo, S, Jessell, MW, Bons, PD, Evans, L and Becker, JK (2010) Numerical simulations of microstructures using the Elle platform: a modern research and teaching tool. J. Geol. Soc. India, 75(1), 110127 (doi: 10.1007/s12594-010-0028-6)
Piazolo, S and 6 others (2012) Substructure dynamics in crystalline materials: new insight from in situ experiments, detailed EBSD analysis of experimental and natural samples and numerical modelling. Mater. Sci. Forum, 715, 502507 (doi: 10.4028/
Piazolo, S, Montagnat, M, Grennerat, F, Moulinec, H and Wheeler, J (2015) Effect of local stress heterogeneities on dislocation fields: examples from transient creep in polycrystalline ice. Acta Materialia, 90, 303309
Prakash, A and Lebensohn, RA (2009) Simulation of micromechanical behavior of polycrystals: finite elements versus fast Fourier transforms. Model. Simul. Mater. Sci. Eng., 17.6, 064010 (doi: 10.1088/0965-0393/17/6/064010)
Poirier, JP (1985) Creep of crystals. Cambridge University Press, New York (doi: 10.1017/CBO9780511564451)
Raabe, D, Sachtleber, M, Zhao, Z, Roters, F and Zaefferer, S (2001) Micromechanical and macromechanical effects in grain scale polycrystal plasticity experimentation and simulation. Acta Mater., 49(17), 34333441 (doi: 10.1016/S1359-6454(01)00242-7)
Randle, V (1993). Microtexture investigation of the relationship between strain and anomalous grain growth. Philos. Mag. A, 67.6, 13011313 (doi: 10.1080/01418619308225356)
Read, WT and Shockley, W (1950) Dislocation models of crystal grain boundaries. Phys. Rev., 78.3, 275 (doi: 10.1103/PhysRev.78.275)
Rigsby, GP (1951) Crystal fabric studies on emmons glacier mount rainier, Washington. J. Geol., 59, 590598 (doi: 10.1086/625914)
Roessiger, J, Bons, PD and Faria, SH (2012) Influence of bubbles on grain growth in ice. J. Struct. Geol., 61, 123132 (doi: 10.1016/j.jsg.2012.11.003)
Roessiger, J and 8 others (2011) Competition between grain growth and grain size reduction in polar ice. J. Glaciol., 57(205), 942948 (doi: 10.3189/002214311798043690)
Roters, F, Eisenlohr, P, Bieler, TR and Raabe, D (2011) Crystal plasticity finite element methods: in materials science and engineering. John Wiley & Sons (doi: 10.1002/9783527631483)
Roters, F and 5 others (2010) Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: theory, experiments, applications. Acta Mater., 58(4), 11521211 (doi: 10.1016/j.actamat.2009.10.058)
Sachs, G (1928) On the derivation of a condition of flowing. Zeitschrift des Vereins Deutscher Ingenieure für Maschinenbau und Metallbearbeitung, 72, 734736
Schulson, EM and Duval, P (2009) Creep and fracture of ice. Cambridge University Press, Cambridge. (p. 432) (doi: 10.1017/CBO9780511581397)
Sellars, CM (1978) Recrystallization of metals during hot deformation. Philos. Trans. R. Soc. London A, 288(1350), 147158 (doi: 10.1098/rsta.1978.0010)
Stipp, M and Tullis, J (2003) The recrystallized grain size piezometer for quartz. Geophys. Res. Lett., 30(21), 2088 (doi: 10.1029/2003GL018444)
Takeda, YT, and Griera, A (2006) Rheological and kinematical responses to flow of two-phase rocks. Tectonophysics, 427.1, 95113 (doi: 10.1016/j.tecto.2006.03.050)
Taylor, GI (1938) Plastic strain in metals. J. Inst. Met., 62, 307324
Thorsteinsson, T (2002) Fabric development with nearest-nieghbour interaction and dynamic recrystallization. J. Geophys. Res., 107(B1), 113 (doi: 10.1029/2001JB000244).
Thorsteinsson, T, Kipfstuhl, J, Eicken, H, Johnsen, SJ and Fuhrer, K (1995) Crystal size variations in Eemian-age ice from the GRIP ice core, central Greenland. Earth Planet. Sci. Lett., 131(3), 381394 (doi: 10.1016/0012-821X(95)00031-7)
Treagus, SH and Lan, L (2004) Deformation of square objects and boudins. J. Struct. Geol., 26.8, 13611376 (doi: 10.1016/j.jsg.2003.12.002)
Treverrow, A, Budd, WF, Jacka, TH and Warner, RC (2012) The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain rate enhancement. J. Glaciol., 58(208), 301314 (doi: 10.3189/2012JoG11J149)
Tullis, J, Dell'Angelo, L and Yund, RA (1990) Ductile shear zones from brittle precursors in feldspathic rocks: the role of dynamic recrystallization. In Ramberg, H, Hemin, AK and Mancktelow, NS eds The brittle-ductile transition in rocks, 6781 (doi: 10.1029/GM056p0067)
Upmanyu, M, Srolovitz, DJ, Shvindlerman, LS and Gottstein, G (2002) Molecular dynamics simulation of triple junction migration. Acta Mater., 50(6), 14051420 (doi: 10.1016/S1359-6454(01)00446-3)
Urai, JL, Means, WD and Lister, GS (1986) Dynamic recrystallization of minerals. In Hobbs, BE ed. Mineral and rock deformation: laboratory studies: The Paterson volume, 161199 (doi: 10.1029/GM036p0161)
Van der Veen, CJ and Whillans, IM (1994). Development of fabric in ice. Cold Regions Sci. Technol., 22(2), 171195 (doi: 10.1016/0165-232X(94)90027-2)
Vollmer, FW (1990) An application of eigenvalue methods to structural domain analysis. Geol. Soc. Am. Bull., 102(6), 786791 (doi: 10.1130/0016-7606(1990)102%3C0786:AAOEMT%3E2.3.CO;2)
Wang, Y, Kipfstuhl, S, Azuma, N, Thorsteinsson, T and Miller, H (2003) Ice-fabrics study in the upper 1500 m of the Dome C (East Antarctica) deep ice core. Ann. Glaciol., 37(1), 97104 (doi: 10.3189/172756403781816031)
Weikusat, I, Kipfstuhl, S, Faria, SH, Azuma, N and Miyamoto, A (2009) Subgrain boundaries and related microstructural features in EDML (Antarctica) deep ice core. J. Glaciol., 55(191), 461472 (doi: 10.3189/002214309788816614)
Weiss, J and Marsan, D (2003) Three-dimensional mapping of dislocation avalanches: clustering and space/time coupling (2003). Science, 299, 89 (doi: 10.1126/science.1079312)
Wen, J, Huang, Y, Hwang, KC, Liu, C and Li, M (2005) The modified Gurson model accounting for the void size effect. Int. J. Plast., 21(2), 381395 (doi: 10.1016/j.ijplas.2004.01.004)
Wilson, CJL (1986) Deformation induced recrystallization of ice; the application of in situ experiments. In Hobbs, BE ed. Mineral and rock deformation: laboratory studies: The Paterson volume, 213232 (doi: 10.1029/GM036p0213)
Wilson, CJL and Zhang, Y (1996) Development of microstructure in the high-temperature deformation of ice. Ann. Glaciol., 23, 293302
Wilson, CJL, Peternell, M, Piazolo, S and Luzin, V (2013) Microstructure and fabric development in ice: lessons learned from in situ experiments and implications for understanding rock evolution. J. Struct. Geol., 61, 5077 (doi: 10.1016/j.jsg.2013.05.006)
Woodcock, NH (1977) Specification of fabric shapes using an eigenvalues method. Geol. Soc. Am. Bull., 88(9), 12311236 (doi: 10.1130/0016–7606(1977)88%3C1231:SOFSUA%3E2.0.CO;2)
Zaffarana, C, Tommasi, A, Vauchez, A and Grégoire, M (2014) Microstructures and seismic properties of south Patagonian mantle xenoliths (Gobernador Gregores and Pali Aike). Tectonophysiscs, 621, 175197 (doi: 10.1016/j.tecto.2014.02.017)