Introduction

Many classical studies of polar ice microstructure report an evolution of the mean grain size with depth according to what can be called the ‘three-stage model’ (Reference Gow and WilliamsonGow and Williamson, 1976; Reference Herron and LangwayHerron and Langway, 1982; Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997): in the upper few hundred metres, grain size increases steadily with depth; below a certain intermediate depth (400–700 m), the grain size stabilizes and remains roughly constant; finally, at great depths (approximately the last 300 m before reaching bedrock, where temperature exceeds −10°C (Reference De La Chapelle, Castelnau, Lipenkov and DuvalDe La Chapelle and others, 1998; Reference Duval and HondohDuval, 2000) the grain size significantly increases again. Here we only deal with the upper two regions, where grain size first increases and then stabilizes.

The initial steady increase in grain size is usually explained by static (‘normal’) grain growth (Reference SmithSmith, 1964; Reference Alley, Perepezko and BentleyAlley and others, 1986; Reference Weaire and RivierWeaire and Rivier, 2009), defined as growth that is only driven by the reduction of free energy of the grain boundaries. The increase in grain size, expressed in mean radius, r, from a starting grain size, r ₀, is usually described by (Reference Anderson, Hansen, Juul Jensen, Leffers and RalphAnderson, 1986; Reference Glazier, Gross and StavansGlazier and others, 1987; Reference Weygand, Bréchet, Lépinoux and GustWeygand and others, 1998)

(1)

The growth exponent n has a theoretical value of 2 in ideal static grain growth of grains with isotropic properties (Reference Glazier, Gross and StavansGlazier and others, 1987). In natural systems, the exponent is usually found to be >2. Any other process or factor that influences grain growth tends to increase n, such as anisotropic boundary energies, pinning, etc. (Reference GowGow, 1969; Reference GowGow and others, 1997; Reference Bons, Jessell, Evans, Barr, Stüwe, Koyi and MancktelowBons and others, 2001; Reference DurandDurand and others, 2006). The parameter k is normally treated as a temperature-dependent material property that is a function of only the boundary energy γ(T) and the grain boundary mobility M(T):

(2)

where T is the temperature and the factor k ₀ is generally assumed to be constant. For ideal static grain growth the value of k ₀ is ∼0.5 in three dimensions and ∼1.12 in two dimensions (Reference MullinsMullins, 1989; Reference Weygand, Bréchet, Lépinoux and GustWeygand and others, 1998). Below we show that in practice k ₀ is actually not a constant, but in fact depends on the microstructure (the ideal case being a particular instance). The factor k ₀ itself is usually difficult to determine from experiments or measurements in nature (i.e. polar ice caps). This is because one normally only obtains k, which also includes the surface energy and grain boundary mobility. If k depends on microstructure through the parameter k ₀, one cannot apply k obtained from one study to another situation where the microstructure may be different. In this paper we show that k ₀ varies with microstructure and how ignoring this may lead to erroneous results if applied to polar ice caps.

If static grain growth were the only process operating in polar ice, the grain size should increase steadily with the age of the ice, and hence with depth. The observation in several ice cores that grain size stabilizes at a certain depth suggests that another process operates which balances the increase in grain size (Reference AlleyAlley, 1992; Reference De La Chapelle, Castelnau, Lipenkov and DuvalDe La Chapelle and others, 1998; Reference DurandDurand and others, 2006). If this other process leads to a reduction of grain size, a balance between grain-size increase and decrease will be reached at some point. The process usually invoked to explain the grain-size reduction process is polygonization or rotational/continuous recrystallization (Reference Urai, Means, Lister, Hobbs and HeardUrai and others, 1986; Reference AlleyAlley, 1992; Reference Alley, Gow and MeeseAlley and others, 1995; Reference Duval and CastelnauDuval and Castelnau, 1995; Reference Faria, Ktitarev and HutterFaria and others, 2002).

Rotational recrystallization is a deformation-driven process. Deformation by dislocation creep introduces dislocations in the crystal lattice, which can accumulate in planar zones or tilt walls that define regions within a grain with small differences in their lattice orientations. The lattice within these regions or subgrains within a grain thus rotate relative to each other. Progressive rotation of the subgrains with ongoing strain eventually leads to the formation of high-angle grain boundaries, and the subgrains they bound become real grains (Reference ReadRead, 1953; Reference Duval, Ashby and AndermanDuval and others, 1983). Rotational recrystallization can be regarded as a process that effectively splits grains into two or more grains (Reference MathiesenMathiesen and others, 2004; Reference Placidi, Faria and HutterPlacidi and others, 2004). Each split increases the number of grains, N, in a volume by 1. The increase in N, and hence decrease in grain size, thus depends on the split rate f per grain:

(3)

The parameter f may depend on many factors, most importantly on strain rate and hence on deviatoric stress (e.g. through Glenn’s flow law; Reference AlleyAlley, 1992). However, in a first approximation it is usually assumed that the strain rate is approximately constant within the upper part of the core where our calculations apply (Reference Lipenkov, Barkov, Duval and PimientaLipenkov and others, 1989; Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997; Reference Montagnat and DuvalMontagnat and Duval, 2000). The split rate of a grain probably also depends on the size and deformation history of that grain. A split rate proportional to grain size was, for example, assumed by Reference MathiesenMathiesen and others (2004) and Reference Placidi, Faria and HutterPlacidi and others (2004), while Reference MorlandMorland (2009) studied the effect of ice flow history. However, the simplest (but not necessarily realistic) assumption is that f is a constant, not depending on grain size or any other factor. This simplification is permissible here, since this paper is mainly concerned with the influence of microstructure on growth rate, and we do not intend to model a particular ice core. For this case, a simple analytical solution exists for the stable grain size. Assuming that the grain growth exponent n is 2 in Equation (1), one derives (see Appendix)

(4)

Here a is a geometrical factor relating the mean grain radius, r, to the number, N, of grains in a volume. For illustration, by applying this equation to the North Greenland Icecore Project (NorthGRIP) ice-core data (Fig. 1), one obtains a growth constant of k ≈ 5.0 × 10⁻³ mm² a⁻¹ and a split rate of f ≈ 1.5 × 10⁻³ a⁻¹ or once every 650 years. These numbers are within the range of those reported in the literature (Reference GowGow, 1969; Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997; Reference SvenssonSvensson and others, 2003; Reference MathiesenMathiesen and others, 2004). The question, however, is whether the values obtained are realistic and meaningful.

Fig. 1. Fit of analytical model (Equation (4)) to the average grain diameter as a function of age as observed in the NorthGRIP ice core (squares; data from fig. 3 in Reference MathiesenMathiesen and others, 2004). Fit parameters are k = 5.0 × 10⁻³ mm² a⁻¹ and a split rate of f = 1.54 × 10⁻³ a⁻¹ or once every 650 years.

Numerical Simulations

We used the numerical modelling platform Elle (Reference Jessell, Bons, Evans, Barr and StüweJessell and others, 2001; Jessell and Bons, 2002; Reference Bons, Koehn and JessellBons and others, 2008) to simulate the process of grain growth and grain splitting. The Elle software was developed to simulate the micro-structural evolution in materials such as rocks. It has been applied to the simulation of a range of processes, such as static grain growth in anisotropic polycrystals or partially molten rocks (Reference Bons, Jessell, Evans, Barr, Stüwe, Koyi and MancktelowBons and others, 2001; Reference Becker, Bons and JessellBecker and others, 2008), dynamic recrystallization (Reference Piazolo, Bons, Jessell, Evans, Passchier, de Meer, Drury, de Bresser and PennockPiazolo and others, 2002, Reference Piazolo, Jessell, Prior and Bons2004) and strain localization (Reference Jessell, Siebert, Bons, Evans and PiazoloJessell and others, 2005). The main distinguishing features are (1) that it uses a two-dimensional (2-D) image of the actual microstructure, and (2) that it uses operator-splitting to allow a range of different processes to operate on, and modify the microstructure. This means that simultaneously operating processes (such as grain growth and grain splitting) are modelled as isolated individual processes that sequentially modify the microstructures in very small increments.

The microstructure is defined by a contiguous set of polygons that are themselves defined by boundary nodes that link straight boundary segments (Fig. 2). The polygons typically represent individual grains. Changes in the microstructure are achieved by (1) changing the properties of polygons or boundary nodes, (2) changing the position of boundary nodes, which implies a change in shape of the polygons, and (3) creating, removing or reordering boundary nodes and segments. A change in shape can be the result of deformation, for which the finite-element code, Basil, is available in Elle (Reference Houseman, Barr, Evans, Bons, Koehn and JessellHouseman and others, 2008). A change in shape can also be the result of the movement of boundaries (grain boundary migration), for example in the case of grain growth.

Fig. 2. Basic structure of the Elle model. The model consists of polygons which represent grains, and these polygons are in turn defined by boundary nodes (a) that are connected by straight boundary segments. Only boundary nodes with two or three neighbours are allowed in the model. The boundary nodes can move (b) and their movement is determined by the curvature of the boundary of the polygon at that point. Grains are split by the introduction of a new straight boundary that links two existing nodes (c).

The movement of grain boundaries is modelled by sequentially selecting each boundary node, and applying a small incremental displacement that depends on the driving force for migration and the intrinsic boundary mobility. In this study we test the validity of Equation (4) by combining a static grain growth routine that moves grain boundaries, and a split routine that divides grains into two daughter grains.

The normal grain growth routine simulates ideal isotropic growth (without grain boundary energy anisotropy). For each time-step, the routine goes through the list of all boundary nodes and calculates the local radius of curvature, r _c, using the node and its immediate grain boundary neighbours. The velocity, v, of the node in the direction of the centre of the curvature is calculated using

(5)

The node is then moved over a distance Δx for a small time increment Δt. This routine results in ideal growth with a linear increase in mean grain area A, implying a growth exponent of n = 2, and k ₀ = 1.22 (Figs 3 and 4a). This would be the growth exponent as expected from theory (Reference Humphreys and HatherlyHumphreys and Hatherly, 1996). However, growth exponents measured in natural ice may deviate from that value due to other processes not taken into account here.

Fig. 3. (a) Growth curves for models of pure static grain growth. For Mγ = 3.2 × 10⁻³, 6.4 × 10⁻³ and 3.2 × 10⁻² mm² a⁻¹, the average grain area increases linearly with time. (b) Plot of k values measured from simulations as a function of the set value of Mγ. The slope of 1.22 is the value of k ₀.

Fig. 4. Results of numerical simulations with Elle. (a) Static grain growth only, for 6500 years and Mγ = 3.2 × 10⁻³ mm² a⁻¹. (b) Simulation with same starting aggregate and settings as for (a), but with splitting at a constant f = 1.54 × 10⁻³ a⁻¹ added, which leads to the establishment of a stable grain size after ∼4000 years, and a different microstructure compared to static grain growth. Size of box is 72 mm × 72 mm.

The effect of rotational recrystallization was implemented by randomly splitting each grain with a probability of f every time-step for each grain. This probability determines the rate of grain-size reduction by splitting. For this, each grain is selected in turn, and a random number generator determines whether the grain will be split. If so, one of its nodes is randomly selected and a new boundary is constructed across the grain, in a random orientation. Each time, the program checks whether the intended split will cause topological problems, such as intersection of the new boundary with another boundary or that a tiny grain has insufficient available nodes to split between. As a result, some splits are cancelled and a set value of f of 1.54 × 10⁻³ a⁻¹ results in an effective split rate of 1.52 × 10⁻³ a⁻¹, meaning that on average 1.3% of attempted splits are cancelled when a steady state has been established.

As expected, a stable grain size is established as a result of the combination of growth and splitting (Figs 4b and 5). For Mγ = 3.2 × 10⁻³ mm² a⁻¹ (k = 3.90 × 10⁻³ mm² a⁻¹) and f = 1.52 × 10⁻³ a⁻¹), the average stable grain diameter is 3 mm². To compare this result with the analytical model, we must rewrite Equation (4) for the 2-D case:

(6)

Fig. 5. Evolution of the average grain diameterwith time. Static grain growth (Mγ = 3.2 × 10⁻³ mm² a⁻¹) results in a linear increase of grain diameter (dotted line) (Fig. 4a). Adding a constant split rate (f = 1.54 × 10⁻³ a⁻¹) for all grains (Fig. 4b) results in the establishment of a stable average grain diameter (dash-dot line). Applying the same settings to an initially large grain microstructure (dashed line) results in the same steady state as for the initially small grain microstructure. For comparison the data from the NorthGRIP core (Fig. 1) have been plotted as well (squares) along with their fit (solid line).

The average stable grain area predicted by the analytical model (Equation (6)) is similar to the value obtained with the Elle simulation, although the stable state is only reached after ∼4000 years in the simulation. To achieve stabilization of the grain size after ∼2000 years, as in the case of the NorthGRIP data, one has to roughly double both k and f. The discrepancy between the analytical model (Equation (4)) and the numerical simulation can be explained by considering the microstructure (Fig. 4). Static grain growth produces a regular foam texture. The frequency distribution of grain diameter has a maximum at about the average grain area (Fig. 6), and the normalized grain-size distribution is time-invariant (for steady-state growth). When a stable grain size is reached due to a balance between grain boundary migration and splitting, the grain size distribution changes significantly, with an increase of the frequency of very small grains, but also an increase in grains much larger than the average.

Fig. 6. Normalized frequency distributions of grain diameter. Solid line is the average of 16 simulations of only static grain growth (Fig. 4a). Dashed line is for eight simulations after a steady state has been reached by the competition of static grain growth and splitting (Fig. 4b).

The change in microstructure changes the growth behaviour. The relatively abundant small grains have a high boundary curvature and quickly disappear. Yet many new small grains constantly appear because in the model every grain has the same chance of being split, independent of its size. The effect of the widening of the grain-size spectrum is an increase in the growth rate that balances the split rate in Equation (4). This can be seen if one stops the splitting when a stable grain size has settled but grain growth is allowed to continue (Fig. 7). The initial growth rate is over three times higher than the stable growth rate that is reached after the mean grain area has about quadrupled. This implies that the factor k ₀ is not a constant, but a function of the microstructure. For the stable foam texture that results from static grain growth only, k ₀ is 1.22. When the microstructure is the result of a competition and random splitting, the effective value of k ₀ increases to 4.2 (an increase by a factor of 3.5).

Fig. 7. (a) Grain growth experiment (Mγ = 3.2 × 10⁻³ mm² a⁻¹) where splitting (f = 1.54 × 10⁻³ a⁻¹) is turned off after 6500 years. The dashed line shows the growth rate of k ₀ = 1.22, which is achieved 4000 years after splitting is stopped, at which stage a foam texture has been established. Just after stopping the splitting, the growth rate is much higher, corresponding to k ₀ = 4.2. (b) Detailed plot of the experiment in Figure 7a after 6500 years (splitting has been stopped). Equation (1) has been fitted to the experimental curve, giving apparent k and n values that are incorrect: n _app is 2.79 instead of 2 and k _app ± is 6.29 × 10⁻³ instead of 3.90 × 10⁻³ mm² a⁻¹.

Discussion

The modelling in this paper is in no way intended to argue that the microstructure and grain size of the upper hundreds of metres of polar ice is determined by a balance of static grain growth and a constant grain-splitting rate. For this reason, we do not attempt to fit the results of the numerical simulations to obtain a growth constant or an average split rate of once every so many years. The dynamics of rotational recrystallization are much more complex (Reference Faria and KipfstuhlFaria and Kipfstuhl, 2004; Reference Weikusat, Miyamoto, Faria, Kipfstuhl, Azuma and HondohWeikusat and others, 2011) than can be grasped by a simple constant split rate that is equally applied to all grains.

The intention of this paper is to show one of the pitfalls of numerical simulations that do not include the effect of microstructure. The simple analytical model of growth versus splitting produces a curve that can be fitted to data from ice cores. At first sight, it appears that the use of a simple splitting constant, f, would be the most problematic simplification. However, our simulations show another simplification that is rarely considered, namely lack of coupling between the growth ‘constant’, k ₀, and f. The parameter k ₀ is determined by the microstructure. As the microstructure is a variable, k ₀ is not a constant, but a variable as well. This observation is of importance because many models that incorporate grain growth, assume k ₀ to be constant (Reference Cotterill and MouldCotterill and Mould, 1976; Reference Randle, Ralph, Hansen, Hansen, Juul Jensen, Leffers and RalphRandle and others, 1986; Reference Montagnat and DuvalMontagnat and Duval, 2000). The numerical simulations show that changing the grain-size distribution from normal to approximately log-normal increases k ₀ by a factor of ∼3.5. Clearly, other factors may influence k ₀, such as grain boundary morphology and grain shape.

The simulation of static grain growth shows that the resulting grain-size distribution is relatively narrow. A normal distribution of measured grain diameters is predicted for static grain growth (Reference Humphreys and HatherlyHumphreys and Hatherly, 1996). However, grain diameter distributions in ice are usually lognormal, even at relatively shallow depths (Reference Arnaud, Gay, Barnola and DuvalArnaud and others, 1998), for example at 115 m depth in the NorthGRIP core (Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997; Reference SvenssonSvensson and others, 2003). This indicates that the microstructure of ice is already strongly affected by processes other than only static grain growth, well above the transition to a stable grain size. This observation supports the suggestion by various authors (Reference KipfstuhlKipfstuhl and others, 2006, Reference Kipfstuhl2009; Reference Durand, Perrson, Samyn and SvenssonDurand and others, 2008; Reference Weikusat, Kipfstuhl, Azuma, Faria, Miyamoto and HondohWeikusat and others, 2009a, Reference Weikusat, Kipfstuhl, Faria, Azuma and Miyamotob) that dynamic recrystallization and other processes (Reference Arnaud, Barnola, Duval and HondohArnaud and others, 2000; Reference Faria, Freitag and KipfstuhlFaria and others, 2010) already commence at relatively shallow depth.

The observation that k ₀ is dependent on the microstructure may have consequences for the interpretation of grain growth experiments to determine the growth exponent n. If the experiment is started with a non-equilibrium microstructure, k ₀ may initially be much higher. As the microstructure stabilizes to that characteristic of static grain growth, k ₀ decreases (Fig. 7). If the initial phase of microstructural equilibration is included in an analysis where k is assumed to be constant, one would erroneously obtain an exponent n that is larger than the real value. For example, the applicable value for k ₀ in a polar ice cap would be different from one obtained in a static grain growth experiment, because the microstructure, and hence grain growth in nature, is influenced by additional factors, such as dynamic recrystallization, presence of impurities and bubbles (Reference Cuffey, Thorsteinsson and WaddingtonCuffey and others, 2000).

Conclusions

We simulated the process of pure static grain growth and grain growth in competition with another process, namely splitting grains at a constant rate. The numerical simulations show that the growth parameter k ₀, normally taken to be a constant, is in fact a function of the microstructure. When the microstructure is only affected by static grain growth, k ₀ is 1.22. The change in microstructure resulting from additional splitting increases k ₀ by a factor of ∼3.5.

The numerical simulations show that the log-normal grain-size distributions observed in polar ice at shallow depth (≥100 m) are not in accordance with the expected distributions for static grain growth. At least one other process must operate to widen and skew the distribution towards a log-normal distribution. This supports the idea that dynamic recrystallization already operates and influences the microstructure at shallow depth.

The growth exponent and grain boundary properties (surface energy and mobility) are usually determined from experimental growth curves. If the microstructure changes during these experiments, k ₀ should not be assumed constant. Making this assumption leads to an overestimate of the growth exponent n.

Acknowledgement

We gratefully acknowledge funding by the German Research Foundation (DFG) project BO-1776/7.