Introduction
Crystal size is a key textural parameter of glacier ice that can influence its microstructural and bulk properties and carry information about its deformation history. In polar ice cores, depth profiles of mean crystal size are routinely measured and studied alongside climatic records and changes in ice fabric (e.g. c –axis orientation). In several deep cores, including those recovered at the GRIP, NGRIP and GISP2 sites in Greenland, Byrd Station in West Antarctica and the summit of Law Dome in East Antarctica, these profiles show that crystal size increases through the top few hundred metres; below this, it grows more slowly to reach approximately steady values in the middle part of the core, before rising again (sometimes after a dip) to large values near the ice-sheet base (Fig. 1).
Measured depth profiles of mean crystal size in five polar ice cores. (a) Crystal width and height in the Greenland Ice Core Project (GRIP) core, from Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others (1997); (b) crystal width, height and vertical area in the GRIP core from Reference Svensson, Durand, Mathiesen, Persson, Dahl-Jensen and HondohSvensson and others (2009); (c) crystal width, height and vertical area in the NorthGRIP (NGRIP) core, Greenland, from Reference SvenssonSvensson and others (2003); (d) crystal width in the Greenland Ice Sheet Project 2 (GISP2) core, measured separately by Reference Alley and WoodsAlley and Woods (1996) and Reference GowGow and others (1997); (e) horizontal crystal area in the Byrd core, West Antarctica, from Gow and Williamson (1976); (f) horizontal and vertical crystal areas in the Dome Summit South (DSS) core at Law Dome, East Antarctica, from Reference Li, Jacka and MorganLi and others (1998). In each profile, each point plots the mean crystal size derived from measurements on a thin section, and ‘S’ marks the region of steady crystal size identified by us or from the original authors' description (if available). The methods of measurement of crystal sizes are given in Table 1.
Research of past decades until several years ago (e.g. Reference AlleyAlley, 1992; Reference Alley, Gow and MeeseAlley and others, 1995; Reference GowGow and others, 1997; Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997; De La Reference De La Chapelle, Castelnau, Lipenkov and DuvalChapelle and others, 1998; Reference Duval and HondohDuval, 2000; Reference SvenssonSvensson and others, 2003; Reference DurandDurand and others, 2006a; also summary by Reference Cuffey and PatersonCuffey and Paterson, 2010) has led to the hypothesis that the crystal-size profiles are shaped by three main mechanisms, as follows. At shallow depths, the dominant mechanism is temperature-controlled normal grain growth, in which grain-boundary migration causes larger crystals to grow at the expense of smaller crystals, so that the mean crystal size increases with time (age) (e.g. Reference GowGow, 1969). Deeper down, while normal grain growth continues to operate, stored strain energy in the crystals causes them to split by polygonization at an increasing rate. Also called ‘rotation recrystallization’, this mechanism starts occurring in the shallow subsurface (Reference Durand, Persson, Samyn and SvenssonDurand and others, 2008) but is apparently important at intermediate depths, where its effect in reducing the mean crystal size is thought to balance normal grain growth to cause the region of steady grain size in the profiles (Reference Alley, Gow and MeeseAlley and others, 1995; Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997). At still greater depths, where the ice temperature exceeds about –10°C, a third mechanism, known as migration recrystallization, becomes dominant (Reference Duval and CastelnauDuval and Castelnau, 1995). Driven by release of stored strain energy, it causes accelerated grain-boundary motion and possibly nucleation of new grains (e.g. Reference Hamann, Weikusat, Azuma and KipfstuhlHamann and others, 2007) and is considered responsible for the dramatic increase in grain size on the profiles near the bed. Note that much of the understanding of these mechanisms has come from interpretation of data on crystal fabric as well as size, and profiles shaped very differently from those in Figure 1 have been measured from the Vostok and EPICA ice cores. The low strain rate and temperature in these cores keep recrystallization rates low, so that crystal size is influenced strongly by impurities (bubbles, solutes, micro-particles) in the ice, which leave their imprint on the profiles over a timescale of several glacial–interglacial cycles (Reference DurandDurand and others, 2006a).
This three-stage or ‘tripartite’ hypothesis has recently been criticized from several angles, as reviewed by Reference Faria, Weikusat and AzumaFaria and others (2014). First, the ‘classical’ parabolic law found for the rate of normal grain growth in field samples of ice and firn (Eqn (1) below) – and used in the hypothesis to describe the growth mechanism – turns out to be different from the normal grain growth law established for pure bubble-free ice samples, which shows a much larger growth rate and steeper temperature dependence (Reference Azuma, Miyakoshi, Yokoyama and TakataAzuma and others, 2012). Hence the classical law describes not pure normal grain growth (coarsening driven by grain boundary energy alone), but its outcome complicated by the retardation effects of impurities, which the law does not resolve. Numerical simulations also suggest that the growth rate of deforming crystals depends on their evolving microstructure and size distribution (Reference RoessigerRoessiger and others, 2011), echoing earlier simulations in which grain flattening is found to alter the rate (Reference Bons and UraiBons and Urai, 1992); nor does the classical law distinguish these factors. A second criticism concerns migration recrystallization. High-resolution microstructural analysis of the EPICA Dronning Maud Land (EDML) ice core has shown that strain-induced grain boundary migration occurs in the firn (Reference KipfstuhlKipfstuhl and others, 2009) and throughout the ice at all depths (Reference Weikusat, Kipfstuhl, Faria, Azuma and MiyamotoWeikusat and others, 2009) in this core. This means that the classical growth law may also be corrupted by signals of migration recrystallization. For the tripartite hypothesis, these realizations oppose its ideas that (1) pure normal grain growth occurs at shallow depths, with minor influence by impurity and strain-induced effects, and (2) migration recrystallization is limited to the lowest part of the ice column. Establishing a sound physical model of migration recrystallization for ice thus seems to be essential; this is challenging because it can affect crystal size and shape in complex ways (e.g. Reference Hamann, Weikusat, Azuma and KipfstuhlHamann and others, 2007).
In this paper we enrich this picture by exploring the transition from growing to steady crystal size in the upper part of the profiles in Figure 1, down to the bottom of their steady regions (labelled ‘S’). Our main goal is to highlight the roles of crystal dislocation density and non-equiaxed crystal shape in the evolution. We do this with a mathematical model formulated to describe the competing effects of grain-boundary migration driven crystal growth and polygonization on the mean crystal size, and we test this model by simulating the profiles. Matching the model to different datasets allows us to estimate the polygonization rates in the five cores. We learn that the problem is not addressable by a single differential equation for the rate of change of crystal size: the dislocation density must also be tracked, because these variables are coupled.
Our model requires, for one of its competing terms, a description of crystal growth. As we do not know how to quantify the rate at which migration recrystallization affects this growth via strain-induced effects (accelerated grain-boundary motion and nucleation), we adopt the classical parabolic law as an empirical descriptor of this growth; the reasons and limitations of this choice are elaborated below. Therefore we emphasize that what we study is a simplistic ‘toy’ model that does not resolve the effects of migration recrystallization and impurities and other potentially important processes (which may, for instance, cause the rapid variations on the profiles), and that our focus is on model development: we are not proposing a universal model of crystal-size evolution and dynamic recrystallization in ice. Using the classical law raises a caveat for our numerical results, but does not upset our conclusion about the coupling, nor preclude future model adaptations. When formulating the model, we consider some possibilities and difficulties of including migration recrystallization in it.
In this connection, sophisticated numerical simulations accounting for diverse recrystallization processes are now used to study polar ice deformation at the crystal scale (review by Reference MontagnatMontagnat and others, 2014). A second goal of our paper is to offer benchmark results showing how well or badly a simple model can be applied to multiple ice cores. This might aid research seeking analytical understanding of these simulations, comparing the properties of different cores or assessing models of different complexity.
Previous authors have addressed various aspects of the observed growth-to-steady transition in the crystal-size profiles. Our model uses their insights but is formulated in an original way to suit our purpose. Reference MathiesenMathiesen and others (2004) proposed a differential-equation model for simulating the changing size distribution of crystal populations across the transition, in which grain growth is represented as a diffusion process and polygonization as a transfer of number density between grains of different sizes; their model yielded an excellent match with the crystal sizes in the NGRIP core down to 1000 m depth. Here we model the mean size of crystals rather than their size distributions, because size distribution data are lacking for some of the cores. Our model brings together three concepts for the first time: (1) By tracking the mean height and width of crystals with separate equations, it describes their flattening by deformation in the ice column. (2) It calculates the rate of grain-boundary migration driven growth as being direction-dependent, as is necessary given (1). (3) It relates the polygonization rate of crystals to the mean dislocation density within them, and this leads to equations for the coupled evolution of crystal size and dislocation density. For this part of the model, we draw on Reference Montagnat and DuvalMontagnat and Duval’s (2000) instrumental work on dislocation dynamics, which elucidates one way of the coupling (how crystal size influences changing dislocation density); we contribute a description of the reverse coupling. As explained later, such two-way coupling features in an earlier model by Reference DurandDurand and others (2006a) for the crystal-size profile of the EPICA ice core, but we make different assumptions from theirs.
Because the ice in the depth ranges of interest mostly dates to the Holocene period (Table 1), we solve our model in a forward sense with constant temperature and strain rate as inputs and do not reconstruct the history of these variables at the core sites.
Information on the five ice cores studied in this paper, their crystal-size measurements, and summary data from their regions of steady crystal size (S in Fig. 1). Indices (i)_(iv) identify repeated crystal-size datasets from the GRIP and GISP2 cores with their different authors and measurement methods
The Model
We consider how crystal (grain) size varies with depth z in the ice sheet, specifically beneath ice-dome centres, where the ice experiences uniaxial vertical compression at strain rate 
(defined to be positive when compressive). Our model may be modified for non-axisymmetric situations at ice divides and flanks; this is not pursued currently. We formulate equations with age t as the independent variable and convert t to z later. Two key variables are the mean height Dz
and mean width Dx
of the crystals, where x denotes the radial direction. Thus the aspect ratio a = Dx
/Dz
quantifies the degree of crystal flattening. Also, by idealizing the shape of the crystals as ellipsoidal, we define their mean horizontal area as 
and mean vertical area as 
these are the areas that result from projecting the crystals on the horizontal and vertical planes, respectively.
Three processes interact in the model: (1) grain growth, driven by grain boundary migration arising from the tendency for crystals to lower their free energy, (2) geometric deformation of the crystals under bulk vertical compression and (3) polygonization, the splitting of crystals when deformation causes dislocations in them to form subgrain boundaries that evolve into grain boundaries; as said before, polygonization reduces the mean crystal size. We consider each process in turn and build it into the model.
‘Classical’ grain growth law
In many earlier descriptions, it is supposed that normal grain growth causes the overall (non-directional) crystal size D to increase with age t according to
Or d (D2)/dt = K, where the growth rate
depends on temperature T (K), and R is the gas constant (8.314 J K–1 mol–1). For constant K, integration gives D 2 = Kt + constant. Best fit of this parabolic growth law to grain-size data from near-surface firn and ice at polar sites with different temperatures (e.g. Reference GowGow, 1969; Reference DuvalDuval, 1985; Reference Cuffey and PatersonCuffey and Paterson, 2010) yields the activation energy Q ≈ 40–60 kJ mol–1 and the rate constant K 0 ≈ 107 mm2 a–1 (the precise values used in our simulations are discussed later). For convenience we call the law with these parameters the ‘classical’ law. In contrast, through laboratory experiments, Reference Azuma, Miyakoshi, Yokoyama and TakataAzuma and others (2012) determined that with the same equations and for pure bubble-free ice, the growth rate K is several orders of magnitude larger and Q = 110–120 kJ mol–1 (fig. 5 of their paper).
We choose the classical law for formulating our model for the following reasons. Ice cores do not contain pure ice; by starting with the pure normal growth rate, one would need to quantify how impurities alter it. Although models exist for the drag effect of bubbles and micro-particles on grain boundary mobility (e.g. Reference Alley, Perepezko and BentleyAlley and others, 1986; Reference DurandDurand and others, 2006a; Reference Azuma, Miyakoshi, Yokoyama and TakataAzuma and others, 2012), their use requires suitable records of these factors for the cores, which we do not have here. Crucially we still lack a model for the effect of migration recrystallization on the right-hand side of Eqn (1) (we have more to say about this below). Thus current knowledge simply does not enable us to write a fully explanatory model for the grain-boundary migration driven growth rate. The alternative starting point – the classical law – is plausible because it describes the growth rates measured at numerous field sites. (Explaining why these rates deviate from the pure rate is important, but not our goal.) To avoid confusion with the pure process, in this paper we refer to the process modelled by the classical law as ‘classical grain growth’.
This choice and its level of empiricism have shortcomings. We cannot treat migration recrystallization nor impurities in our model; nor will we try to predict the effects of crystal microstructure and dislocation density on K below. An intrinsic assumption in our model simulations is that the growth rate, which may already encapsulate such effects via the classical parameters, does not change significantly with depth (the simulations assume constant temperature, so constant K). How valid this is can only be gauged once we can reliably model all such effects on crystal size. This part of the model awaits improvement, but does not compromise our main insights and conclusion.
In order to define our modelling approach further, we draw a strong analogy between our empirical description of crystal growth and the parameterization of basal sliding in glacier flow models, where simple (e.g. power) sliding laws are assumed to allow different aspects of glacier and ice-sheet motion to be studied, when the governing processes of sliding (e.g. subglacial water and sediment dynamics) are incompletely understood and still being researched. The parameters in such laws are also constrained by observations.
Anisotropic grain growth
When the crystals are non-equiaxed (a ≠ 1), we expect them to coarsen in the vertical and radial directions at different rates, because the migration velocities of grain boundaries depend on their curvature. Herein, ‘anisotropic’ grain growth refers to this directional phenomenon, not to the dependence of the motion of grain boundaries on the relative orientation between them and crystal lattices.
We now seek corresponding growth laws where the directional coarsening rates are functions of Dz and Dx . Unfortunately, the theory of Reference HillertHillert (1965) explaining the parabolic growth in Eqn (1) considers the overall grain size and, as far as we know, has not been reformulated for anisotropic grain shape. Lacking a firm theoretical basis here, we motivate a model for the coarsening rates by using recent computational evidence of Reference Di Prinzio and NaselloDi Prinzio and Nasello (2011). In their three-dimensional (3-D) Monte Carlo simulation of flattened grains experiencing normal grain growth, these authors observed that the overall size D still obeys Eqn (1), and that Dz (the size along the axis of flattening) obeys an equation of the same form:
They went on to build a model for the changing elongation ratio of crystals, with the aim of using this ratio to date the ice.
Here we adopt Eqn (3) for Dz and use it with Eqn (1) to derive a separate evolution equation for the mean crystal width Dx. Define the overall grain size as the geometric mean of the two directional grain sizes:
then we are
Substituting for the rates dD /dt and dDz /dt from Eqns (1) and (3), and casting the result in terms of Dz and Dx , yields
Equivalently we can write
in which the geometric factor
measures how much the growth of Dx deviates from parabolic. For a ≈ 1, as is typical for crystals found in ice cores at domes and divides, the deviation is small because then g (a) ≈ 1–2(a – 1)2/3 ≈ 1 by Taylor series expansion. Our simulations below for the five cores would not be much different if we made this approximation, but for completeness we use Eqn (7).
Deformation
Vertical compression distorts most crystals by flattening them. Given the axial symmetry, we account for this flattening by modifying Eqn (3) and Eqn (6) to
and
where the new terms describe a purely geometrical effect, not any strain-induced effects on grain growth (which could influence the terms in K, but which we said we are not modelling). The assumption here, used also by Reference Di Prinzio and NaselloDi Prinzio and Nasello (2011) and made in two earlier studies of crystal size by Reference Madsen and ThorsteinssonMadsen and Thorsteinsson (2001) and Reference SvenssonSvensson and others (2003), is that the strain rate felt by the crystals on average – when considering a large random sample of them – equals the bulk strain rate, despite the different orientation of their slip systems and the presence of grain boundaries.
Polygonization and the role of dislocations
Deformation can cause dislocations to pile up at internal tilt walls (subgrain boundaries) in a crystal, dividing it into subgrains that are initially misoriented from each other at low angles. Polygonization occurs when further rotation of the subgrains causes a subgrain boundary to increase misorientation angle and form a grain boundary (Reference PoirierPoirier, 1985; Reference AlleyAlley, 1992). Hence our model of the polygonization process involves dislocation dynamics.
Reference De La Chapelle, Castelnau, Lipenkov and DuvalDe La Chapelle and others (1998) proposed that in polycrystalline ice with a mean grain size D, the evolution of the mean density of (immobile) dislocations ρ can be described by
where b is the Burgers vector (4.5 × 10–10 m) and α 0 is a factor of order 1. The idea is that dislocations created during deformation and acting as carriers of strain become immobile after travelling a mean distance ~ D and add to the dislocation density (a work-hardening effect). At the same time, dislocations are removed as migrating grain boundaries sweep past them (a recovery process). The two terms on the right-hand side of Eqn (10) model the production and recovery rates respectively, and predict both rates to decrease with the grain size. The first term utilizes the Orowan relation (Reference PoirierPoirier, 1985), and α 0 in the second term can be raised to mimic higher dislocation density near grain boundaries. Note that the definition of ρ as a mean dislocation density does not imply that the actual dislocation density is homogeneous across or within crystals (see our discussion later of migration recrystallization). Also, Eqn (10) is a highly simplified model. To explain the behaviour of polycrystalline materials, metallurgists have built sophisticated models of dislocation storage, which, for instance, distinguish different kinds of dislocations and relate their motion to microstructural changes, or include recovery processes such as dislocation climb or annihilation (e.g. Reference BergströmBergström, 1970; Reference NesNes, 1997; Reference Kocks and MeckingKocks and Mecking, 2003).
Reference Montagnat and DuvalMontagnat and Duval (2000) added a second recovery term to Eqn (10) to represent the loss of dislocations to the polygonization process. As explained later, this enabled them to calculate the mean dislocation density ρ along sections of ice cores where steady grain size is maintained by competition between grain growth and polygonization. By using measured profiles of D along several deep cores (including GRIP and Byrd) and estimates for at these sites as inputs to their modified equation, they showed that predicted values of ρ are consistent with values measured by X-ray diffraction on core samples (~ 1011 m–2). Thus they concluded that grain boundary migration and polygonization are fundamental mechanisms underlying the creep rheology of polar ice because these processes accommodate dislocation slip and counteract strain hardening. We then we have use a key idea behind their calculation to develop our model below.
Both Eqn (10) and its modified form by Reference Montagnat and DuvalMontagnat and Duval (2000) involve crystal size, but neither of them describes how crystal size evolves. In order to link the changing dislocation density ρ to the rates of change of Dz and Dx , the crux for us is to link it first to the rate of crystal splitting by polygonization. Weikusat and others (2011), Reference Hamann, Weikusat, Azuma and KipfstuhlHamann and others (2007) and Reference KipfstuhlKipfstuhl and others (2006) have used microscopy and diffraction techniques to study the properties of subgrain boundaries in deformed poly-crystalline ice at high resolution. Their observations reveal the immense complexity of the sub-crystal processes leading to polygonization. Here, in order to form a simple plausible model of the polygonization rate without resolving these processes, we conjecture that dislocation removal by polygonization occurs at a rate proportional to ρ when averaged over many crystals; i.e. the denser the dislocations, the more of them (per unit time) participate in forming new grain boundaries. Therefore we suppose
in which D is now expressed in terms of directional grain sizes, and P, the polygonization rate factor, encompasses influences such as temperature, microstructural properties and impurity effects. From a statistical thermodynamical viewpoint, one may expect an Arrhenius temperature dependence, P ∝ exp(–Qp /RT), with Qp being the activation energy for limiting processes (e.g. lattice diffusion for subgrain boundaries to reorganize into grain boundaries). In Eqn (11) we have multiplied the dislocation production rate term of Reference De La Chapelle, Castelnau, Lipenkov and DuvalDe La Chapelle and others (1998) and Reference Montagnat and DuvalMontagnat and Duval (2000) by a constant 1/ β (>1), because most dislocations travel a distance limited by a chord of the crystal cross section and not the full grain diameter D. We use β = π /4 (assuming an elliptical section) in our modelling, but acknowledge that β could be smaller if most dislocations initiate within the crystals and become immobile before reaching grain boundaries.
How fast is new grain boundary area created for a given rate of consumption of dislocation density by polygonization? Here we use Reference Montagnat and DuvalMontagnat and Duval’s (2000) idea. At a subgrain boundary with misorientation angle , the dislocation spacing h is given by
(Reference PoirierPoirier, 1985). Reference Montagnat and DuvalMontagnat and Duval (2000) assumed that subgrain boundaries transform into grain boundaries when θ reaches a small critical angle θ
c = 5°, so that h ≈ b /θc. Hence, if d ρ
– of dislocation length per unit volume is involved in polygonization, then the new area of grain boundary per volume is 
(cf. Eqn (10) of Reference Montagnat and DuvalMontagnat and Duval, 2000). Our polygonization term in Eqn (11) implies d ρ
– = Pρ dt, and consequently
We next consider the impact of this area increase on the directional grain sizes. For an overall grain size D, the specific grain boundary area (area per unit volume; m–1) is
where C ≈ is a dimensionless constant related to grain shape. In our model, Sv may be conceived as the sum of the area of essentially horizontal grain-boundary surfaces Sv ,H and the area of essentially vertical grain-boundary surfaces Sv, V, with
and c 1 ≈ 2 and c 2 ≈ 1. Reference Durand, Persson, Samyn and SvenssonDurand and others (2008) studied the distribution of misorientation angles between neighbouring grains in the upper part of the NGRIP core and the distribution of their c –axis orientations, and deduced that polygonization occurs at all depths and is isotropic in that core despite bulk deformation. Therefore we apportion the area generation rate in Eqn (13) isotropically, to write
with the constant f equal to 1/3.
Coupled system
Gathering Eqns (11), (16a) and (16b), with the last two rearranged to incorporate the effects of grain growth and deformation from Eqns (8) and (9), we have
where g is given by Eqn (7). This is our complete model. If we had formulated it for the overall grain size D and ignored deformation, then Eqns (18) and (19) would become
Our model differs from Reference Montagnat and DuvalMontagnat and Duval’s (2000) model in distinct ways. Most notably, those authors treated the grain size as known; the assumption of steady state in Eqn (20) then allowed them to determine the dislocation loss rate P ρ = K θ c =bD 3 in their model), which is used in Eqn (17) (without the –correction) to find the equilibrium dislocation density EQ for the region of steady grain size (D EQ) in an ice core:
(summary by Reference Schulson and DuvalSchulson and Duval, 2009, p. 128). In contrast, our model in Eqns (17–19) predicts co-evolving grain size and dislocation density and the depth profiles of both of these quantities. By modelling grain size directionally, we also account for textural anisotropy caused by the deformation.
A further difference is that while Reference Montagnat and DuvalMontagnat and Duval (2000) supposed (following De La Reference De La Chapelle, Castelnau, Lipenkov and DuvalChapelle and others, 1998) that polygonization becomes active only below a certain depth in the ice column (400 m at Byrd and 650 m at GRIP), polygonization is modelled here as occurring simultaneously with grain growth and deformation at all depths. The notion that, within a continuum-modelling context, polygonization does not switch on abruptly at depth but can occur at a low rate in the shallow subsurface has long been held by one of us (T.H.J.), on the ground that the continuum process refers to splitting of crystal collections in a statistical sense rather than of individual crystals. A continuous rate description for polygonization is consistent with the findings of Reference Durand, Persson, Samyn and SvenssonDurand and others (2008) that polygonization operates at shallow depths at NGRIP, and is used in the models of Reference MathiesenMathiesen and others (2004) and Reference RoessigerRoessiger and others (2011) and in earlier considerations by Reference Jacka and LiJacka and Li (1994) and Reference Placidi, Faria and HutterPlacidi and others (2004).
In studying the crystal-size profile of the EPICA Dome C core, Reference DurandDurand and others (2006a) also devised a model that includes coupling between grain-size and dislocation-density evolution. While these authors considered impurity effects in detail, they ignored deformation and anisotropic grain growth and treated the loss of dislocations to polygonization differently. Rather than express this loss via a rate factor (P), they posited a critical value of dislocation density ρ c = 2θc /bD, based on the assumption that the subgrain size equals the grain size, below which polygonization is turned off (equivalent to setting Pρ = 0 in Eqn (17)); otherwise ρ is locked at the critical value (implying a dislocation loss rate Pρ that conspires with the other terms in Eqn (17) to let ρ track ρ c). In contrast, our model assumes no such discontinuity and calculates the dislocation loss rate at all depths. We will see below that a constant P –value suffices for it to match measured crystalsize data.
Three other observables in an ice core may be predicted by our model: the mean horizontal and vertical crystal areas 
and 
respectively and the aspect ratio a (=Dx
/Dz
). Differentiating these variables with respect to t and using Eqns (18) and (19) yields the evolution equations
and
if c 1 ≡ 2 and c 2 ≡ 1, the expressions in the two pairs of square brackets reduce to (a + 1)/3 and (a – 1)/3 respectively (we assume f = 1/3).
Migration Recrystallization
Although we have chosen not to model migration recrystallization, it is worth considering how Eqns (17–19) might be extended to treat this process, and the obstacles. During dislocation creep, a difference in the dislocation density between two neighbouring crystals (associated with their stored strain energies) causes a driving force on the motion of their common grain boundary – additional to the driving force due to grain-boundary energy – that promotes grain-boundary migration toward the crystal with higher strain energy (Reference PoirierPoirier, 1985). For a deforming crystal population, this mechanism can change the motion of grain boundaries by accelerating some of them, and also nucleate crystals. Reference Weikusat, Miyamoto, Faria, Kipfstuhl, Azuma and HondohWeikusat and others (2011, p. 469) estimated that, assuming a dislocation density difference of 1012 m–2 (an order of magnitude larger than our simulated mean dislocation densities below), the strain-induced energy exceeds the grain-boundary energy for a nominal value of the grain-boundary curvature observed in the EDML ice core.
In our model, the obvious terms that could be changed to incorporate this mechanism are those three containing K, because K encapsulates the driving force (e.g. they could be made dependent on dislocation density). But the mean density ρ does not suffice for this consideration: we need to know the variation of dislocation density across crystals (let us call this 
) and preferably within crystals too. Making arbitrary estimates of (e.g. from the mean) does not help; its underlying processes must be modelled. One way to do this is to modify Eqn (17) to track the dislocation density on a per crystal basis, by removing the averaging on ρ, Dz
and Dx
, so that the dislocation production rate (hence dislocation density) in differently sized crystals will differ. Thus we reformulate the model for a population of crystals to follow their distributions of ρ, Dz
and Dx
as probability density functions and follow the statistics of neighbour–neighbour relationships. This approach demands many local considerations about the geometrical dependences of processes, which is why recent efforts have tackled the problem with numerical simulations (utilizing fast Fourier transform, cellular automata or finite-element techniques) that resolve the grain-boundary network (Reference MontagnatMontagnat and others, 2014). We think that a continuum population approach is feasible and can yield valuable insights that complement the simulations, because an example exists in the work of Reference MathiesenMathiesen and others (2004).
Another point we wish to make is that because strain-induced grain boundary motion enlarges some crystals and shrinks others and participates in crystal nucleation, whether it increases or decreases the mean crystal size of a sample, and how fast relative to other processes, cannot be judged from its visible occurrence alone without proper calculations (such as those mentioned above). As far as we know, such calculations have not been made for the cores in Figure 1, nor most other cores, and are much needed now. Research on the importance of migration recrystallization should also differentiate between its impacts on mean crystal size, ice fabric and rheology, and be careful about generalizing the relative importances of processes from one core to another. What we can infer safely from this discussion is that dislocation density and crystal size are pervasively coupled in the crystal-size problem.
Preliminary Analysis of the Model
The model will be applied to the upper part of the ice column where is near-constant, and the ice near-iso-thermal, so that the rates K and P are assumed constant. As initial conditions (at t = 0 and on the ice-sheet surface), we impose small Dz, Dx
(≡ Dz
) and ρ to represent fine-grained ice with equiaxed grains that have undergone little deformation. We ignore complex processes in the firn and use the same model there.
With constant , K and P, the model’s behaviour can be anticipated before we evaluate it with data. From the interactions described by Eqns (17–19), we expect the dislocation density ρ, grain lengths Dx
and Dz
and areas A
H and A
V to increase initially (at shallow depths) and equilibrate at large t (at depth). In this model, near-surface grains grow primarily by (classical) grain growth because they contain few dislocations to drive fast polygonization; but their small size means that deformation creates dislocations at a high rate that outstrips dislocation removal by polygonization, so that ρ increases. As we descend into the ice column (age increases), polygonization speeds up, with two consequences. The removal rate of dislocations increases to rival their production, which is slowed by a larger grain size. Also, polygonization drives grain-size reduction at a rate that rises to offset the rate of grain growth (the former rate increases, and the latter rate decreases, with grain size). The outcome is coupled evolution of grain size and dislocation density toward steady values, with the grains flattening under vertical compression.
This behaviour may be deduced from the equation for each variable, where a growth term dominates when the variable is small, and a large negative term counteracts growth when the variable is large. Notice while polygonization provides such stabilization in each equation, there is additional stabilization in Eqn (17) from the term representing grain boundary migration, and in Eqns (18) and (23) from the thinning (-) term. Setting d/dt = 0 in Eqns (17–19) yields simultaneous equations for the equilibrium values of ρ, Dz
and Dx
, but since their coupling precludes a closed-form solution, we solve them numerically below. In the Appendix, we analyze the isotropic version of the model mathematically to show that its solution stabilizes towards equilibrium but may oscillate for some model parameters, with ρ or D overshooting the equilibrium before tending to it. Slight overshooting in ρ is observed in two of our simulations (Fig. 2b and c).
Modelled depth profiles of crystal size for the five ice cores in our study, down to the bottom end of their steady grain-size regions, as defined in Figure 1. The rows refer to the (a, b) GRIP, (c) NGRIP, (d) GISP2, (e) Byrd and (f) Law Dome DSS ice cores. Also shown are crystalsize measurements from Figure 1; except for the GISP2 crystal-width data of Reference GowGow and others (1997) (d, second panel), these measurements have been multiplied by 1.5 to correct for the sectioning effect. In each panel, the curves depict simulated profiles of mean crystal width Dx
, height Dz
, horizontal area A
H, vertical area A
V and aspect ratio a, and crystal dislocation density ρ. These profiles have been optimized to fit different sets of measurements, as explained by the keys at the left-hand end. In (a), the points for A
H derive from measurements of Dx
via . In (a-c) the points for a derive from the measurements via Dx
/Dz
.
Equation (24) shows that grain flattening is driven exclusively by vertical compression (second term on the right-hand side) and occurs via positive feedback, but that two other processes oppose this by causing elongated grains to become equiaxed over time. Anisotropic grain growth alone (the first term on the right of the equation) implies
whereas polygonization alone (the last term) implies
(if we take c 1 ≡ 2 and c 2 ≡ 1); both rates here are positive when a < 1 and negative when a > 1, so these processes drive the aspect ratio towards 1. That grain-boundary migration driven grain growth has this effect is evidenced by the simulations of Reference Di Prinzio and NaselloDi Prinzio and Nasello (2011) and a cellular automaton model of Reference Raghavan and SahayRaghavan and Sahay (2009). That polygonization can have a similar effect is intuitive when we consider how it splits flattened grains into daughter grains whose shapes will be rounder on average.
The model highlights contrasting effects of the vertical compression on the mean horizontal and vertical crystal areas. The signs of the -terms in Eqns (22) and (23) indicate that compression causes radial spreading of crystals to increase A
H, but in any vertical cut causes material loss into the normal direction and net shrinkage of A
V. This has implications for crystal areas measured from thin sections at shallow depths. Assuming small, equiaxed grains (small A
H ≈ A
V ≈ A, a ≈ g ≈ 1) and negligible polygonization rate there, the growth rates in Eqns (22) and (23) are approximately
These expressions show that apparent grain-growth rates inferred from crystal sizes measured on horizontal and vertical thin sections will diverge by 1: 5
A, and be enhanced and reduced (respectively) from the non-deforming growth rate; thus estimates of K (and K
0 and Q) derived from such measurements may have errors that require strain-rate correction, especially at high-accumulation sites where is large. The errors are negligible (a few per cent of K) at the GRIP, NGRIP, GISP2 and Byrd core sites and slightly larger at Law Dome. The considerations here may be repeated for pure shear beneath linear ice divides. In this case, vertical compression enhances dA
H/dt, but its effect on dA
V/dt depends on the azimuth of the vertical cut (this rate is unaffected if the cut is perpendicular to the divide axis and reduced if the cut parallels the axis). Even though our model does not account for the effects of impurities and migration recrystallization on grain growth, these results underline the importance of determining crystal size in all three principal directions from ice-core measurements.
Testing the Model with Ice-Core Data
We now examine how well the model can reproduce the measured crystal-size depth profiles in Figure 1; in so doing, we provide benchmark simulations for the five ice cores. Table 1 summarizes key information about these cores, including the crystal-size measurements made on them and the methods used, and the mean and standard deviation of crystal sizes in the steady region on each profile. A further crystal-size profile from the ice core at Siple Dome, West Antarctica, could be used, but we decide not to do so, given complications at that site: the profile itself does not show a clear steady region (Reference Di Prinzio, Wilen, Alley, Fitzpatrick, Spencer and GowDiPrinzio and others, 2005), and the temperature and strain rate in the ice, which vary strongly with depth today (Reference Gow, Engelhardt and HondohGow and Engelhardt, 2000; Reference ElsbergElsberg and others, 2004), might have changed during the Holocene due to increasing accumulation rates and thinning of the dome (Reference WaddingtonWaddington and others, 2005).
Optimization procedure
We fit the model to the measurements in each dataset in two steps. Since a prerequisite for a good match is that the model simulates the right steady crystal size, we first optimize this match, by integrating Eqns (17–19) numerically to find the equilibrium values of Dx and Dz (and hence of A H and A V) and by tuning the polygonization rate factor P. Generally, the smaller the steady crystal size measured, the higher is the P- value needed, if other factors remain constant. Age-to-depth conversion is not required in this step, nor is the outcome dependent on initial conditions. We then resimulate the model with P at its optimal value, this time adjusting the initial crystal sizes Dx (t = 0) and Dz (t = 0) (which are kept equal) until the solutions Dx (t) and Dz (t) capture the observed grain growth profile at shallow depths and its curved approach to equilibrium as well as possible, when these solutions are converted to depth solutions using the age–depth scale of the core (Table 1). In this second step, the initial mean dislocation density ρ (t = 0) is set at the low value of 1010 m–2, which is typical of annealed ice (Reference Montagnat and DuvalMontagnat and Duval, 2000). This step puts the model to a different test because P is fixed and we do not alter other parameters.
In both steps, the model is solved as a forward problem down to the bottom of the steady grain-size region. We do not reconstruct the histories of ice temperature, accumulation rate and age–depth scale of the cores (e.g. Reference MorlandMorland, 2009). The model parameters are taken to be θ c = 5° (following Reference Montagnat and DuvalMontagnat and Duval, 2000), α 0 = 1, β =π /4, f = 1/3, c 1 = 2 and c 2 = 1, and for each core we apply the strain rates and temperatures listed in Table 1. The surface temperature is used at GRIP, NGRIP and GISP2 because the ice at these sites is near-isothermal through the depth interval of the simulations. The ice temperature varies more with depth at Byrd and Law Dome, so for these sites we use the mean temperature of the interval (in brackets in Table 1). Note that the assumption of unaxial compression is more approximate for NGRIP, GISP2 and Byrd, as these sites do not lie at dome summits.
Because the data from the five cores describe different crystal-size variables, and sometimes multiple measurements of a variable have been made on the same core, we can choose different combinations of variables for the model to match. When investigating each core, we optimize the model once if data exist for only one variable (e.g. A H for Byrd); we make separate optimizations if both crystal length and area datasets or multiple datasets on the same variable exist (e.g. the crystal width/height data for GRIP of Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others (1997) and of Reference Svensson, Durand, Mathiesen, Persson, Dahl-Jensen and HondohSvensson and others (2009)). Moreover, if data exist for both measures of length (Dx and Dz ) or both measures of area (A H and A V), we match the model to their steady values in a single optimization by minimizing the mean square error of fit to both measures. This pattern of application leads to three model optimizations for the GRIP core, two for the NGRIP and GISP2 cores and one for the Byrd and Law Dome cores, and thus nine sets of model-simulated depth profiles and nine estimates of the polygonization rate factor P.
Before these optimizations, we multiply most of the crystal-size data by a factor to correct for the ‘sectioning effect’: thin-section samples of a core rarely intersect crystals where they are widest, so the measured mean of crystal length or area underestimates the true mean. In this regard, different methods of measurement also suffer from different biases (e.g. Reference JackaJacka, 1984; Reference Durand, Gagliardini, Thorsteinsson, Svensson, Kipfstuhl and Dahl-JensenDurand and others, 2006b). We correct all of the records in Figure 1 except the crystal-width record of Reference GowGow and others (1997) for GISP2 (Fig. 1d, crosses), because none of those records have been treated for the sectioning effect, and because Dx, Dz, A H and A V in our model refer to true mean crystal dimensions. For simplicity we use the same factor of 1.5 in all corrections. This value derives from stereological calculations assuming monosized spheres and holds for both crystal length and area (Reference UnderwoodUnder-wood, 1970, table 4.1). We do not vary its value to account for the effect of changing aspect ratio of the crystals, because observations of their true 3-D shape are lacking and the simulated values of the aspect ratio a below are close to 1 (Fig. 2). Nor do we reconcile the different methods of measurement in this paper, as doing this properly requires reapplying them to the same (and many) thin sections and comparing the results.
We leave the GISP2 crystal-width data of Reference GowGow and others (1997) unadjusted because these authors had effectively corrected them for the sectioning effect, by measuring the 50 largest crystals in each thin section (footnote to Table 1); multiplying the data by 1.5 would overcorrect them. Although Reference Durand, Gagliardini, Thorsteinsson, Svensson, Kipfstuhl and Dahl-JensenDurand and others (2006b) pointed out that Gow’s correction method is not always consistent as it samples one end of the grain-size distribution, Figure 1d suggests that it works reasonably well at GISP2: throughout the core, the mean crystal widths measured by Reference GowGow and others (1997) are indeed roughly twice those measured by Reference Alley and WoodsAlley and Woods (1996) with another method (which need the correction factor).
To estimate K 0 and Q in Eqn (2), we use the grain growth rates and temperatures compiled by Reference Cuffey and PatersonCuffey and Paterson (2010, table 3.1) for shallow firn and ice at 13 different polar sites. To reiterate, the growth law with these ‘classical’ parameters cannot resolve the effects of strain-induced grain boundary migration and impurities on the mean grain size, and the classical parameters may already be influenced by these effects. Fitting Eqn (2) to these data yields K 0 = 8.78×106 mm2 a–1 and Q = 42.4 kJ mol–1. However, we correct the value of K 0 here too, for two reasons. First, the growth rates compiled by Cuffey and Paterson refer not to K in our model, but to dA /dt, which is approximately equal to π K /4 at shallow depths (Eqns (27) and (28)); thus a first correction factor is 4/ π. Second, most of the compiled growth rates are based on measured grain sizes that have apparently not been corrected for the sectioning effect; this implies a further correction factor of 1.5. (More precisely, two of the growth rates have been corrected, while a few others may derive from crystal length squared (dD 2/dt) rather than area and need a larger correction factor of 1.52 for the sectioning effect, but the level of information in the sources used by Cuffey and Paterson makes this uncertain.) Here we disregard the subtle non-uniformities in these data and adopt 1.5 × 4/ π as a fixed combined correction factor, and hence the corrected rate constant K 0 = 1.68 × 107 mm2 a–1, for our modelling.
Modelled depth profiles
Figure 2 plots the optimal depth profiles simulated by our model alongside the crystal-size measurements, which are shown with the same symbols as in Figure 1 and have been treated for the sectioning effect, as described above. Since the model is meant to capture trends in the measurements, not their rapid variations, we do not report r –square values of fit, which would reflect these variations. In terms of organization, Figure 2a and b show the results of model optimizations made to the GRIP data of Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others (1997) and Reference Svensson, Durand, Mathiesen, Persson, Dahl-Jensen and HondohSvensson and others (2009) separately, while Figure 2c–f show the results for NGRIP, GISP2, Byrd and Law Dome. In Figure 2b and c, the results of different matches made to crystal length and to crystal area are distinguished by curve thickness, as indicated by the key on the left. In Figure 2, we generally show crystal width and height in the first column, area in the second, aspect ratio in the third, and dislocation density in the fourth. Exceptions appear in Figure 2d (for GISP2), where the first two panels are both devoted to crystal length, and in Figure 2f (for Law Dome) where they are devoted to area.
Figure 2 shows that the model can simulate the measurements’ overall trajectories successfully, although less so in the optimization for Law Dome. Putting aside this case and considering the other eight optimizations first, the simulated curves of crystal size and dislocation density show the expected initial growth and equilibration. Importantly, when the model is optimized to match the mean steady grain size(s) (be this Dx, Dz, A H, A V or their combinations), it also fits the initial growth trend of the measurements and their curved approach to steadiness, including the depth scale of this approach; this match is achieved while P is kept constant with depth in each simulation. (In Figure 2a we verify the excellent match for GRIP by plotting A H against age instead of depth.) The simulated aspect ratios show that crystals flatten progressively with depth, with final values of a in the range 1.1–1.5. Moreover, for the GRIP and NGRIP cores (Fig. 2a–c), where data for both Dx and Dz exist, the model is able to match their mean steady values simultaneously, although the simulated aspect ratios underestimate the measured ratios slightly. The fluctuations or ‘noise’ on most of the measured profiles, not reproducible by the model, suggest physical controls on the grain size unaccounted by us, such as impurities.
Figure 2b and c show that for the GRIP and NGRIP datasets of Reference SvenssonSvensson and others (2003, Reference Svensson, Durand, Mathiesen, Persson, Dahl-Jensen and Hondoh2009), the model optimized for Dx and Dz overestimates A V, and accordingly the model optimized for A V underestimates Dx and Dz . This discrepancy looks uniform and concerns the magnitude, not form, of the curves. We trace its cause to the use of the same correction factor on measured crystal area and lengths. Svensson’s data of Dx, Dz and A V satisfy the relationship A V = DxDz /4 approximately (to within 15%), but after we multiply all of them by 1.5, A V becomes 26–27% smaller on average than what DxDz /4 predicts; this is why the optimizations for Dx and Dz together and for A V yield different results. We have not tried using different correction factors for crystal area and length to remove the discrepancy, given the lack of information (other than a better fit) to justify their choice. Instead, with Svensson’s datasets, we treat the two optimizations as separate and report both of their polygonization rate factors below.
The simulated dislocation densities reach equilibrium values of the order of 1011 m–2 and vary from core to core. As expected, increasing α 0 to 2 or 3 in the model lowers these densities because grain boundary migration removes dislocations faster (this was fully analysed by Reference Montagnat and DuvalMontagnat and Duval, 2000). Then the optimal polygonization rate factor P changes correspondingly, but the model’s ability to match the observed profiles and our evaluations later about the spread of P are unaffected.
We turn to the model optimization for the DSS core at Law Dome in Figure 2f. For this core, Reference Li, Jacka and MorganLi and others (1998) measured A H and A V by counting all grains within a known area on each thin section (Table 1), so large errors or biases in their data are unlikely. However, their measurements exhibit a strong scatter that is visible even in logarithmic scale; it is therefore difficult to judge how plausible the fit between model and data is. It may in fact be questioned whether the measured crystal areas in the ‘steady region’ identified by us (400–700 m depth; Fig. 1f) represent a steady state. In this regard, Figure 2f shows that the simulated profiles of crystal size and dislocation density have not reached equilibrium in this region. If one interprets the measured crystal areas there to be still increasing with depth, then the model would be capturing roughly the right trend. Even so, Figure 2f shows another problem. The rapid growth of the shallow measurements (at <200 m depth) cannot be modelled for any prescribed initial grain size. The model underestimates the grain-growth rate in this region, which according to Reference Li and JackaLi and Jacka (1999) is ≈3.84 × 10–2 mm2 a–1 for A H, when uncorrected for the sectioning effect. This rate is enhanced by about three times compared to the rate predicted by the classical growth law (best-fit Eqn (2); again ignoring correction), consistent with what Reference Li and JackaLi and others (1998) and Reference Li and JackaLi and Jacka (1999) have reported before. The geometric effects of vertical compression that we discussed in relation to Eqns (27) and (28) are limited at Law Dome and cannot explain the observed large enhancement in the growth rate of A H; besides, an enhanced growth rate is observed for A V as well. Reference Li, Jacka and MorganLi and others (1998) suggested that these enhancements are due to dynamic recrystallization processes resulting from high compressive strain rate at the site. A likely candidate is migration recrystallization; as discussed before, future work could explore how to model its physics.
Polygonization rate factor P
Table 2 lists the nine values of P found from the optimizations. Figure 3 plots them against the ice temperature T and also 1/T. The error bars in P derive from model fit to the upper- and lower-bound steady grain sizes in Table 1. As some of the uncertainty with the model optimization for the Law Dome DSS core is accounted by the error bar, we include its P –value tentatively in this discussion.
Polygonization rate factors P (a–1) retrieved by optimizing our model to match the steady crystal sizes listed in Table 1 for the five ice cores. The indices (i)–(iv) follow those in Table 1. For all but the crystal-width data at GISP2 measured by Reference GowGow and others (1997), the steady crystal sizes are multiplied by 1.5 to correct for the sectioning effect prior to matching
Polygonization rate factor P in the five ice cores of our study (see labels under the plot) against ice temperature T and its reciprocal variable 1/T. Each estimate of P is found by matching the model in Eqns (17–19) to the mean crystal size in the steady crystalsize region of the corresponding core (corrected for the sectioning effect in most cases, as explained in the text). Plot symbols indicate the crystal-size variables being matched in each optimization. The indices (i)–(iv) follow those in Tables 1 and 2.
The results in Figure 3 highlight the interesting question of what governs the polygonization rate in polar ice. As mentioned before, P might increase with temperature following an Arrhenius dependence, P ∝ exp(–Qp
/RT), where Qp
is activation energy. The P –values found for the different cores and datasets vary considerably in the range 10–5 to 10–2 a–1, although a constant (depth-invariant) value of P has enabled the model to match the measured crystalsize profiles in each simulation. If we regard the P –values for GRIP, NGRIP and GISP2 as a single cluster for central Greenland (these sites have similar temperatures), and the Byrd and Law Dome P –values as two other clusters, then their general trend across the plot suggests that P decreases with temperature: this is opposite to the expected Arrhenius dependence. The strain rate is highest at Law Dome, second highest at GISP2, and lowest at Byrd and NGRIP (Table 1), so we do not see a clear pattern of control by And if we increase the model parameter α
0 or θ
c in the optimizations, then all P –values in Figure 3 will shift up, but their pattern remains similar. This means that unknown factors must control the polygonization rate factor P, whether or not temperature plays a role. One possibility is the evolving microstructure within crystals, which might be accountable partially through the dislocation density; another is the effect of impurities on subgrain boundary formation. Other recovery mechanisms beside polygonization may also be significant in our cores, and they are not distinguished by the model and are lumped into the rate P.
Conclusions
We have explored a continuum model of crystal-size evolution that combines several ingredients: crystal flattening under compression, directional grain growth, polygonization, and the build-up and removal of crystal dislocations. The model is not a complete description because it does not resolve impurity and migration-recrystallization effects. Simulations show its ability to match the upper part of the crystal-size profiles from the GRIP, NGRIP, GISP2 and Byrd ice cores (including both crystal width and height data from the first two cores). Given recent observational evidence from other ice cores of the importance of migration recrystallization, our study might seem controversial for revisiting the competition between grain growth and polygonization – an idea from the ‘tripartite’ hypothesis. However, the competition is studied here as a subset of the full interactions and thus a building block towards the general theory of crystal-size evolution and dynamic recrystallization in ice. For those engrossed in debates surrounding the validity of the tripartite hypothesis, our finding that the simulations did not universally fail to fit the profiles means that two mutually exclusive possibilities stand:
The impacts of migration recrystallization and impurities, or of their variations with depth, on the rate of change of mean crystal size are limited in comparison with the impacts of the modelled processes, within the depth ranges of the cores concerned.
These impacts significantly control the rate of change of mean crystal size on the observed profiles, but the nature of these controls is such that the profiles are reproducible by a model formulated on incomplete physics.
Consequently, more extensive calculations are required to quantify and disentangle the depth-varying impacts of normal grain growth, impurities, migration recrystallization and polygonization on the crystal size in the cores studied, and generally in all existing cores. The impacts of these processes on ice fabric and rheology need to be quantified separately, and we caution against qualitative generalizations about the relative dominance of processes from one core to another.
In their review, Reference Faria, Weikusat and AzumaFaria and others (2014) proposed a conceptual model for the evolution of the mean (isotropic) crystal size D by the equation
where F is a rate function governed by the whole suite of recrystallization processes (and is hence difficult to determine at present), and the form of F behaves as an ‘attractor’ in the 3-D state space of D, and T, so that F = 0 defines the steady-state grain size. (We have written Eqn (29) in a mathematical notation different from the one used by Reference Faria, Weikusat and AzumaFaria and others (2014), after consulting the first author.) In the current paper, we have explored some of the processes determining F, and our results support a dynamical-system view of the problem. Moreover, we now know that Eqn (29) is not an adequate description of crystal-size evolution, because the mean crystal size and mean dislocation density (ρ) evolve in a strongly coupled manner. For polar ice in axially symmetric deformation, three coupled first-order differential equations are necessary to track the evolution of two directional crystal sizes and dislocation density (e.g. Eqns (17–19) or (22–24)); this introduces two more state variables into the system besides D. Under non-axially-symmetric configurations, mean crystal sizes in three orthogonal directions need to be tracked. And under axially symmetric configurations where crystal flattening is small (Dz
≈ Dx
≈ D), the two grain-size equations reduce approximately to a single evolution equation for D. However, the system of co-evolving and D cannot be further reduced (Fig. 2; Appendix). Finally, we expect the inclusion of migration recrystallization to complicate the coupling (e.g. by adding a measure of the spatial variation of dislocation density or of stored strain energy as a state variable).
One of our motivations is to evaluate data from multiple ice cores in a comparative study. Most previous simulations of crystal size investigated its variation in individual cores, although the model of Reference MathiesenMathiesen and others (2004) has been fitted to the GRIP and NGRIP profiles in separate publications (Reference MathiesenMathiesen and others, 2004; Reference Svensson, Durand, Mathiesen, Persson, Dahl-Jensen and HondohSvensson and others, 2009). With new data from future cores, more comprehensive comparisons will be possible.
It was not straightforward for us to work with different ice cores due to varying practices of measuring and reporting crystal sizes. This calls for a much-needed standardization of the measurement methods. Existing core samples could be remeasured with different methods too, to increase the compatibility between datasets.
An outstanding question is what governs the polygonization rate. Whereas some authors described this rate in terms of the average time crystals take to subdivide (e.g. De La Reference De La Chapelle, Castelnau, Lipenkov and DuvalChapelle and others, 1998), we have quantified it by the rate of loss of dislocation density during polygonization, normalized by the instantaneous dislocation density: the resulting ratio is our rate factor P. A constant P enables our model to match individual crystal-size profiles – not just their steady values, but also their growth and approach to steadiness. And as the temperature T and strain rate are constant in each simulation, the triplet (T, , P) completely determines the modelled profiles in the depth range of interest. But while P should depend on temperature, the pattern of P –values found for the different cores (Fig. 3) suggests that additional factors operate. Microstructural analysis, which has been used to map various sub-crystal-scale mechanisms at high resolution, may shed light on these factors.
Acknowledgements
We thank the Antarctic Climate and Ecosystems Cooperative Research Centre in Hobart, Australia, for hosting and financially supporting a visit by F.N. to collaborate with T.H.J. in August 2012. Critical reviews by S. Faria (editor), P. Bons and N. Azuma helped us improve the manuscript.
Appendix
Here we examine the local stability of the isotropic grain-size evolution model in Eqns (17) and (20)
to understand how its solution approaches equilibrium. To simplify the algebra, it is useful to define the constants
Setting d/dt to zero in Eqns (A1) and (A2) and solving for D and ρ yields the equilibrium grain size
and the equilibrium dislocation density
Writing D =D EQ + d and _ = _ EQ + r, expanding Eqns (A1) and (A2), and linearizing, leads to the system of equations
Where
Let the solutions ∝ d and r be /exp(σt), and then non-trivial solutions have the eigenvalue _ which satisfies
or (after expanding and doing the algebra)
As this quadratic equation for σ has positive coefficients, its roots are either real and negative, or complex with a negative real part. Thus d and r decay exponentially, implying D and ρ reach equilibrium stably decay is oscillatory if
in which the constant is
For 0.5 < α0 < 4.5, which holds true for the values of _ 0 used in this paper (_ 0 ¼ 1) and experimented by Montagnat and Duval (2000) (α 0 ¼1–3), the condition in Eqn (A10) for oscillations to occur is satisfied if
Outside this range of α0, the condition is never satisfied (for B 3 > 0).
Finally we relate B 3 to the model parameters. By substituting for the equilibrium grain size D EQ from Eqn (A4) and defining
Eqn (A11) becomes
which is a decreasing function that has a 1/B
4-type singularity as B 4 → and tends to 0 as B
4→ ∞. Thus, ultimately the parameter that governs whether the decay oscillates is 
For sufficiently small P, B
4 will be large enough and B
3 small enough to cause oscillations.
