Predicting the Kinetics of Ice Recrystallization in Aqueous Sugar Solutions.

The quality of stored frozen products such as foods and biomaterials generally degrades in time due to the growth of large ice crystals by recrystallization. While there is ample experimental evidence that recrystallization within such products (or model systems thereof) is often dominated by diffusion-limited Ostwald ripening, the application of Ostwald-ripening theories to predict measured recrystallization rates has only met with limited success. For a model system of polycrystalline ice within an aqueous solution of sugars, we here show recrystallization rates can be predicted on the basis of Ostwald ripening theory, provided (1) the theory accounts for the fact the solution can be nonideal, nondilute and of different density than the crystals, (2) the effect of ice-phase volume fraction on the diffusional flux of water between crystals is accurately described, and (3) all relevant material properties (involving binary Fick diffusion coefficients, the thermodynamic factor of the solution, and the surface energy of ice) are carefully estimated. To enable calculation of material properties, we derive an alternative formulation of Ostwald ripening in terms of the Maxwell-Stefan instead of the Fick approach to diffusion. First, this leads to a cancellation of the thermodynamic factor (a measure for the nonideality of a solution), which is a notoriously difficult property to obtain. Second, we show that Maxwell-Stefan diffusion coefficients can to a reasonable approximation be related to self-diffusion coefficients, which are relatively easy to measure or predict in comparison to Fick diffusion coefficients. Our approach is validated for a binary system of water and sucrose, for which we show predicted recrystallization rates of ice compare well to experimental results, with relative deviations of at most a factor of 2.

Introduction

When freezing foods or biomaterials, one generally forms a mixture of many micrometer-sized ice crystals dispersed in a freeze-concentrated liquid product, or matrix. During storage, these systems are generally not in a global equilibrium state but slowly evolve in time by a number of physical processes, collectively known as “recrystallization”.[1−6] Whatever the mechanism leading to recrystallization, the effect is a system whose characteristic length scale (e.g., the average crystal size) grows in time, leading to a coarser microstructure and an accompanying loss of product quality. To gain better control over quality losses requires a fundamental understanding—and predictive description—of the physical mechanisms and driving forces determining the rate of recrystallization. In this work we aim at developing such a predictive description, for a model system containing water, ice, and sugars.

If storage conditions are (approximately) isothermal, the main mechanism leading to recrystallization is Ostwald ripening, or coarsening,[7−9] which refers to the growth of large crystals at the expense of smaller ones so as to minimize the total surface energy of the system. Much of the theoretical understanding of this process dates back to the work of Lifshitz–Slyozov and Wagner (LSW),[10,11] in which was predicted the existence of a stationary regime for long coarsening times, characterized by a distribution of particle sizes that evolves self-similarly upon rescaling by a characteristic length scale (e.g., the average particle radius or critical radius). The characteristic length scale raised to some power was shown to scale linearly with time, the dynamics of the system thus determined by e.g.with ⟨a⟩ being the average particle size, K being the coarsening rate, and the integer n depending on whether growth is limited by surface attachment kinetics (n = 2) or diffusion through the bulk matrix (n = 3). In addition to a general dependence on several material parameters (e.g., the binary diffusion coefficient and surface energy), the diffusion-limited case requires the coarsening rate to carry an additional dependence on the volume fraction of dispersed-phase particles in the system, as neighboring particles influence each others’ diffusion field. Due to its assumption of a growth rate of particles within an infinite medium, the LSW approach is strictly valid only in the limit of zero dispersed-phase volume fraction and is therefore unable to capture this effect. Many extensions of LSW theory to finite volume fraction have been developed;[12−23] all predict the same temporal law (eq ), but with a coarsening rate K that increases with volume fraction. Other assumptions known to limit the application of LSW-type theories are the implicitly assumed spherical particle geometry, fixed particle positions, and absence of spatial correlations between particles. Various extensions that (partially) alleviate these assumptions have been developed; for a review, see ref (9).

Available experimental results for coarsening of ice in sugar solutions[24−27] (and in many real food systems or biosystems) qualitatively confirm eq with n = 3, implying diffusion-limited growth. Any attempts at arriving at a quantitative agreement between theory and experiments have been less conclusive, however, and show several contradictory results.[24,27] In the work of Sutton et al.,[24] for example, a modification of the LSW theory to finite dispersed-phase volume fraction was shown to result in coarsening rates that deviate from experiments by as much as 1 order of magnitude, while in the work of Budke et al.,[27] calculations based on the LSW theory were shown to compare reasonably well to experimental coarsening rates (at low volume fraction). We find the poor predictions in the work of Sutton are largely caused by the use of some simplifying assumptions on the thermodynamics of the system, stemming from the original treatment of LSW (i.e., matrix phase diluted in precipitating species, matrix phase modeled by ideal solution, and equal phase densities). While the approach of Budke is based on these same assumptions, it also involves incorrect use of the water self-diffusion coefficient instead of the binary (Fick) diffusion coefficient of water–sugar, leading to a fortuitous cancellation of errors in the predicted coarsening rate. We should note the use of incorrect diffusion coefficients is no exception but appears in several other applications of Ostwald ripening theory throughout the literature.[26,28]

Here we show that a rigorous application of Ostwald ripening theory leads to very satisfactory (within a factor of 2) predictions of experimental coarsening rates over a large and relevant range of ice-phase volume fractions. By rigorous, we mean an approach that includes a careful estimation of material properties, is not hampered by the above-mentioned thermodynamic oversimplifications, and uses an accurate model for the volume-fraction effect. The fact that the sugar solutions surrounding the ice crystals are generally nonideal, and for many cases nondilute, poses some difficulty on the estimation of material properties. To alleviate some of the difficulties, we propose an alternative formulation of Ostwald ripening theory in terms of the Maxwell–Stefan[29,30]—instead of the Fick[29,31]—approach to diffusion. The merit of doing this is a cancellation of the thermodynamic factor (describing the nonideality of the solution), which is generally a difficult property to measure or predict. Moreover, we show that Maxwell–Stefan diffusion coefficients can be estimated from self-diffusion coefficients using the Darken equation,[32] which further simplifies the approach in comparison to one where Fick diffusion coefficients have to be estimated.

The various thermodynamic simplifications that have found common use in describing Ostwald ripening have clearly made it difficult to decide on which form of theory to use for a particular system of interest. To clarify this, we find it necessary to rederive some of the basics of Ostwald ripening theory in section (and in greater detail in the Supporting Information of this work). In addition to the traditional LSW approach, various extensions of LSW theory to finite volume fraction are discussed. Discrimination between the various models is made by comparing the volume-fraction-dependent part of the coarsening rate to results obtained from phase-field (PF) simulations or multi-particle diffusion (MPD) simulations from the literature.[18,33,34] In section , we introduce the advocated Maxwell–Stefan approach for describing diffusion and show how this approach leads to a simplified description of Ostwald ripening. Section focuses on the calculation of material properties. In section , theoretical predictions for the coarsening rate in a water–sucrose system are compared to experimental results from the literature.[26,27] Additionally, we calculate the typical time to grow crystals of a certain size and discuss the validity of several commonly used approximations. Our findings are summarized in section .

Theory of Ostwald Ripening

In the following discussion, we limit ourselves to curvature-driven, bulk-mass-transfer limited coarsening of spherical crystals in a liquid matrix. The matrix phase is assumed to be isothermal and of uniform pressure. The surface energy is assumed to be isotropic and independent of crystal size. For a more complete analysis, we refer the reader to Ratke and Voorhees[9] or the review article by Voorhees.[35]

System Definition

We assume a system comprising a large number of spherical ice particles (phase β) dispersed within a bulk matrix phase (α) of temperature T and uniform pressure pα. The particles are pure in water (species 1), while the matrix phase is a binary mixture of water and sucrose (species 2). The composition of the matrix phase is defined by the mole fraction xα = (x1α, x2α). The phases are separated by a sharp interface s (of negligible thickness in comparison to the crystal radii), at which any densities (of total mass, species mass, moles, internal energy, etc.) are assumed to change discontinuously. Mass and molar densities are denoted as ρ and c, respectively. The particle phase is considered to be spatially uniform. Given that the system is not in a global equilibrium state, the concentration profiles within the matrix phase are nonuniform and depend on the spatial coordinate vector r = (x,y,z).

The distribution of crystals within the system is described by a local particle density n(r,a,t), which, under multiplication with dadr, defines the number of crystals of radius between a and a + da in a volume element between r and r + dr. Assuming that the crystals are distributed uniformly throughout the system, the particle density reduces to n(r, a, t) = n(a, t), which can be factorized further to the uniform number density n(t) = N(t)/V (with N being the number of crystals and V the volume of the system) and a particle-size distribution (PSD) f(a,t) as n(a,t) = n(t)f(a,t). The PSD is a probability distribution and normalizes to unity.

General Considerations

During coarsening, the flux of crystals from size class a to a + da is given by n(a,t)⟨ȧ⟩, where ȧ ≡ da/dt is the radial growth rate of a single crystal of radius a, and ⟨ȧ⟩ denotes a statistical average over all crystals of this size. Accordingly, the time evolution of the particle density must fulfill the continuity equationwhere the right-hand side is zero due to the absence of particle nucleation and accretion.

The microscopic growth rate is related to the material fluxes relative to the moving interface by mass (or molar) balance. Due to the imposed spherical shape of the crystals, only the radial component of these fluxes needs to be considered. Considering a flux as positive when directed along the radial coordinate r starting from the center of mass of a particle, we derive (Appendix A)where Jα(a) is the diffusional part of the radial flux of species i, evaluated at the matrix side of the interface. If both phase densities are equal, the denominator reduces to the difference in species density between the two phases, and eq takes on the form generally found in the literature on particle coarsening.[9,36] The form presented here explicitly accounts for differences in phase densities, which for ice in concentrated sucrose solutions can be non-negligible.

For binary systems, the species mole fractions at the matrix side of the particle surface xα(a) are fully determined by imposing the condition of local thermodynamic equilibrium. Due to the effects of curvature, the conditions for equilibrium are slightly different in comparison to those for two bulk phases. As shown in Appendix B, the difference can be captured accurately using the following analytical approximation—usually referred to as the Gibbs–Thomson equationwhere the subscript eq was introduced to denote the value of a property at bulk phase equilibrium. The parameter lc, is a capillary length scale, calculated aswhere σ is the ice matrix surface energy, Vmβ is the molar volume of ice, Δx = xβ – xα is the composition difference between ice and solution (i.e., miscibility gap), R is the gas constant, and Γ is the thermodynamic factor of the matrix phasewhich, for a binary mixture, is independent of the index i. For an ideal solution, the thermodynamic factor reduces to unity. We note the sign of the capillary length as defined by eq depends on whether i is chosen as the precipitating species (Δx > 0, l > 0) or nonprecipitating species (Δx < 0, l < 0).

The diffusional flux is usually calculated on the basis of Fick’s law. Although (in section ) we will show that Ostwald ripening is actually more conveniently described in terms of the Maxwell–Stefan approach to diffusion, for the purpose of review we will for now stick to the notation used in the literature. Defining the molar-averaged velocity of the mixture as the reference velocity to which diffusional fluxes are calculated, Fick’s law is written as[29,31]where D is the binary Fick diffusion coefficient, which is the same for any of the two species in the mixture.

A combination of eqs , 4, and 7 allows the microscopic growth rate of a crystal to be expressed in terms of its surrounding composition profile. Assuming the curvature-induced shift in equilibrium mole fractions is negligible in comparison to the miscibility gap (Δx ≫ (2l/a)xα), one obtainsAlthough it might seem the effects of curvature are hereby neglected altogether, it is important to realize that the curvature-induced shift in equilibrium composition (eq ) still determines the boundary condition for the composition profile at the crystal surface and thereby the microscopic growth rate.

The main difficulty for any theoretical approach aimed at describing Ostwald ripening is to define the second boundary condition: namely, the composition at some position in the matrix. The reason is that this composition depends on the local environment of a crystal, which is generally not known. In the majority of theoretical work, the local environment is treated approximately, by using some statistically averaged composition profile for each crystal size class. Different theories vary, however, in how the averaging is performed.

Limit of Zero Volume Fraction: LSW Theory

In the analysis of LSW, all particles are assumed to interact with an infinite medium of mean-field composition xα = (1/Vα)∫xα(r)dr at r → ∞. Effectively, this implies a vanishing volume fraction of dispersed-phase material.

The LSW theory also assumes quasi-static conditions, negligible convective flow, and composition-independent diffusion coefficients and total density, so that the continuity equation reduces to the Laplace equation, leading toThe gradient at the surface of a particle follows asInserting eq for the mole fraction at the surface, and introducing the supersaturation θα = (x – xα)/xα, this is rewritten aswhere the term in brackets can be interpreted as a local supersaturation at the (curved) surface of a particle. This result suggests that, for any distribution of particles undergoing Ostwald ripening, there will be a critical size class with particles of radius a* = 2lc,/θα for which the averaged growth rate equals zero. Introducing the critical radius in eq and substituting the result in eq leads to the averaged growth rate of particles within size class a, aswith the dimensional prefactorRewriting to the averaged volumetric growth rate B ≡ ⟨da3/dt⟩ = 3a2⟨ȧ⟩ allows eq to be expressed solely in terms of the rescaled particle radius z = a/a* asThis result clearly shows that particles with size greater than the critical particle radius (z > 1) grow, while particles with smaller size (z < 1) shrink. Thereby the main characteristic of the Ostwald ripening process is captured.

By combining the growth rate (eq ) with a mass balance constraint for the precipitating species, LSW deduced the existence of an asymptotic solution to the continuity equation (eq ) for long coarsening times, according to n(z,t) = feq(z)n(t). The self-similar solution for the rescaled PSD, here denoted as feq(z), could be calculated analytically. In the Supporting Information it is shown that the PSD takes the formWithin the asymptotic regime, the temporal evolution of quantities such as the average particle size, the critical radius (or supersaturation), and the number of particles was shown to obey several scaling laws, the proportionality constants of which, called coarsening rates, follow naturally from the theory. The coarsening rate for the scaling of the average particle size cubed with time (eq ) is obtained aswith ξ the dimensional prefactor of eqs and 14. Please see the Supporting Information for a detailed derivation.

Extensions to Finite Volume Fraction

The extension of LSW theory to finite volume fraction has a rather long history,[9,12−22,35] and to review all of its aspects is outside the scope of this paper. We here focus only on the key aspects and results; for a detailed review, the reader is referred to the Supporting Information.

Basically, there are two main techniques used to connect the averaged growth rate to the volume fraction of dispersed-phase particles: (1) cell models,[12,13,15,18,22] which assume each particle to grow within a cell of finite size (connected to volume fraction) instead of an infinite medium, and (2) source-sink models,[14,16,17,19−21,35] which assume each particle to interact with an ensemble of sources and sinks of the precipitating species, the strength of which is connected to volume fraction through rigorous cluster expansion methods[16,17,19−21] or more empirical techniques connecting to macroscopic constraints of the system of interest.[14] The growth and coarsening rate obtained from any of these models can generally be written as a simple extension of the LSW result, with the additional dependence on volume fraction included as a multiplicative factor. For the coarsening rate, one obtainswith ξ being the dimensional prefactor of eqs and 14.

In this work we analyze three different models for calculating the volume-fraction-dependent factor g(ϕ): the cell model of Marsh and Glicksman (MG)[22] and the source-sink models of Brailsford and Wynblatt (BW)[14] and Marqusee and Ross (MR).[17] All models require numerical evalution of g(ϕ), as the growth rates generally carry an implicit dependence on the PSD and vice versa (see the Supporting Information for details). We find the numerically calculated values are well represented by the following correlations in dispersed-phase volume fractionwhich will be used throughout this work.

Model Discrimination

In Figure , the volume-fraction dependent part of the coarsening rate, as calculated using eqs –20, is compared to results obtained from multi-particle-diffusion (MPD) simulations[18,33] and phase-field (PF) simulations[34] from the literature. In contrast to the theories leading to eqs –20, these simulation methods rigorously solve for the composition profile within the matrix (albeit in the Laplacian approximation eq ) and thus lead to a more complete description of the local environment of the crystals.[9] These simulations thus provide a valuable reference for model discrimination. We find that the MG theory results in the closest representation of the effect of volume fraction on the coarsening rate; therefore, this will be our model of choice. For volume fractions lower than ϕ = 0.6, the BW theory also leads to reasonable comparison to simulations. The MR theory is only accurate at low volume fractions.

Figure 1

Volume-fraction-dependent part of the coarsening rate g(ϕ) = K(ϕ)/KLSW as obtained from various theories (MG, Marsh–Glicksman;[22] BW, Brailsford–Wynblatt;[14] MR, Marqusee–Ross[17]) in comparison to results obtained from phase-field simulations of Kim[34] and multi-particle-diffusion (MPD) simulations of Voorhees–Glicksman[18] and Akaiwa–Voorhees.[33] The red line is a correlation to simulation results as obtained from eq .

Simplification in Terms of Maxwell–Stefan Diffusion Coefficients

To predict coarsening rates in real systems, one requires accurate measurement or prediction of the material properties in eqs and 13. Given that the required material properties involve surface properties (surface energy), dynamic properties (binary Fick diffusion coefficient), and properties describing complicated solution thermodynamics (thermodynamic factor), this is generally a very difficult task, and any simplifications leading to cancellation of some of these properties are highly desirable. We find a considerable simplification is possible by introducing the Maxwell–Stefan (MS) approach for diffusion.

The MS approach assumes a gradient in chemical potential (divided by temperature) as the driving force for isothermal–isobaric diffusion (in accordance with nonequilibrium thermodynamics[37]), leading to the following equation for the diffusional flux with respect to a molar-averaged reference frame[29,30]where Đ is the MS diffusion coefficient. As both frameworks describe the same diffusional flux, inspection of eqs and 21 shows that Fick and Maxwell–Stefan diffusion coefficients are related, viawhere Γ is the thermodynamic factor as introduced in eq . The usefulness of introducing the MS approach in this context lies largely in the fact that insertion of eq in the dimensional prefactor of eq leads to cancellation of the thermodynamic factor (since Γ/Γeq ≈ 1), leading towhereis the capillary length for coarsening in an ideal, nondilute solution (for which by definition Γ = 1).

With this simplification (which can be considered exact), the problem of estimating D and Γ has shifted to estimating Đ. As we show in the next section, this can be done merely on the basis of knowledge of the self-diffusion coefficients within the mixture, which constitutes a second advantage of using the MS approach.

Estimation of Material Properties

Bulk Properties

The density of ice is approximated by its respective value at T = 273.15 K as ρβ = 916.7 kg/m3.[38] The density of the matrix is approximated by that of an ideal mixture, according towhere wsα is the mass fraction of sucrose in the matrix and vw* and vs* are the specific volumes (per unit mass) of pure (as denoted by the asterisk) water and sucrose, respectively. The latter are approximated by their values at T = 273.15 K[38,39] as vw* = 999.8425 kg/m3 and vs* = 1592 kg/m3. Calculated densities of the solution were compared to experimental values at T = 293.15 K;[38] for the entire range of ws percentage relative deviations are within 0.8%. For the systems analyzed in this paper, deviations are expected to be even smaller.

The mass fraction of sucrose at bulk-phase equilibrium is obtained by correlating available experimental data for the freezing-point depression of ice in aqueous sucrose solutions[38,40,41] aswith a0 = 1.195 and a1 = 5.781 and Tm being the melting temperature in degrees Celsius. In principle, the experimental data used for this correlation include melting temperatures down to Tm = −30 °C. However, independent fits of ws,eqα and Tm showed consistent results for a0 and a1 only when data for Tm < −14 °C were discarded. The data points at these low temperatures therefore seem unreliable, and we chose not to incorporate them in the correlation. Equation can thus be considered accurate for Tm > – 14 °C and ws,eqα < 0.62.

Surface Energy

The surface energy of ice immersed in a liquid is a difficult property to measure experimentally. For ice in pure water some consensus seems to have been reached, however, that, for a flat interface (small curvature) at T0 = 273.15 K, the orientationally averaged surface energy σ*(T0) ≈ 29 mJ/m2.[42,43] The agreement of this value with results obtained from molecular simulations based on state of the art force fields for water[44−47] puts some confidence in its use. As molecular simulation results[46,47] indicate the difference in surface energy between the different ice planes is only marginal, it seems justified to use only the orientationally averaged value. Available experimental and simulation work seems to agree that the surface energy linearly increases with temperature, with the slope ∂σ*/∂T ≈ 0.2 mJ/(m2 K).[42] The surface energy of ice in pure water may thus be approximated aswith T either in degrees Celsius or Kelvin.

For ice in aqueous solutions of sucrose, no reliable experimental or simulation data for the surface energy is available. Here we derive an approximate result by borrowing a concept from a paper by Warkentin et al.,[48] in which the experimentally observed reduction of the homogeneous nucleation temperature of ice due to addition of solutes is explained in terms of classical nucleation theory (CNT). Warkentin et al. argue that, for a critical nucleus to form, one first needs to clear a region of solution from solutes that is large enough to contain the nucleus. The probability for this to happen is proportional to exp(−ΠV/kBT), where Π = −(RT/Vwα) ln aw is the osmotic pressure of the solution, aw and Vwα are the activity and partial molar volume of water in the solution, respectively, and V is an exclusion zone of sizewhere a and as are the radius of the nucleate and the Stokes radius of a solute molecule, respectively. The ratio of the molar volumes of pure water and ice Vw*/Vmβ appears because a smaller volume of water needs to be cleared from solutes in comparison to the volume that the nucleus will occupy. We here approximate the partial molar volume of water in solution by the molar volume of pure water Vwα ≈ Vw*, which for solutions of sucrose is a very accurate approximation. In the framework of CNT, the Gibbs free energy of nucleation becomeswhere Δv = cβ(μwβ(T) – μwα(T)) is the Gibbs free energy difference per volume of ice formed from pure water at a temperature T. The mechanical work associated with the expansion of water when turned into ice is proportional to the atmospheric pressure and is neglected. Expanding the product (a + as)3 and collecting the same powers of a, we obtainwhich suggest the surface energy of ice in an aqueous solution of sucrose can be approximated asCombination with eqs and 27 allows estimation of the surface energy of ice along the saturation line, for which we show some results in Figure . These results were calculated on the basis of an estimate of the Stokes radius of sucrose as = 4.4 Å (obtained from the molar density of pure sucrose) and the correlation of the water activity of Zobrist et al.[28,49] for calculating the osmotic pressure. A near-linear dependence on the melting temperature is observed, which is well correlated asWe predict a small but significant effect of sucrose on the surface energy of ice, with an increase of around 10% for a 50 mass % solution.

Figure 2

Surface energy of ice along the melting line of ice in an aqueous solution of sucrose σsat, as obtained from eqs and 26 (black, solid line). The red dashed line results from a linear fit in melting temperature Tm, given by eq .

Warkentin et al. showed that the Stokes radii of various solutes (including the disaccharide trehalose) obtained from measured critical cooling rates (which is the required cooling rate to prevent nucleation and form a glassy system instead) and eq compare quite well to experimental data. Accordingly, we expect eqs and 32 to be reasonable approximations as well.

Diffusion Coefficients

The binary Maxwell–Stefan diffusion coefficient Đ can be estimated based on self-diffusion coefficients, using the Darken or Vignes equations, respectively, as[29]Application of eqs and 34 to Ostwald ripening requires self-diffusion coefficients of the species within the mixture, along the melting line of ice. Here we obtain these by extrapolating the Vogel–Fulcher–Tammann (VFT) fits of the experimental data of Girlich et al.[50] to the melting line calculated by eq . We obtainwith Tm in degrees Celcius. On the basis of eq , these correlations should not be used for Tm < −14 °C or ws,eq > 0.62.

A comparison of both correlations to experimental data of Girlich et al. is presented in Figure . For comparison we included the correlation for the self-diffusion coefficient of water by Zobrist et al.[28] We find good agreement.

Figure 3

Self-diffusion coefficients of water (top) and sucrose (bottom) in binary water–sucrose solutions. Solid lines represent Vogel–Fulcher–Tammann (VFT) fits to the experimental data of Girlich et al.[50] (symbols), whereas dotted lines denote self-diffusion coefficients at the melting line of ice, as obtained from the correlation of Zobrist et al.[28] (red) and eqs and 36 (blue). The crosses denote (from top to bottom) ws = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6. The experimental data for the self-diffusion coefficient of water correspond to (from top to bottom) ws = 0.1, 0.2, 0.3, 0.4, 0.7, whereas for sucrose the data correspond to ws = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7.

As experimental Maxwell–Stefan diffusion coefficients of water–sucrose are not available, validation of eqs and 34 by comparison to experiments can only be done indirectly, by first transforming to the Fick diffusion coeffient as D = ĐΓ and then comparing the result to experimentally obtained Fick diffusion coefficients (which are available for T = 298 K).[51] We note that such a comparison should be handled with caution, as both measured Fick diffusion coefficients and calculated thermodynamic factors are usually prone to large uncertainties. In Figure , we compare theoretical predictions to experimental data. Two different activity models[28,49,52] were analyzed for calculating the thermodynamic factor from eq . Self-diffusion coefficients were obtained by correlating the experimental data of Girlich et al.[50] The overall agreement to experimental Fick diffusion coefficients is quite reasonable. The Darken equation proves slightly superior over the Vignes equation; therefore, this will be our model of choice.

Figure 4

Comparison between experimental and predicted Fick diffusion coefficients of binary water–sucrose mixtures of mole fraction x and T = 298 K. Experimental data is from Irani et al.[51] Predicted values were obtained using eq , on the basis of either the Darken (eq ) or Vignes (eq ) equation to calculate the Maxwell–Stefan diffusion coefficient Đ and the water activity models of either Starzak et al.[52] or Zobrist et al.[28] to calculate the thermodynamic factor Γ. The self-diffusion coefficients as required by the Darken or Vignes equation were obtained from a power law in xs that was correlated to the experimental self-diffusion coefficients of Girlich et al.[50]

We note that although the Darken equation is usually put forward as empirical, it can in fact be derived from statistical mechanics, using the relation among Đ, Onsager’s phenomenological coefficients, and the Green–Kubo relation, respectively.[53] The full statistical–mechanical result contains two contributions. The first contribution, which constitutes the Darken equation, includes only the effect of velocity autocorrelations (i.e., correlations between the velocity of a single molecule at different times). The second contribution, which is neglected if the Darken equation is used, contains an additional cross-term that accounts for velocity cross-correlations (i.e., correlations between the velocity of different molecules of same or different type at different times). The cross-term scales as x1(1 – x1) and thus reduces to zero for very dilute or concentrated solutions. Although the relatively accurate predictions obtained using the Darken equation at T = 298 K suggest the cross-term omitted in eq is small, it is uncertain whether this remains so for temperatures as low as the saturation temperature. As no data of binary (Fick or Maxwell–Stefan) diffusion coefficients is available for such low temperatures, and as this type of data is generally difficult to obtain experimentally, it could be of value to try and obtain it from Molecular Dynamics simulations (e.g., along the lines of the work of Batista et al.[54]). In such simulations the cross-term can be calculated directly,[53] leading to a quantitative measure for the applicability of the Darken equation to water–sugar systems. We should add that this requires development of accurate force fields, which is by no means a trivial task for these kinds of systems.[54−56]

Results and Discussion

Comparison to Experiments

In Figure , we compare theoretical calculations of the coarsening rate of ice in aqueous solutions of sucrose to experimental data of Budke et al.[27] and Hagiwara et al.[26] Theoretical calculations are based on the MG theory (eqs , 18, 23, and 24), with material properties obtained from the correlations listed in section , specifically eqs and 26 for the matrix density and saturated mole fraction of sucrose (when converted using molar masses), eq for the surface energy of ice, and eq with eqs and 36 to calculate the binary Maxwell–Stefan diffusion coefficient. Considering the many assumptions underlying the theory, and the difficulties in estimating material properties, the agreement is very reasonable, with relative deviations smaller than a factor of 2 for all cases studied. It is important to note that no adjustable parameters were used; the shown results are predictions.

Figure 5

Comparison of predicted coarsening rates K of ice in aqueous solutions of sucrose to experimental results of Budke et al.[27] and Hagiwara et al.[26] for varying volume fractions of ice ϕ. The data point of Hagiwara in the top figure is for the temperature T = −5.8 °C. Solid lines are results obtained from the MG theory[22] combined with the Darken equation for estimating the MS diffusion coefficient. Dashed lines are obtained by using eq to correct for minor inaccuracies in the description of the volume-fraction effect resulting from the MG theory. Dotted lines additionally involve use of the empirical modification of eq to correct for any inaccuracies in the MS diffusion coefficient resulting from the Darken equation.

Let us first discuss the intermediate- to high-volume fraction regime. For temperatures T = −6, −8 °C theoretical predictions overestimate experimental results, while for the lowest temperature studied T = −10 °C, a slight underestimation of experimental coarsening rates is observed. In principle, one would expect a theory for diffusion-limited Ostwald ripening such as that used here to lead to an underestimation of the coarsening rate, as various effects known to enhance coarsening (e.g., accretion, crystal movement, spatial correlations between crystals)[9] or anticipated to enhance coarsening (e.g., small convective fluxes in the matrix induced by the density difference between the phases) are not explicitly accounted for.

One possible reason for the observed overestimation (instead of the expected underestimation) of coarsening rates at T = −6, −8 °C is the overestimation of the volume-fraction effect due to using the MG theory (see Figure ). To correct for this, we correlated the PF and MPD simulation results from Figure (dashed lines) asAs shown in Figure , however, application of this correction does not fully resolve the issue.

A second reason for the overestimation could be inaccuracies in the Darken equation, caused by neglecting the contribution due to velocity cross-correlations between water and sucrose. In fact, as a good rule of thumb, the Darken equation overestimates the MS diffusion coefficient when the thermodynamic factor is larger than unity, while it underestimates the MS diffusion coefficient when the thermodynamic factor is smaller than unity.[57,58] For the water–sucrose solutions considered here, available experimental data[39,40] and thermodynamic models/correlations[28,49,52,59] all agree that the thermodynamic factor is larger than unity; therefore, use of the Darken equation is expected to lead to an additional overestimation of the coarsening rate. To account for inaccuracies in the Darken equation, Moggridge[57] proposed an empirical modification of the formwhich was shown to compare very well to the experimental Fick diffusion coefficients of a large number of binary mixtures. As the relation between MS and Fick diffusion coefficients of eq is exact, a more appropriate way of writing the above modification is in terms of the MS diffusion coefficient, asApplying this equation to the water–sucrose mixtures under consideration leads to a further reduction of the coarsening rate toward experimental values (dotted lines in Figure ). The required thermodynamic factor was obtained by correlating available experimental data for the water activity of saturated sucrose solutions of Lerici et al.[40] and Norrish (see Mathlouti[39]), aswhich should not be used for ws,eq > 0.7. We note the results based on eqs and 40 should be interpreted merely as an indication of the inaccuracies underlying the Darken equation, as the thermodynamic factor is very sensitive to any experimental inaccuracies in the activity and eq was not tested directly on water–sugar systems.

The small remaining overestimations of the experimental coarsening rate are most probably caused by inaccuracies in the correlation for the self-diffusion coefficient of sucrose along the saturation line (see blue lines in Figure ; ws = 0.46 and ws = 0.52), as this involved quite severe extrapolation of experimental data to lower temperatures.

In summary, when relieved from small inaccuracies introduced by the Darken equation, experimental self-diffusion coefficients, and the approximate description of the volume-fraction effect, theoretical predictions will most probably underestimate experimental coarsening rates in the intermediate- to high-volume-fraction region for all three temperatures studied. This is the expected behavior for a theory of the form used here.

Let us now turn our attention to the region of very low volume fractions, where experimental coarsening rates remain overestimated. We believe the overestimation at low volume fractions to be an artifact of the experimental setup,[26,27] which consisted of a microliter sample kept between two small glass plates. Probably to ensure good temperature control of the sample, the spacing between plates was chosen to be very narrow, on the order of ∼10 μm. A side effect of such a small spacing is that any large crystals in the sample must have touched the plates and grew in an effective 2-D system. As argued by Marqusee,[60] the difference between 2-D and 3-D coarsening is most notable in the low-volume-fraction region (or area fraction for 2-D); therefore, the observed discrepancy between theory and experiments in this region could be due to the fact the system is in fact in a crossover regime between 2-D and 3-D coarsening.

Evolution of Average Crystal Size during Storage

Given the satisfying comparison between theory and experiments as presented in the previous section, we now calculate some results that could be useful as a guide for preserving frozen foods or biological materials. A relevant quantity in this context is the time it takes for crystals to grow to a size large enough to lead to degradation of product quality. As this size may differ from product to product, we calculate the time it takes to multiply the average crystal radius by a factor n, according toThe equilibrium volume fraction of ice ϕ (which is needed to calculate K and thus t) is related to the initial mass fraction of sucrose in the system w0 (where initial means before formation of any ice) and temperature T through the mass fraction ϕm = 1 – w0/ws,eq(T), according toResults for an initial weight fraction of sucrose w0 = 0.28 and initial radius ⟨a⟩0 = 30 μm are shown in Figure for n = 2, 3, 4, 5, 6, as a function of storage temperature. For the largest part, the growth time increases nearly exponentially with decreasing temperature, a result caused by the temperature dependence of the self and MS diffusion coefficients. The observed minimum results from a competition between an increasing volume fraction and decreasing dimensional prefactor (predominantly due to the diffusion coefficient) as the temperature is lowered. For temperatures close to 0 C the effect of volume fraction is stronger, causing the growth time to first slightly decrease.

Figure 6

Time it takes to multiply the average crystal radius by a factor n, for a system of initial average crystal radius ⟨a⟩0 = 30 μm and an initial mass fraction of sucrose (before formation of ice) w0 = 0.28. Solid lines from bottom to top correspond to n = 2, 3, 4, 5, 6.

Given the range of validity of some of the correlations developed in section , we only show results for T > −14 °C. Therefore, the diagram of Figure is probably more appropriate for storage of frozen foods (e.g., in commercial or home freezers) than for the cryostorage of biological tissue (which usually involves lower temperatures).

To allow easy calculation of the growth time t for arbitrary ⟨a⟩0, w0, and n, we correlated the dimensional prefactor of eq (in units of μm3/min) in terms of the melting temperature Tm (in degrees Celsius) aswith a0 = 2.8425, a1 = 0.50262, a2 = 0.034348, a3 = 0.0016574, a4 = 3.3152 × 10–5, and −14 °C < Tm < −1.5 °C. Together with eq , and any of the correlations given in eqs −20, and 37 for the volume-fraction dependence of the coarsening rate, this equation fully determines K, and thus t via eq .

Analysis of Common Approximations

Let us finally turn to analyzing the effect of some commonly made approximations. The validity of any (non-volume-fraction-related) approximations added to the theory can be tested by their effect on the dimensional prefactor ξ (see eqs and 24). For such an analysis it proves convenient to write the prefactor as a product of three contributions, aswithwhere the last term can be simplified in terms of the mole fraction of the precipitating species (water) as

A common approximation,[24,27] which stems from the original analysis of LSW, is the assumption the matrix phase is diluted in the precipitating species (i.e., 1 – x1,eq ≈ 1, D ≈ Đ). While for precipitation of metals or pharmaceuticals from (solid) solution, this is generally a quite reasonable assumption, for concentrated aqueous solutions of sugars it is not. For example, for a 50 mass % sucrose solution, the error introduced in ξ (and thus in the coarsening rate K) is around a factor of 1/10.

Another common approximation is equal phase densities. The error introduced by this (through ξ1) is much smaller, but for the systems considered here, it is still significant. For a 50 mass % sucrose solution, the error introduced is approximately a factor of 1.5.

Finally, in some studies it is assumed[27] (or claimed[26]) that the self-diffusion coefficient of pure water should be used instead of the binary diffusion coefficient to calculate ξ2. For a saturated 50 mass % sucrose solution, the error introduced by this is easily 1 order of magnitude. We stress that this is not merely a reflection of the fact that the solution considered is nondilute. Even for an infinitely diluted solution, use of the self-diffusion coefficient of pure water to approximate the binary Maxwell–Stefan diffusion coefficient is simply incorrect. In fact, in that limit, the binary Maxwell–Stefan coefficient reduces to the self-diffusion coefficient of sucrose infinitely diluted in water, as indicated by eq (which at infinite dilution becomes exact). What has probably added to the confusion is the fact that, when the self-diffusion coefficient of pure water is used in combination with the incorrect assumption of a matrix diluted in the precipitating species, many of the errors cancel and the obtained coarsening rates compare quite well to experimental results (see e.g. the results for zero volume fraction in the work of Budke et al.[27])

We conclude that a predictive description of coarsening of ice in sugar solutions requires the full dimensional prefactor from eqs and 5 or eqs and 24 to be considered. In principle, all material properties in this prefactor need to be estimated carefully. Although the prefactor is most sensitive to errors in the bulk-equilibrium value of the mole fraction (see eq ), this property can usually be obtained accurately without much effort. The largest problem lies in accurate estimation of the binary diffusion coefficient.

Summary and Conclusion

We examined whether a theory for bulk diffusion-limited Ostwald ripening can act as a predictive tool for estimating isothermal recrystallization rates of ice in a model system for frozen foods or biomaterials. For a binary water–sucrose system, good predictions of recrystallization rates of ice are obtained, with relative deviations to experimental values not exceeding a factor of 2.

To obtain predictions of this quality, it proved necessary to rectify several misconceptions that seem to have entered the literature on this subject. A detailed discussion on this was given in section . A first conclusion that can be drawn from this is that in order to describe coarsening in nonideal, nondilute solutions of nonequal phase densities, one should use the dimensional prefactor as derived in this work (eqs and 13 or eqs and 24). The simplified form as presented in the work of LSW (and in the majority of papers citing LSW) is simply inadequate for describing such systems. Second, we should stress that the binary diffusion coefficient cannot be replaced by the self-diffusion coefficient of pure water, even if the solution is dilute. There are various examples in the literature where such an approach is assumed—or even claimed—to be correct, resulting in large errors in predicted coarsening rates.

To facilitate estimation of material properties, we introduced the concept of Maxwell–Stefan diffusion coefficients. On the basis of this concept, an alternative, but exact, formulation of Ostwald ripening theory was proposed, which does not require calculation of thermodynamic factors. Additionally, using the Darken equation to estimate binary Maxwell–Stefan diffusion coefficients, an approximate reformulation was proposed, which requires merely self-diffusion coefficients of the species in solution as input. This greatly simplifies the approach.

Our calculations indicate that the characteristic time for crystal growth by Ostwald ripening is rather long: e.g., for a system of initial (before formation of ice) mass fraction of sugar w0 = 0.28 and initial average crystal radius of ⟨a⟩0 = 30 μm, on storage at −12 °C, it takes approximately 2 years to grow to an average radius of 90 μm. On the basis of this insight, we conclude that the relatively large crystals sometimes found in frozen foods such as ice cream might be due to other driving forces for recrystallization not taken into account by the model presented in this paper. We expect temperature fluctuations to play a role here.[61−64] It might be that our assumption of a truly isothermal system is insufficient, as freezer temperatures typically fluctuate to some extent, due to either opening–closing of the freezer or periodic freeze–thaw cycling to thaw off ice from freezer walls. The effect of temperature fluctuations on the recrystallization rate will be the subject of a subsequent work.