Stacking disorder in ice I

Traditionally, ice I was considered to exist in two well-defined crystalline forms at ambient pressure: stable hexagonal ice (ice I h ) and metastable cubic ice (ice I c ). However, it is becoming increasingly evident that what has been called cubic ice in the past does not have a structure consistent with the cubic crystal system. Instead, it is a stacking-disordered material containing cubic sequences interlaced with hexagonal sequences, which is termed stacking-disordered ice (ice I sd ). In this article, we summarise previous work on ice with stacking disorder including ice that was called cubic ice in the past. We also present new experimental data which shows that ice which crystallises after heterogeneous nucleation in water droplets containing solid inclusions also contains stacking disorder even at freezing temperatures of around (cid:2) 15 1 C. This supports the results from molecular simulations, that the structure of ice that crystallises initially from supercooled water is always stacking-disordered and that this metastable ice can transform to the stable hexagonal phase subject to the kinetics of recrystallization. We also show that stacking disorder in ice which forms from water droplets is quantitatively distinct from ice made via other routes. The emerging picture of ice I is that of a very complex material which frequently contains stacking disorder and this stacking disorder can vary in complexity depending on the route of formation and thermal history.

In the past, metastable ice I had been identified as having a cubic crystal structure with space group Fd% 3m, whereas the stable structure of ice I is hexagonal with space group P6 3 /mmc. 36Both of these forms of ice I are made up of layers composed of sixmembered puckered rings of hydrogen-bonded water molecules and only differ in the way the layers are stacked on top of each other (see Fig. 1).In ice I h , each layer is a mirror image of the previous layer; whereas in ice I c each successive layer is shifted a distance equal to half the diameter of a six-membered ring. 36,37owever, it is becoming increasingly clear that the ice I phase which forms in many experiments does not fit either of these two well-defined crystal structures.Instead, metastable ice I is typically made up of a combination of both cubic and hexagonal stacking sequences which together do not poses cubic nor hexagonal symmetry (Fig. 2) and it has the trigonal space group P3m1. 9,12,13,30In order to emphasise the fact that metastable ice I is neither cubic nor hexagonal, and not a simple mixture of the two, Malkin et al. 30 recommended calling this material stacking disordered ice I (ice I sd ).We use the term stacking disordered ice I or ice I sd in this paper for metastable ice I which contains stacking disorder and reserve the term cubic ice (ice I c ) for the hypothetical form of ice I with a cubic structure as envisaged by Ko ¨nig 24 more than 70 years ago.
In the mid 1980s Kuhs et al. 14 were the first to appreciate that their samples of ice I, which they referred to as ice I c , contained appreciable amounts of stacking faults.Specifically, they suggested that features in their neutron diffraction patterns were due to a small probability of finding hexagonal sequences in a structure dominated by cubic sequences.These hexagonal sequences were thought of as stacking faults in a cubic structure, but this early study lacked any quantitative analysis of the density of these stacking faults.At the time, Elarby-Aouizerat et al. 38 suggested that while stacking faults could account for some non-cubic features, they could not account for the all the non-cubic features in the diffraction patterns.With the development of a detailed stacking fault model of ice, Hansen et al. 12,13 were able to quantitatively show that all the non-cubic features in ice I formed from ice V and ice IX could be accounted for with stacking disorder.They also showed that correlations between stacking faults was different dependent on the route of formation.Morishige and Uematsu 16 also showed through modelling diffraction patterns that ice in mesopores contained significant stacking disorder.After posing the question of 'is it cubic?',Moore and Molinero 39 arrived at the conclusion that metastable ice I contains comparable fractions of cubic and hexagonal layers.Subsequently, Malkin et al. 30 showed that ice resulting from the homogeneous freezing of water droplets at around À40 1C was made up of a fully randomised mixture of 50% cubic and hexagonal sequences.More recently, Kuhs, et al. 9 showed that ices generated in multiple ways, including through vapour deposition and recrystallization of gas hydrates, all contains stacking disorder to varying degrees and possesses a trigonal, rather than cubic, space group.In summary, there is compelling evidence that ice I has a propensity to be stacking disordered.
In this perspectives article we initially review the diffraction data in the literature, which shows that the non-cubic features of the diffraction patterns associated with stacking disorder are ubiquitous to all metastable ice I that we are aware of.The data also suggests that there are strong differences in the nature of the stacking disorder observed depending on method of generation.0][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58] We present new experimental data for ice formed from water droplets in which freezing was initiated heterogeneously.In order to do this we built on our knowledge of heterogeneous ice nucleation, [59][60][61][62][63] in order to select materials which nucleate ice over a range of temperatures between 237 and 263 K.These results are unique in that they allow us to systematically probe the evolution of ice structure as a function of ice nucleation temperature.In addition we also present new data for ice I sd formed via the recrystallization of ice II.We then model this diffraction data to quantify the number density and nature of cubic and hexagonal stacking sequences, and summarise this information on a so-called ''stackogram''.This new type of plot helps to visualise the differences in the nature of stacking disorder in ice I sd made through different routes.

The ubiquity of stacking disorder: evidence from diffraction patterns
Diffraction patterns provide information on long range order (and disorder) in crystal structures.A diffraction pattern of a welldefined crystalline material contains sharp Bragg peaks, the position of which can be related to distance between adjacent planes in a crystal structure known as d-spacings.The Bragg condition is met in an ordered crystal where the long range order leads to constructive interference satisfying Bragg's law and giving rise to sharp diffraction peaks.If for some reason, such as the introduction  of stacking disorder, the long range order in a crystal is interrupted this influences the resulting diffraction pattern.
The calculated diffraction patterns of well-defined ice I h and ice I c are shown in Fig. 3 for both neutron and X-ray diffraction.Together with these patterns we also plot literature diffraction patterns for ice I which was formed in a variety of ways.Table 1 summarises the experimental conditions and routes through which this ice was formed.These patterns clearly do not correspond to hexagonal ice and in the past it has been generally concluded that this ice is therefore cubic.However, comparison with the perfect cubic pattern also reveals a poor match.The literature patterns clearly contain a peak at a d-spacing of 3.86 Å which is in the same position as the most intense peak in the hexagonal ice pattern.These experimental patterns of metastable ice I also have a region of raised intensity between 3.43 and 3.86 Å d-spacing and in addition the relative ratio of the peak intensities does not match the calculated patterns.Inspection of high resolution patterns also reveals other less obvious inconsistencies. 12,30To varying degrees the experimental patterns in Fig. 3 all contain features inconsistent with the cubic crystal system and to the best of our knowledge there is no diffraction pattern of ice which is a good match to that expected for well-defined cubic ice.Instead, metastable ice I always contains features associated with stacking disorder.

Stacking disorder in ice from a computational perspective
The propensity of ice I to grow with stacking-disorder has been observed in a large number of molecular dynamics and Monte Fig. 3 A selection of (A) X-ray and (B) neutron diffraction patterns of metastable ice I from the literature.Many of these patterns were digitised from the printed plots in the respective papers, and the quality of some of the older reproductions is relatively poor.Nevertheless, signatures of stacking disorder are apparent in all of these patterns, many of which were presented as ice I c .6][57] The study of the structure of ice during nucleation and growth through molecular simulations faces three main challenges.First, the accuracy of the water models in describing the ice-liquid equilibrium temperature for ice I, and the relative stability of the cubic and hexagonal ice polymorphs.Of the most popular fully atomistic models of water, only those in the TIP4P family (TIP4P, 67 TIP4P-Ew, 68 TIP4P/2005, 69 TIP4P/ice 70 ) predict that ice I is the most stable crystal at ambient pressure.70,71 Of these, only TIP4P/ice predicts a melting temperature for ice I h close to the experimental value, 272.2 vs. 273.15K; the others underestimate T m by at least 20 K. [69][70][71] To our knowledge, the relative stability of hexagonal (I h ) and cubic (I c ) ice I polymorphs has not been established for any of these models.A significant issue related to the use of fully atomistic models is their computational cost, which limits the sampling time and size of simulation systems that can be modelled. The coarse-rained monatomic water model mW was developed to alleviate the costs of simulations of liquid water and ice, while keeping an accurate description of the structures and phase transitions between them.72 mW models water as a single particle that interacts through short-range twobody and three-body interactions; the latter encourage tetrahedral configurations that mimics the effect of hydrogen bonds at less than 1% of the computational cost of atomistic models.72 Hexagonal ice is the most stable crystal phase of mW water at ambient pressures, with a melting temperature of 274 K. 72 Cubic ice is marginally less stable than hexagonal ice for the mW model.40,65,72 The second significant challenge arises from the rare nature of ice nucleation events, which calls for very long simulation trajectories or the use of advanced simulation methods.Ice nucleation by long brute force simulations has been achieved for mW water under a broad range of supercooled conditions 40,49,52,[73][74][75][76][77][78][79][80][81][82] and for TIP4P water supercooled at negative pressures.83 Among the advanced methods to sample rare events, umbrella sampling, 84,85 metadynamics [86][87][88] and forward flux sampling [89][90][91][92] have been successfully used to nucleate ice at various degrees of supercooling and with several water models.50,[93][94][95][96][97][98][99][100][101][102][103] The ice polymorphs obtained using different sampling methods are not always identical, even when the model and simulation conditions are the same.94,95 Differences in outcomes may have two origins: first, some are equilibrium sampling methods (e.g.umbrella sampling and metadynamics), while others are intrinsically non-equilibrium (e.g. forward flusampling); second these methods have been implemented with different order parameters used to bias or identify the formation of ice.This leads to the third challenge for the study of the structure of ice formed in molecular simulations: the need for methods to identify water molecules as liquid, cubic ice and hexagonal ice.Several algorithms have been proposed for this purpose.Brukhno and co-workers introduced the first set of order parameters that distinguished cubic and hexagonal ice, however these were not rotationally invariant hence they cannot be used to identify arbitrarily oriented crystallites.50 Moore and co-workers then introduced CHILL, 40 a rotationally invariant algorithm based on the correlation of bond-order parameters, 104,105 which distinguishes cubic ice, hexagonal ice and liquid water from the number of staggered and eclipsed OÁ Á ÁO bonds.40 Variants of the latter approach using bond-order parameters were also introduced by Li et al., 96 and Sanz et al.Unbiased simulations of homogeneous ice nucleation with the coarse-grained mW water model 52,107,108 reveal that the ice embryos already have cubic and hexagonal ice sequences.39,40,[42][43][44][45] The same was found in the nucleation of ice using the atomistic TIP4P/2005 and TIP4P water models combined with advanced sampling methods that do not bias the structure of the embryo.30,46,50 Ice embryos need to have at least B200 water molecules for stacking to be discernible, but they contain mixed cubic and hexagonal ice sequences from their inception (Fig. 4).39 The fraction of cubic stacking sequences (cubicity) of the crystallites nucleated from deeply supercooled water increases with their size, plateauing at about 60%.39,43 The ratio of cubic to hexagonal sequences is similar for ice nucleated from bulk water, solutions with salt, and water confined in nanopores and in nanoparticles.39,40,[42][43][44]47,49 The stacks of cubic and hexagonal sequences in the ice nucleated and grown from deeply supercooled water are very short and their order seems to be random.30,39,50 These stacks reorganize as the crystallites consolidate in the process of growth.39 The evolution of the stacking disordered structures into hexagonal ice is outside of the time scale accessible to brute-force molecular simulations.
The free energies of the stacking disordered and hexagonal ices have been reported to be within B100 J mol À1 , 18,25,39,40,51,52,107,[109][110][111] and calculations indicate that the free energy barrier for the creation of a critical embryo that is purely cubic or purely hexagonal is lower than for the creation of a stacking-disordered embryo. 4544]47,[49][50][51][55][56][57] It is possible to grow ice I h from a hexagonal ice seed and ice I c from a cubic seed using Umbrella Sampling simulations that very slowly bias the size of the embryo, sampling equilibrium configurations for each embryo size. 45These results suggest that stacking disorder has a kinetic origin: it is controlled by the non-equilibrium process of ice growth.A recent study, however, suggests that stacking disorder is thermodynamically favored. 656][57] A configuration of a large-scale molecular dynamics simulation of ice grown at 270 K from liquid mW water is shown in Fig. 5.The simulation started with a slab of hexagonal ice exposing the basal planes to the supercooled liquid; 41 of the 81 new layers formed have cubic order.Interestingly, the cubicity of the ice grown at 270 K is lower than for ice grown at higher supercooling. 51,57A recent study of the equilibrium interface between ice and liquid in TIP4P/2005 water shows extensive reconstruction of the interface of the originally hexagonal ice slab, resulting in a stacking disordered ice interface with 60% cubic ice and 40% hexagonal ice in contact with liquid water. 66A systematic study of ice growth at different temperatures and a characterization of the relative stabilities of cubic and hexagonal ice for the water models employed would be needed to assess the temperature dependence of the cubicity of ice grown in simulations.
Ice tends to grow in a layer by layer fashion in the direction perpendicular to the stacking plane, (111) of ice I c or (0001) of ice I h . 51,55,57Recent studies 57 reveal that hexagonal and cubic arrangements, which only differ by a small displacement of the molecules with respect to the underlying layer (Fig. 1), occur with similar probability on the basal surface of ice, and their competitive formation and dissolution slows down the growth rate of ice in the direction perpendicular to the basal plane. 57his in-layer competition does not occur in the prismatic and secondary prismatic planes, 51,[55][56][57]112,113 which only produced stacking faults in simulations at very high supercooling. 53 Formtion of stacking disordered ice can be considered a case of extensive cross-nucleation [114][115][116] between ice polymorphs.Cross-nucleation involves the nucleation of one crystal structure (polymorph) on the face of another, and usually favours the faster growing polymorph, irrespective of whether it is the most stable crystal.116,117 Cross-nucleation between polymorphs that The lines connect pairs of molecules with cubic or hexagonal order, which only differ by a small displacement of the molecules with respect to the underlying layer, see Fig. 1.Pairs of molecules with cubic and hexagonal ice order are distinguished by the number of eclipsed and staggered intermolecular bonds within the pair, identified with the CHILL algorithm: 40 molecules in hexagonal ice have three staggered and one eclipsed bond, and molecules in cubic ice have four staggered bonds. We-defined stacks are not observed in embryos with less than B200 molecules.have a common plane, as is the case for (111) of ice I c and (0001) of ice I h , is very facile as it does not require the nucleation of an interfacial nucleation layer to seamlessly connect the two polymorphs.117,118 Clathrate hydrates provide other examples of cross-nucleation between water crystal polymorphs. 117,118The main differences between cross-nucleation in ice and in clathrates is that not only the stabilities of ice I h and the hypothetical ice I c polymorphs are almost identical, but that they apparently also have similar rates of growth.
Although cubic and hexagonal ice layers can stack seamlessly, in-layer competition between hexagonal and cubic order can result in the formation of defects when domains with distinct order propagate in the same layer. 57Large-scale simulations of ice growth revealed the presence of lines of defects, coupled 5-and 8-membered water rings, in the boundary between domains of cubic and hexagonal order coexisting within a single layer. 53It has been proposed that these defects facilitate the nucleation of stacking faults. 53A recent study, however, suggests that the 5-8 defects are not the origin but the result of in-layer competition between cubic and hexagonal order in the newly formed layer of ice. 57The density of these defects should increase with supercooling as the rate of nucleation of new layers on the growing ice surface becomes faster.The role that these defects and their annealing play on the long-term decrease of stacking disorder in ice is an open question that warrants further investigation.
Heterogeneous nucleation of ice has been only recently achieved in molecular simulations. 48,58,119The two surfaces studied, graphitic carbon and kaolinite, produce distinct preferential orientation of the critical ice embryos on the surface.Ice embryos expose the basal plane to the graphite surfaces, resulting in the growth of stacking disordered ice with the stacking plane parallel to the surface. 48Kaolinite, on the other hand, favours the formation of hexagonal ice embryos that expose the prismatic plane to the mineral surface. 58These embryos, however, are finite in size and expose all faces to liquid water.This implies that, upon growth, they can develop stacking faults.The experimental results in the next section show that ice nucleated by kaolinite presents significant stacking disorder, supporting the evidences from molecular simulations that the non-equilibrium growth process plays an important role in the development of stacking disorder in ice.

Heterogeneous freezing of water droplets
In the past, the structure of ice crystallised after homogeneous nucleation of pure water droplets at around À40 1C, 27,30 and solution droplets below À40 1C, 6,[27][28][29]31,120 have been investigated. Howeve, to date the structure of ice which crystallises following heterogeneous nucleation in water droplets has only been investigated through molecular simulations, 48,58 and has not been subject to laboratory experiments.In this section, we present experiments where the aim was to quantify the degree of stacking disorder in ice which crystallises in droplets where nucleation was induced by well characterised heterogeneous ice nucleants.This builds on our previous research on the quantification of the efficiencies with which various solid materials nucleate ice when immersed in water droplets.59- 62,121 These experiments involved cooling water droplets containing solid nucleating agents at a controlled rate, measuring the freezing temperature and subsequently recording a diffraction pattern of the ice which crystallised.The technique is similar to that described previously, 26,27,30,31 but with the exception that we now work with droplets containing nucleants.Fine powders of a variety of minerals and other solids were suspended in ultra-pure water (18.2 MO). The aqueous suspensions were mixed with paraffin oil (Fisher Scientific) and lanolin (Aldrich Chemical Company) in order to create water-in-oil emulsions with droplets of a volume median diameter, d vm , of 17 AE 3 mm.This was significantly larger than the 0.9 mm used by Malkin et al. 30 and ensured that each droplet contained a representative number of sub-micron particles.The droplets sizes were then determined using an optical transmission microscope with a 10Â objective.The X-ray diffractometer (Bruker D8 Advance, Cu Ka) used in these experiments was configured in standard 2y reflection geometry and was equipped with an Anton Paar TTK450 temperature control stage.
In order to determine the temperatures at which the water droplets froze, the diffraction angle (2y = B401) of a strong reflection associated with all types of ice I was continually monitored upon cooling at 30 K min À1 .This peak is insensitive to the phase of ice, only varying in intensity by 1.9% between ice well-defined I h and ice I c and being no larger for ice I sd (according to our calculated patterns, which is consistent with our measurements) and therefore provides a useful proxy for the amount of ice crystallised.The cooling rate of 30 K min À1 was chosen on the basis that we wanted to quantify stacking-disorder in ice which initially crystallised, and slower cooling rates allow more time for recrystallization of stacking-disorder.The area under the peak at B401 was normalised to the peak area at 173 K where all the droplets were frozen and the peak area was at its maximum.Plots of the fraction of pure water droplets (0.9 and 17 mm) and droplets containing kaolinite (B17 mm) frozen, as a function of temperature during cooling, are shown in Fig. 6A.
It is clear from Fig. 6A that the addition of kaolinite particles induces freezing at higher temperatures than droplets containing no particles (the median freezing temperature increases by 1 to 5 K depending on the weight fraction of kaolinite).This means that the mode of nucleation shifts from homogeneous to heterogeneous.In addition, the solid lines in Fig. 6A are the freezing curves predicted on the basis of literature parameterisations; 60,121 the excellent agreement indicates that the presence of the oil and surfactant did not influence nucleation.The fraction frozen curves for droplets containing a range of other ice nucleants that triggered heterogeneous nucleation with median freezing temperatures as high as 257 K are shown in Fig. 6B.
Diffraction patterns of the sample between 2y = 20-701 were recorded at 173 K, a temperature at which stacking-disorder was assumed to not change on the timescale of the diffraction This journal is © the Owner Societies 2015 measurement.The range of 2y covers all of the strong ice I c and ice I h reflections, and the diffraction pattern for frozen droplets containing 0.1 wt% kaolinite is compared with the diffraction pattern for ice resulting from homogeneous freezing (for 0.9 mm droplets 30 ) in Fig. 7. Also shown are the calculated ice I h and ice I c patterns.Despite ice nucleating heterogeneously on kaolinite, the ice I sd produced is almost indistinguishable to that seen after homogeneous nucleation, which was previously shown by Malkin et al. 30 to have fully random stacking.Utilising molecular simulations, Cox et al. 58 found that the prismatic face of hexagonal ice nucleates on kaolinite.Our results suggest that even if nucleation of one phase or another is preferred, the growth of stacking disordered ice subsequent to nucleation is controlled by kinetic factors during crystal growth.
A fully random stacking sequence is consistent with a twodimensional nucleation (or layer-by-layer growth) mechanism, 122 where each successive layer nucleates independently and can either stack in a cubic or hexagonal fashion.Our result indicates that the probability of nucleating a cubic or hexagonal sequence is equal.This also implies that the crystal structure of the macroscopic frozen droplet is independent of the structure of the initial critical cluster that nucleates.This is important, because it is often assumed that nucleation of one particular phase will define or serve as a template for the resulting crystal growth. 1 This has implications for the interpretation of the analysis of nucleation data using classical nucleation theory (CNT).
Several authors used experimentally determined nucleation rates to establish ice-water interfacial energies on the basis that ice I c nucleates. 32,121,123,124Murray et al. 121 justified this assumption on the basis of X-ray diffraction data reported in a separate study, 27 which clearly contain features consistent with stackingdisorder.Huang and Bartell 32 used electron diffraction data to identify their nanometre sized ice particles as ice I c , but also note that there were 'imperfections in the longer-range internal order', which may also be consistent with stacking-disorder.In order to employ classical nucleation theory to derive the icewater interfacial energy, these authors used the best available estimate of the thermodynamic properties of ice I c . 32,121,123,124iven the intrinsic quantitative inaccuracies of CNT, the additional assumption that the critical cluster has the same structure as the final macroscopic crystal is in question, these interfacial energies should be taken with some caution.
Previous studies have shown that ice I sd transforms to ice I h increasingly rapidly as temperature is elevated. 8,13,27For droplets suspended in an oil emulsion, where the only recrystallization mechanism is a solid-to-solid transformation, as opposed to a vapour or solvent mediated route, 31 the time required for recrystallization increases rapidly below about 238 K. 27 Hence, during  cooling at 30 K min À1 droplets which froze below this temperature did not have the opportunity to transform to the stable ice I h , locking the ice in the metastable structure in which it initially crystallised.We hypothesise that ice which initially forms at higher temperatures is also stacking-disordered, but freezing at higher temperatures may provide an opportunity for the ice to anneal.In order to explore this hypothesis, solid particles of different types with the capacity to nucleate ice at different temperatures were introduced into the droplets.The resulting diffraction patterns after complete crystallization are shown in Fig. 8.The patterns are arranged in order of lowest freezing temperature at the top and a pattern of fully annealed ice I h at the bottom for comparison.It is evident from the observations that on increasing freezing temperature the peaks unique to ice I h at 26, 33 and 441 all increase in intensity, the broad feature between 22 and 271 disappears, all peaks become increasingly sharp and the relative intensities of the peaks generally approaches that of ice I h .This evolution of the diffraction patterns is consistent with higher freezing temperatures allowing the initially stacking-disordered ice to anneal leading to ice which is structurally closer to ice I h .

Modelling of stacking disorder in ice I
In order to quantitatively characterize the stacking disorder in ice I, we employed the DIFFaX (Diffracted Intensities from Faulted Xrystals, v1.813) computer program for calculating powder diffraction patterns of crystals containing stacking disorder. 125IFFaX uses a general recursion algorithm and requires information about the structure of the layer, the stacking probabilities and the symmetry relationships between the stacked layers.The calculated pattern is convolved with a profile function in order to account for finite crystallite size and instrument broadening effects.DIFFaX can be configured to take into account varying complexities in stacking disorder.The simplest form of stacking disorder in ice is where cubic and hexagonal sequences are randomly arranged.For example, Malkin et al. 30 showed that ice resulting from homogeneous nucleation in pure water droplets was fully random with 50% cubic and 50% hexagonal stacking.More complex stacking can also be possible where there is correlation or memory in the stacking sequence.For example, once in a hexagonal stacking sequence the probability of transitioning back to cubic may be small, but when it does there may also be a low probability of transitioning back again.This would result in extended clusters of cubic and hexagonal sequences.In this example, the correlations are between nearest neighbour layers (1st order memory effects; Kuhs et al. 9 define this as an interaction range, s, of 3), but there can also be correlations between sequences further apart (2nd, 3rd, etc. order memory effects; s = 4, 5 etc.).Kuhs et al. 9 and Hansen et al. 12,13 found that the ice they made through a variety of routes contained 2nd order (s = 4) memory effects.In what follows, we explore the possibility of complex memory effects.
Stacking disordered ice with no memory effects (s = 2) can be described in DIFFaX using a single stacking probability, F c , which indicates the likelihood of cubic stacking.The probability of hexagonal stacking, F h , is simply 100 À F c .To model 1st order (s = 3) memory effects, it is necessary to use two independent stacking probabilities, F cc and F hc which define the probabilities of cubic stacking after a previous cubic or hexagonal stacking events, respectively.The two hexagonal stacking probabilities, F ch and F hh , are obtained from 100 À F cc and 100 À F hc .In this  3.The MCDIFFaX fits here included 2nd order memory effects.The results were indistinguishable from the 1st order model indicating 2nd order memory effects were not significant.
paper we have modified the original terminology used by Malkin et al. 30 so that we can describe higher order memory effects (the relationship between the different notations for stacking probabilities is summarised in Table 2).
Since both, F cc and F hc , can vary independently between 0 and 100% there are a wide range of possible ice I sd structures.The relationship between ice structure, and F cc and F hc is illustrated in the ''stackogram'' shown in Fig. 9.If the two stacking probabilities F hc and F cc are equal then there is no 1st order memory in the stacking sequences of the ice structure and we refer to this as random stacking (red line in Fig. 9) which can be described with F c only.The top right of the stackogram (Fig. 9) represents perfect ice I c , while the bottom left describes perfect ice I h .In the case where F hc and F cc are not equal then there is a memory in the stacking sequences and there are a number of distinct regimes which can be defined.Anywhere to the right of the random line is a regime with the same close-packed planes tending towards the strictly alternating stacking sequence in the third dimension (polytype) at the bottom right of the stackogram (hexagonal, cubic, hexagonal, cubic etc.).As will be seen later no known ice sample falls in this half of the stackogram.The region above the random line is characterised by extended sequences of either cubic or hexagonal stacking, or extended regions of both.The top left of the diagram is for a sample consisting of a 50% mixture of bulk cubic and bulk hexagonal crystallites.
In order to illustrate the sensitivity of the diffraction pattern to stacking disorder and 1st order memory effects we have plotted an array of calculated diffraction patterns in Fig. 9.We focus on the region around the 111 peak of cubic ice where the pattern is particularly sensitive to stacking disorder.The red line indicates a combination of stacking probabilities that produces random stacking (F hc = F cc ).Above this line F hc o F cc whereas, below F hc 4 F cc .One striking feature of this plot is the strong sensitivity of the diffraction pattern to the values of F hc and F cc , which means that finding a unique fit to the diffraction pattern is highly likely.This is illustrated in Fig. 10 where the 1st order memory model was used to produce patterns for many combinations of F hc and F cc , and the goodness of the fit is expressed as w 2 .There is a clear single minimum and no local minima.Malkin et al. 30 allowed F hc and F cc to vary independently, hence they made no a priori assumption about which type of stacking was present in their ice (i.e.random or 1st order memory effects).They found that for ice which nucleated homogeneously from pure water there was no memory effect and the resulting ice was fully stacking disordered (i.e.F hc and F cc = 50%).
The diffraction pattern is remarkably sensitive to even small degrees of stacking disorder.This is illustrated in Fig. 11 where droplets containing silver iodide froze around 257 K and the peaks associated with hexagonal character are clearly present.However, the pattern is clearly inconsistent with well-defined hexagonal ice and instead we find, F hc = 1% and F cc = 10%.Note that the fitted diffraction pattern from ice annealed at just below 0 1C has a F hc and F cc = 0% as expected for ice I h .This indicates that we can detect stacking fault probabilities on the order of only 1% in ice I h .
Kuhs et al. 9 concluded that the ice they made via decomposition from gas hydrates and vapour deposition contained 2nd order memory effects (s = 4).Hence, it is prudent to analyse  the new data with a DIFFaX model capable of resolving second order memory interactions.In order to resolve 2nd order memory effects (s = 4), a model with four independent stacking probabilities (F ccc , F hcc , F chc and F hhc ) is required in DIFFaX.The relationships between these stacking probabilities and the parameters defined by Hansen and co-workers, 9,12 are given in Table 2.
In the model employed by Hansen et al. 12,13 and Kuhs et al. 9 the best fit to their diffraction data was found using a leastsquare minimisation approach.To obtain the best possible fit to our diffraction data we have written a new computer programme (MCDIFFaX) which embeds DIFFaX in a least-square environment and uses a Monte Carlo algorithm to search for the best values of the various stacking probabilities, lattice constants, peak profile parameters and zero shift. 126Our fitting strategy is to start off with F c = 0.5 and no memory effects.Once the w 2 has converged, first and second order memory effects are introduced successively as the programme runs to see if this is required to improve the quality of the fit to the data.
In addition to the least-square approach we have used a grid sampling approach where the four stacking probabilities in the 2nd order memory model were systematically varied independently across their entire ranges from 0 to 100%.This test allowed us to verify the approximate location of the global minimum.This was initially done in large steps, for example 10%.Once a rough indication of the best values was obtained, a higher resolution run (e.g.steps of 1%) was performed over a smaller range of numbers to improve the accuracy of the result.The resultant structure file is compared to the observed structural results and scored using a delayed cross correlation function.The cross correlation function allows the identification of all possible stacking ratio fits.The maximum value of the cross correlation is used as a score of how similar the results are, where score = 1 À maximum correlation.This produces a 4D array of scores between 0 and 1.The dimensions of the array are determined by the number of steps in (F ccc , F hcc , F chc and F hhc ).The first minimum found in this 4D array is the lowest score present.An algorithm was developed to locate all of the minima.The current minimum value is initially set to the lowest minimum score and the cell (set of parameters, F ccc , F hcc , F chc and F hhc ) where the lowest minimum score occurs is marked as having been visited.All not yet visited cells are then searched to see if they contain a score between the current minimum value and a small positive increment to it.Any cells found that are within Fig. 11 Comparison of the experimental X-ray diffraction patterns of frozen water droplets containing silver iodide and a DIFFaX prediction for stacking probabilities of F hc = 1 or 2%, F cc = 10% and ice I h .Inset is a characteristic ice I h peak at B33.51 which is sensitive to the stacking ratio.
This journal is © the Owner Societies 2015 two neighbouring cells in any direction of an already visited cell are assumed to belong to that minimum peak.If a cell is found that is further away and not touching any neighbouring cells is considered to be another (second) minimum peak and is recorded and marked as visited.The current minimum value is then incremented by the small positive increment and the process repeated until the current minimum value exceeds 1. Fitting all the diffraction patterns using this technique showed that there was a single minimum in the parameter space and confirmed the result of the least squared routine.
The 2nd order memory MCDIFFaX fits to the diffraction patterns are shown in Fig. 8.We also fitted a 1st order memory model to the same data and the resulting stacking probabilities are shown together with those for the 2nd order memory model in Table 3.The correlations between the various stacking probabilities in Table 3 are illustrated in Fig. 12.In the case where 2nd order memory effects were not deemed important there should be no dependence on the structure of the sequence two stacking events away, i.e.F hcc = F ccc and F chc = F hhc .Hence, it is clear from Fig. 12A that there is no substantial 2nd order memory effect.The only pattern in which there may be a minor 2nd order memory effect is for droplets containing microcline which froze around 259 K and for which we found F hcc = 20% and F ccc = 24%.The DIFFaX fit to the experimental homogeneous nucleation diffraction pattern in ref. 30 showed no indication of 1st order memory effect, this has been corroborated by re-analyzing the data with MCDIFFaX.In contrast, Kuhs et al. 9 and Hansen et al. 12 found strong differences between their equivalent parameters indicating complex stacking patterns in ice made in a variety of different ways.Given, that Kuhs et al. 9 and Hansen et al. 12 show strong 2nd order memory effects in the ices recrystallized from high pressure phases, the possibility of even higher order memory effects should not be discounted.
Since there is no significant 2nd order memory in ice nucleated homogeneously or heterogeneously from pure water the first order variables F hc and F cc are therefore adequate to describe this ice.This is demonstrated in Fig. 12B, where the F hc and F cc from the 1st order memory fitting procedure and the same values derived from the 2nd order memory model are compared (the stacking probabilities from the 2nd order memory model are converted to 1st order equivalents using the relationships set out in Table 2).Hence, the remaining discussion of stacking in ice from water droplets is based around the 1st order memory model, and the F hc and F cc parameters.
At median freezing temperatures below B237 K, the ice is fully stacking-disordered (F hc and F cc E 50%) and indistinguishable from ice resulting from homogeneous nucleation in 0.9 mm sized droplets (see Fig. 7 and 8).Fig. 13 shows that as the freezing temperature increases, F hc and F cc decrease and the proportion of cubic sequences decreases.One striking result from this study is that there is still detectable stacking disorder even at the highest studied median freezing temperature of 257 K.In this instance F hc = 1% and F cc = 10%.This corresponds to approximately one cubic sequence (or stacking fault) in 100 layers of what is dominantly ice I h .
At freezing temperatures above 249 K, the stacking probabilities F hc and F cc are no longer equivalent, indicating that the stacking disorder is no longer random (s = 3).This is illustrated in the stackogram in Fig. 14 where the points corresponding to freezing temperatures higher than 249 K fall above the random stacking line.This suggests that domains of hexagonal ice form with higher freezing temperatures, together with smaller domains of stacking disordered ice.The monotonic decrease in the fraction of cubic sequences with increasing freezing temperature is consistent with the hypothesis that the initial ice to crystallize is ice I sd which anneals to stable ice I h if the experimental conditions allow.Hence, the structure of the resulting ice is a function of the temperature at which nucleation occurred, and the time available for re-crystallization, rather than the mode or Table 3 The average stacking ratio and corresponding goodness of fit (w 2 ) or score for each set of at least four heterogeneous nucleation experiments with a cooling rate of 30 K min À1 and a d vm of B17 mm.T onset is the temperature ice was first identified and T f is the median freezing temperature  ) is defined as the sum of the squared difference between model and experimental pattern at each point in 2y divided by the number of degrees of freedom.c Score using a delayed cross correlation function: score = 1 À maximum correlation (http://paulbourke.net/miscellaneous/correlate).Notation conversion and calculating F hc and F cc from MCDIFFaX model ratio can be found in nature of the nucleation event.Slower cooling rates had a similar effect to warmer freezing temperatures, with slower cooling rates leading to fewer cubic sequences.The lack of dependence of stacking-disorder on nucleation mechanism suggests that whenever ice grows from supercooled water in nature or technological applications it initially grows as ice I sd with randomly arranged 50% hexagonal and 50% cubic stacking sequences.Though, in many situations latent heat release may rapidly elevate the temperature to a regime in which stacking-disorder will anneal.These results offer detailed insights into the crystallization process of water, a process important in many fields, 127 and suggest that stacking-disorder plays an important transient role in the crystallization of ice from water in general.

Stacking disorder in ice I made via routes other than from liquid water
The equivalent 1st order memory stacking probabilities F hc and F cc for ice I made via the recrystallization of ice V and IX, 12,13 the decomposition of CO 2 hydrates 9 and via vapour deposition 9 are illustrated in Fig. 14.As discussed above, both Kuhs et al. 9 and Hansen et al. 12,13 report 2nd order memory effects, information which is not possible to plot on a 2D plot, hence the stacking probabilities reported by those authors have been converted to F hc and F cc .The details of the translation the parameters given by Hansen and co-workers to F hc and F cc can be found in Table 2. Additionally, we plot the stacking probabilities for ice from liquid water and new data for ice recrystallized from ice II.It is striking that all of the data in the stackogram (Fig. 14) is above of, or on, the random line (i.e.where F hc r F cc ).This is perhaps consistent with the fact that stacking disorder in ice   This journal is © the Owner Societies 2015 is a result of kinetic limitations in the growth of ice I and the stacking sequences never tend towards more ordered sequences which would be to the right (F hc 4 F cc ) of the random line.In general, ice which forms at higher temperatures tends to contain larger domains of hexagonal sequences.It is also evident that ice made in different ways seems to occupy different parts of the stackogram.Ice from liquid water is either at the random line for ice which nucleates at the lowest temperatures or just above the random line for ice which nucleated at warmer temperatures where domains of hexagonal sequences form.Whereas, ice formed from re-crystallisation of high pressure phases and gas hydrates has F cc than 50%.More work is required to understand what causes the variation in ice I sd generated through different pathways.
The results in Fig. 14 for ice I recrystallized from ice II were performed for this study.The ice II sample was prepared by freezing 1 mL of H 2 O in an indium gasket inside a Specac press die precooled to 77 K, followed by heating of the sample to 233 K at 0.4 GPa.The sample was then rapidly quenched to 77 K, recovered under liquid nitrogen and transferred onto an XRD sample holder under liquid nitrogen.This was subsequently mounted on a cold stage at 98 K of our X-ray diffractometer as described above.Quantitative Rietveld refinement showed that the ice II sample was contaminated with only 2-5% ice V. We then followed the recrystallisation of the ice II as the samples were warmed at rates of between 0.1 to 30 K min À1 in separate experiments.The samples were heated from 98-148 K (in 5 K increments) and then to 168 (in 2 K increments) at the defined rate.In order to record diffraction patterns of the ice at each increment without further transformation of the ice during the period of the measurement, the samples were cooled back to 110 K at each increment.Recrystallization to ice I sd was observed between 150 and 165 K.A selection of the resulting diffraction patterns for ice I sd , immediately after recrystallization, are shown in Fig. 15 together with the MCDIFFaX fits.The resulting F hc and F cc values are listed in Fig. 14 and Table 4.This ice falls well above the random line in a similar region to ice recrystallized from ice V and IX.Intriguingly, ice I sd from ice II is the closest to well-defined cubic ice of all the ices summarised in Fig. 14.
It is apparent from Table 4 and Fig. 14 that the stacking ratio of ice I sd recrystallized from ice II is dependent on the heating rate.There is a monotonic decrease in the fraction of cubic sequences with increasing heating rates, suggesting that increasing the rate of heating promotes the recrystallization of ice II to ice I h or the transformation of cubic sequences to hexagonal sequences.Thus, cubic stacking seems to be favoured kinetically.
The resulting stacking ratios from the 2nd order memory model, in Table 4, identify that there is a 2nd order memory effect since F hhc a F chc and F hcc a F ccc .The probability of appearance of a particular type of stacking sequence is affected by the next-nearest stacking sequences.The 2nd order memory grid sampling programme also identified multiple minima.However, upon further investigation and comparison of the numerically calculated fits with the experimental data the lowest scoring minima was the only reasonable fit; the best fits can be seen in Fig. 15.The ice II to ice I transition is endothermic and therefore there is no associated release of latent heat. 128This is perhaps why the sample forms ice which is furthest of any from the stable hexagonal structure and closest to perfect ice I c .

Symmetry of ice I sd
The six-fold rotational symmetry in crystals of ice I h is related to the six-fold screw axis in its crystal structure.Introduction of cubic sequences into a hexagonal structure disrupts the six-fold screw axis, but does not affect the three-fold rotational symmetry on the same axis.Hansen et al. 13 recognised this and identified the space group P3m1 for ice I containing stacking disorder.This space group is in the trigonal crystal system, a subgroup of the hexagonal crystal family.This was also  recognised in the past by Hallet et al. 129 who also showed numerous images of atmospheric ice crystals with three-fold rotational symmetry.A review of atmospheric trigonal ice crystals reveals that trigonal ice crystals occur over a wide range of atmospheric conditions and have been observed in clouds all around the globe both in the troposphere and stratosphere. 130

Conclusion
The emerging picture of ice I is of a highly variable and much more complex material than was typically thought in the past.All diffraction patterns for what was identified as cubic ice in the past contain features consistent with a substantial degree of stacking disorder and there exists no pattern indicative of welldefined cubic ice (to the best of our knowledge).We argue that the name stacking disordered ice, ice I sd , is the most appropriate name for this material.This name has now been used in the literature for several years, 4,30,58,62,123,124,[131][132][133][134] and we recommend it is used in place of the term cubic ice where the ice is stacking disordered.The term cubic ice should be reserved for well-defined ice I with cubic space group Fd% 3m.One striking feature of ice I sd is the fact that it is highly variable with the memory within the structure of ice I sd being strongly dependent on the route through which the ice was formed.Ice which forms homogeneously from pure water is fully stacking disordered, whereas ice recrystallized from high pressure crystal phases, for example, has a more complex structure in which there are correlations between nearest and next-nearest layers.The varying degrees and nature of stacking disorder raises important questions about the validity and applicability of the physical data for metastable ice I.Much of this data was collected without a record of the nature of the stacking disorder.6][137] It is possible that the broad range of enthalpies of transformation reported in these studies are in part related to the variability in stacking disorder in ice I sd , 11 although it may also reflect variability in other structural defects.In addition, it was very recently found that there are signatures in the Raman spectrum of stacking disorder and that the spectral features depend on the nature of the stacking disorder in ice. 134As a community, we should also consider revisiting measurements of thermodynamic quantities, (e.g.heat capacity and vapour pressure) as well as interfacial energies, optical constants and spectroscopic properties using well-characterized ice in order to understand the relationship between stacking disorder and the physical properties of ice I sd .
The new results presented in this work reveal that stacking disorder can be important in ice crystallized at much higher temperatures than previously thought.We showed that resulting ice heterogeneously nucleated at 257 K contains measurable stacking disorder.At these warm temperatures stacking disorder is likely to be transient, but nevertheless it appears to be an intermediate in the crystallisation of water to ice under conditions pertinent to much of the Earth's atmosphere.At below B200 K, which is relevant for stratospheric clouds, equatorial cirrus clouds and mesospheric clouds, stacking disorder will likely persist for extended periods of time and these clouds may even be made of ice I sd . 6,8,138This is consistent with the appearance of ice crystals in very cold cirrus with threefold symmetry, 139 since ice I sd has a trigonal space group.Intriguingly, ice crystals with three-fold symmetry are also observed in much warmer clouds (see Murray et al. 130 ).
Given what we have learnt about stacking disorder in ice, does well-defined cubic ice exist?We are not aware of any diffraction evidence for well-defined cubic ice.However, there are scattered reports which suggest that well-defined cubic ice may exist.A rare halo at 281 to the sun, called Scheiner's Halo, has been observed and is consistent with octahedral crystals of cubic ice. 140,141Others have suggested that angles between crystallites in snow crystals are consistent with growth off the faces of an initial crystal of cubic ice with an octahedral shape. 142Crystals sampled in the polar stratosphere sometimes have what appears to be 4-fold rotational symmetry consistent with crystals of cubic ice with a cubic habit, 143 and crystals grown from vapour deposition around 200 K in the laboratory have been observed with an apparent cubic morphology. 144Whilst we cannot rule out the existence or possibility of creating true cubic ice, there is an increasing realisation that there is strong propensity for ice I to grow with stacking disorder.

Fig. 1
Fig. 1 Stacking of layers in hexagonal (A) and cubic (B) ice.The vertical is normal to the (0001) basal surface of hexagonal ice, and the (111) surface of cubic ice.Only oxygen atoms are shown, connected by hydrogen bonds.

Fig. 2
Fig. 2 Possible stacking sequence in stacking disordered ice (ice I sd ).Only oxygen atoms are shown which are connected by hydrogen bonds.

Fig. 4
Fig. 4 Development of stacking disorder within an ice embryo in mW water crystallized at 180 K. Adapted from ref. 39.Panels (a) to (e) display the evolution of the cubic (red) and hexagonal (green) ice features within a growing embryo.The lines connect pairs of molecules with cubic or hexagonal order, which only differ by a small displacement of the molecules with respect to the underlying layer, see Fig.1.Pairs of molecules with cubic and hexagonal ice order are distinguished by the number of eclipsed and staggered intermolecular bonds within the pair, identified with the CHILL algorithm:40 molecules in hexagonal ice have three staggered and one eclipsed bond, and molecules in cubic ice have four staggered bonds.Well-defined stacks are not observed in embryos with less than B200 molecules.

Fig. 5
Fig.5Stacking disordered ice grown at 270 K. Water molecules are shown in green if they belong to a hexagonal sequence, red if they belong to a cubic sequence and in gray if they belong to the liquid.The simulation started with a slab of six layers of hexagonal ice (indicated by a black line below the picture) exposing the basal faces to the supercooled liquid.The simulations were performed with a periodic simulation cell containing 110 592 molecules modeled with the mW water model,107 which has a melting temperature of 274 K.107 The CHILL algorithm40 was used to identify cubic and hexagonal ice.Cubic layers represent 40 out of the 81 new layers grown.

Fig. 6
Fig. 6 Fraction of droplets frozen as a function of temperature.(A) Shows curves for droplets of pure water and droplets containing various weight fractions of kaolinite cooled at 30 K min À1 .Median droplet size was held almost constant in these experiments.The solid lines are based on timedependent parameterisations of homogeneous freezing of pure water 121 and heterogeneous freezing by kaolinite. 60(B) Shows curves for a range of solid nucleating agents.Error bars are omitted as they are comparable to symbol size.

Fig. 7
Fig.7Experimental diffraction patterns for water droplets containing 0.1 wt% of kaolinite, cooled at a rate of 30 K min À1 compared with frozen pure water droplets of 0.9 mm from Malkin et al.30 Simulated diffraction patterns using DIFFaX of well-defined ice I h and ice I c are also shown.The gaps in the experimental pattern are where diffraction peaks from the sample support are observed.

Fig. 8
Fig. 8 Experimental and model best fit X-ray diffraction patterns for frozen droplets containing a variety of solid inclusions (d vm B 17 mm), which froze over a range of median freezing temperatures (T f ).(A) Is an example of experimental data (crosses) and MCDIFFaX fit (solid line) for droplets containing microcline.(B) Contains experimental (solid line) and MCDIFFaX fit (dashed line) diffraction patterns for droplets containing a range of ice nucleants.The data was normalized to the largest peak.The fitted stacking probabilities are listed in Table3.The MCDIFFaX fits here included 2nd order memory effects.The results were indistinguishable from the 1st order model indicating 2nd order memory effects were not significant.

Fig. 9 A
Fig. 9 A ''stackogram'' containing an array of calculated diffraction patterns (d-spacing = 3.40-3.85Å) for various values of F hc and F cc .The pattern is for the values of F hc and F cc at the centre of each box.The red diagonal line is where F hc and F cc are equivalent, which represents the random stacking line.Deviations from this line represent stacking with 1st order memory.The nature of ice at each of the corners of the diagram is indicated.

Fig. 10
Fig. 10 Chi-squared (w 2 ) maps of DIFFaX fits to the experimental data as a function of F hc and F cc using a 1st order memory model, for droplets containing: (A) kaolinite 0.05 wt% and (B) microcline 0.8 wt.Both plots show a clear and unique minimum in w 2 .

Fig. 12
Fig. 12 Correlation plots between stacking probabilities for ice resulting from heterogeneous nucleation in water droplets.(A) Shows that there are negligible second order memory effects.(B) Shows that the 2nd order model produces the same F hc and F cc as the 1st order model.The blue dashed line is the 1 : 1 line.

Fig. 13
Fig. 13 Cubicity, F cc and F hc as a function of T f for droplets containing a range of heterogeneous ice nucleants; the diffraction patterns and MCDIFFaX fits are shown in Fig. 8.The fit through the cubicity data is for illustrative purposes only.

Fig. 14 A
Fig.14A stackogram illustrating the stacking probabilities (F hc and F cc ) of ice I resulting from heterogeneous nucleation in water droplets from this study as well as those from: Malkin et al.30 for homogeneous nucleation in water; Kuhs et al.9 who report stacking probabilities in ice I made from CO 2 hydrates and vapour deposition onto a cold substrate; and Hansen et al.12,13 who made ice I sd from ice V and IX.F hc and F cc would both equal 100 for hypothetical perfect cubic ice I c .

Fig. 15 X
Fig. 15 X-ray diffraction patterns of ice I sd recrystallized from ice II (solid lines) together with MCDIFFaX fits including 2nd order memory effects (dashed lines).

Table 1
Experimental details for literature diffraction data shown in Fig.3 14This journal is © the Owner Societies 2015

Table 2
Notation conversion

Table 4
The average stacking ratio and corresponding score for ice I sd recrystallized from ice II at different heating rates a Cubicity is calculated as in Table2.