## Abstract

We assess the morphological stability of a non-equilibrium ice–colloidal suspension interface, and apply the theory to bentonite clay. An experimentally convenient scaling is employed that takes advantage of the vanishing segregation coefficient at low freezing velocities, and when anisotropic kinetic effects are included, the interface is shown to be unstable to travelling waves. The potential for travelling-wave modes reveals a possible mechanism for the polygonal and spiral ice lenses observed in frozen clays. A weakly nonlinear analysis yields a long-wave evolution equation for the interface shape containing a new parameter related to the highly nonlinear liquidus curve in colloidal systems. We discuss the implications of these results for the frost susceptibility of soils and the fabrication of microtailored porous materials.

## 1. Introduction

When a colloidal suspension such as a clay soil freezes, an interesting phenomenon occurs in which the ice segregates into lenticular or dendritic shapes. For example, in soils composed mainly of silt particles, the ice usually forms nearly planar layers, called *ice lenses* (Taber 1929; Watanabe & Mizoguchi 2000), whereas in a wide range of clays and colloidal suspensions, the ice can grow in several directions, producing a remarkable variety of patterns (Taber 1929; Chamberlain & Gow 1979; Brown 1984; ; figure 1). Recently, we proposed a mechanism for ice segregation in freezing colloidal suspensions based on the concept of morphological instability studied in alloys quantitatively applied to the ice–suspension interface (Peppin *et al.*2006, 2007, 2008). The approach relies essentially on treating a colloidal suspension as a ‘two-component thermodynamic system’ in which the solute particles are vastly larger than the solvent, and by exploiting the physical basis for the dynamics and thermodynamics of pre-melting between these particles and the frozen solvent. Here, we extend this approach by incorporating the anisotropic growth kinetics of ice and by developing a weakly nonlinear analysis of long-wavelength modes to reveal evolution equations derived in the alloy setting and others.

Our previous experimental studies focused on using bentonite clay as an easily accessible and broadly relevant colloid and hence we begin in §2 by discussing its physico-chemical properties. The freezing-point depression and particle diffusivity *D*(*ϕ*) are given as functions of particle volume fraction *ϕ*. It is found that over a significant range in *ϕ* the diffusivity is approximately constant, simplifying many aspects of the stability analysis to follow. The interfacial properties of ice are discussed in §3, where we consider the effect of interfacial curvature on the equilibrium freezing temperature, and allow for the non-equilibrium effects of kinetic undercooling and kinetic anisotropy. Moreover, we discuss the mechanism of particle trapping and its basis for the segregation coefficient. At low freezing velocities, the segregation coefficient is zero, allowing an experimentally convenient scaling to be employed in our analysis of colloidal suspensions. Owing to the fact that our analysis of the thermophysical properties of bentonite suspensions allows us to construct interfacial and boundary conditions that are directly analogous to the binary alloy system, we then, in §4, construct a model of the unidirectional solidification of such suspensions. Indeed, it is found that the clay becomes ‘constitutionally supercooled’ under typical freezing conditions, in agreement with experimental results (Brown 1984; Peppin *et al*. 2008). We then analyse the morphological stability of a planar steady-state ice–clay interface in §5 and obtain an equation for the critical parameters at which the interface becomes unstable. The linear analysis predicts the possibility of travelling-wave instabilities, as found in the atomic alloy system by Young *et al*. (1987), and we explicitly find their operative range of wavenumbers in terms of the model parameters. The structure of the linear problem suggested a long-wavelength weakly nonlinear analysis would be fruitful. As such, we derive an evolution equation for the interface shape in this limit. In §6, we discuss our results, describe the conditions that must be satisfied for the stable solidification of a colloidal suspension and suggest future research directions and applications to geophysical and technological systems.

## 2. Physico-chemical properties

### (a) Freezing-point depression

In this section we discuss the freezing-point depression of bentonite, a material for which an extensive range of experimental data is available. Figure 2 shows a system in which a macroscopic portion of ice is in equilibrium with a macroscopic quantity of unfrozen water-saturated bentonite clay. As the temperature is lowered, the ice phase grows and the particle volume fraction of the clay increases. At a critical *breakthrough* temperature *T*_{B} and volume fraction *ϕ*_{B}, ice is able to enter the pore space between particles (Brown & Payne 1990; Shanti *et al.* 2006). For kaolinite clay, the ice-entry temperature is *T*_{B}=−0.77 °C (Brown & Payne 1990). As far as we are aware, the ice-entry temperature has not been measured for bentonite, though Brown & Payne (1990) found no pore ice down to at least −8 °C. We note that here *T*_{B} is considered an equilibrium temperature. In this sense, it is analogous to the eutectic transformation in alloys. Ice can also enter the pore space at temperatures warmer than *T*_{B} via non-equilibrium particle trapping (Wettlaufer & Worster 2006). We discuss particle trapping further in §3*c*.

Here, we study relatively slow freezing velocities and temperatures above *T*_{B}, such that no pore ice is present. We also assume the ice phase is devoid of particles, in which case the freezing-point depression (colloidal liquidus curve) can be obtained by equating the chemical potential of water in the unfrozen colloid (known, e.g., from osmotic pressure measurements) with the chemical potential of pure ice. For temperatures near to the bulk freezing temperature *T*_{m} of pure water, and neglecting the temperature dependence of the osmotic pressure *Π*, this leads to the relation (e.g. Low *et al.* 1968)
2.1
where *ρ*_{w} and *L*_{f} are the density and latent heat of fusion, respectively, of water. Experimental measurements of the osmotic pressure of bentonite can be fit well to an equation of the form
2.2
where *v*_{p}=(4/3)*ΠR*^{3} is the volume of a particle of radius *R* (here 0.5 μm), *k*_{B} is Boltzmann’s constant and
2.3
is the compressibility factor, where *a*_{1}=1×10^{7}, *a*_{2}=2×10^{9}, *a*_{3}=3×10^{9}, *a*_{4}=−8×10^{9} and *ϕ*_{p}=0.64 are fitting parameters (Peppin *et al.* 2008).

### (b) Permeability and diffusivity

To obtain the particle diffusivity, we use a generalized form of the Stokes–Einstein relation
2.4
where *D* is the particle diffusivity, *k* is the permeability, *η* is the fluid viscosity, *T* is the absolute temperature and *P* is the mixture pressure (Russel *et al*. 1989; Peppin *et al*. 2006). Given measurements of the permeability and osmotic pressure *Π*, equation (2.4) can be used to determine the diffusivity as a function of volume fraction. Experimental measurements of the permeability of bentonite can be fit to the equation (Peppin *et al.* 2008)
2.5

Hence, with equations (2.2), (2.3) and (2.5), and using a value of *η*=0.002 kg m^{−1} s^{−1} for the viscosity of water at 0°C (Hallett 1963), equation (2.4) yields the diffusivity of bentonite as a function of *ϕ*, plotted in figure 3. Interestingly, the diffusivity remains approximately constant over a large range in volume fraction before finally diverging near the random close packing limit *ϕ*_{p}. This property of the diffusivity will enable us to simplify the stability analysis in §5.

## 3. Surface properties of ice

### (a) Interfacial curvature

In this section, we discuss the effect of curvature on the ice–colloid interface temperature in the case where the macroscopic radius of curvature of the interface is much larger than the hydrodynamic radius of a particle. The thermodynamic conditions for equilibrium are then given by (e.g. Defay & Prigogine 1966)
3.1
3.2
3.3
where *T*_{i} is the temperature of the ice, *T*_{c} is the temperature of the colloid; *T*_{Ie} is the equilibrium interface temperature, is the chemical potential of the ice, is the chemical potential of water in the colloid, *P*_{i} is the pressure of the ice, *P*_{c} is the pressure of the colloid, *γ* is the ice–colloid surface energy and *K* is twice the mean interfacial curvature (defined to be negative for a convex finger of ice extending into the colloid phase). Not much is known about the concentration dependence of *γ*, though limited experimental evidence indicates that the dependence is weak, or at least bounded as the concentration increases (Schramm & Helper 1994; Dong & Johnson 2004). In the following, we take *γ* to be a constant.

We assume that the densities and specific entropies (entropy per unit mass) of ice and water and that the osmotic pressure *Π*(*ϕ*) of the clay are independent of temperature and pressure over the range of our study, and that there are no anisotropic stresses on the ice. Whence, the chemical potentials can be expanded about the planar ice–water freezing temperature *T*_{m} and atmospheric pressure *P*_{atm} to yield (e.g. Prigogine & Defay 1954; Acker *et al.* 2001)
3.4
and
3.5
where ν_{i}=1/*ρ*_{i}, ν_{w}=1/*ρ*_{w} are the specific volumes of ice and water, respectively, *s*_{i}, *s*_{w} are the specific entropies of ice and water and is the chemical potential of pure ice/water at *T*_{m} and *P*_{atm}. Inserting equations (3.4), (3.5) and (3.2) into equation (3.3) yields
3.6
In obtaining equation (3.6), we used *T*_{m}(*s*_{w}−*s*_{i})=*L*_{f} and we have ignored a term (*P*_{c}−*P*_{atm})(*ν*_{i}−*ν*_{w})/*L*_{f}, which represents the effect of pressure melting on the interfacial temperature and can usually be neglected relative to the other terms (see equation (2) of Dash *et al.* (2006)). For a planar interface *K*=0, and equation (3.6) reduces to equation (2.1).

### (b) Pre-melting, kinetic undercooling and kinetic anisotropy

When the interface between ice and a foreign particle is studied on the nanometre scale, it is found that the ice is not in direct contact with the particle, but rather a thin film of pre-melted water exists between the two surfaces (Dash *et al.* 2006). During solidification, fluid flows around the particle and into the film, enabling the ice to grow. Viscous dissipation in the film leads to a kinetic undercooling (Wettlaufer & Worster 2006), which can be represented on the macroscopic scale by
3.7
where *V*_{n} is the ice velocity in a direction normal to the freezing interface and Δ*T*=*T*_{Ie}−*T*_{I} is the undercooling. The non-equilibrium interface temperature is *T*_{I}, and *T*_{Ie} is the equilibrium (*V*_{n}=0) temperature given by equation (3.6). If the surface properties are anisotropic, as is the case for ice (e.g. Nagashima & Furukawa 1997), then some orientations of the interface will grow faster than others, and we must write
3.8
where *h*(*x*) is the height of the curved interface above the planar location and the subscript *x* denotes differentiation with respect to *x*. In the simplest case (small *V*_{n}), equation (3.8) can be written in a linear form as
3.9
where *κ*_{u} and *κ*_{a} are material-dependent constants accounting for kinetic undercooling and kinetic anisotropy. Solving equation (3.9) for *T*_{I} gives, with equations (3.6) and (2.1),
3.10
where *ϕ*_{I} is the particle volume fraction at *h*. Equation (3.10) gives the ice–colloid interface temperature when the system is not at equilibrium and the interface is curved and moving. In terms of *h*(*x*), the interfacial curvature can be written as
3.11

There is little experimental data on *κ*_{u} or *κ*_{a} for colloidal particles. For 1.1 μm silica microspheres, Watanabe (2002) has obtained *κ*_{u}=0.43 μm s^{−1} °C^{−1}. For bentonite and kaolinite, no measurements of *κ*_{u} or *κ*_{a} are available. Therefore, in the following, we use *κ*_{u}=0.43 μm s^{−1} °C^{−1} and study the effect of a range of *κ*_{a} values.

### (c) Particle trapping and the segregation coefficient

The segregation coefficient *k*_{ s} is defined as the ratio of particle concentration in the ice phase to the particle concentration in the clay phase at the ice–clay interface. In obtaining equations (2.1) and (3.6) it has been assumed that *k*_{ s}=0 (i.e. the ice is assumed pure). This is an excellent approximation when and *T*>*T*_{B}. When *V* ≠0, many theoretical and experimental studies have shown that, above a critical freezing velocity *V*_{c}, the particles are engulfed by the ice (Uhlmann *et al.* 1964; Cisse & Bolling 1971; Rempel & Worster 1999). This suggests that for *V* ≫*V*_{c}. Experimental data indicate that for submicron hydrophilic particles, *V*_{c}≫10 μm s^{−1} (Cisse & Bolling 1971). In the present work, we consider freezing velocities much less than *V*_{c}, so that the particles are completely rejected by the ice and *k*_{ s}=0.

## 4. Unidirectional solidification

### (a) Formulation

Figure 4 illustrates the system to be studied. A layer of colloid of initial volume fraction *ϕ*_{0} is placed to thickness *L*_{0} in a glass cell between two fixed temperature blocks that produce a temperature gradient *G*_{T}, with *T*_{H}>*T*_{C}, where *T*_{C} is below the liquidus temperature *T*_{f}(*ϕ*_{0}) of the bulk suspension. The cell is moved through the blocks at a fixed speed *V* (≪*V*_{c}). All of the particles will therefore be pushed ahead by the growing ice, and if the interface remains stable, a steady state will eventually develop. In the following, we determine a relationship among the parameters *V* , *ϕ*_{0}, *L*_{0} and *G*_{T} that separates the stable and unstable states.

We assume that the sample is contained within thick transparent plates with high thermal conductivity, such that heat flow is mainly through the plates and latent heat effects can be neglected; the so-called frozen temperature approximation (Caroli *et al.* 1982; Davis 2001). In such a case, the cooling blocks induce a constant linear temperature profile
4.1
where *T*_{i} is the temperature of the stable planar interface. From equation (3.10), *T*_{i} can be written as
4.2
where *ϕ*_{i} is the particle concentration at the planar (*z*=0) interface. The first term in equation (4.2) accounts for the equilibrium freezing-point depression while the second term is the non-equilibrium kinetic undercooling.

For simplicity, we assume the clay density *ρ*_{c} is a constant independent of *T* and *ϕ*, though we allow *ρ*_{c} to differ from the ice density *ρ*_{i}. In a frame of reference moving at the pulling speed *V* , the equation describing the conservation of colloid mass is
4.3
where *r*=*ρ*_{i}/*ρ*_{c} and the subscripts *t* and *z* denote differentiation with respect to time and the coordinate parallel to the growth direction. Conservation of mass at the interface *z*=*h*(*x*,*t*) is given by
4.4
where *ϕ*_{I} is the particle concentration at *z*=*h*,
4.5
is the speed of the deforming interface in the normal direction and
4.6
is the unit normal to the interface. The far-field boundary condition on equation (4.3) is
4.7
The temperature at the perturbed interface is given by equation (3.10).

### (b) Steady state

There exists a steady-state solution of the governing equations with a planar interface *h*=0. The temperature field is given by equations (4.1) and (4.2). The steady-state concentration profile can be obtained by integrating equation (4.3) from 0 to *z* using equation (4.4) as a boundary condition to yield
4.8
A boundary condition for equation (4.8) at *z*=0 can be obtained by using global mass conservation in the form
4.9
where *ϕ*_{i} is the volume fraction at the planar interface. As noted in §2*b*, over a large range in *ϕ* the diffusivity is approximately constant, in which case equation (4.9) can be solved to give *ϕ*_{i}=*rϕ*_{0}*L*_{0}*V*/*D*, and equation (4.8) then yields
4.10

Figure 5 shows predicted profiles of volume fraction, temperature and freezing temperature for a bentonite clay having initial volume fraction *ϕ*_{0}=0.1, density ratio *r*=0.9, freezing velocity *V* =0.05 μm s^{−1} in a temperature gradient *G*_{T}=1 °C mm^{−1} and three different initial layer thicknesses. The upper and lower solid lines in each figure show the concentration and temperature profiles, respectively, while the dash-dotted line is the freezing temperature *T*_{f}(*ϕ*). As the initial layer thickness increases, the steep concentration gradient leads to a situation where, for a certain distance ahead of the freezing interface, the temperature is below the freezing temperature, and the suspension is thermodynamically (constitutionally) supercooled. The concentration profile predicted by using the constant diffusivity *D*=6.4×10^{−11} m^{2} s^{−1} from figure 3 is shown by the dotted curves in figure 5. When the volume fraction everywhere in the system is lower than 0.5, the comparison is very good. In the stability analysis to follow, we treat cases where *ϕ* is less than 0.5, so that *D* can be approximated as a constant.

## 5. Stability analysis

### (a) Linear stability

Now we carry out a linear stability analysis of our governing system of equations (4.3)–(4.7) that admit the steady-state solution (4.10) with a planar interface *h*=0. We scale all lengths in the system with δ_{ϕ}=*L*_{0}*ϕ*_{0}, velocities with *V* , time with δ_{t}=*L*_{0}*ϕ*_{0}/*V* and concentrations with δ_{ϕ}*G*_{ϕ}=−*r*δ_{ϕ}*ϕ*_{i}*V*/*D*, where . We define dimensionless variables, denoted by the superscript *, as follows:
5.1
There are five associated non-dimensional groups in the system: the Peclet number
5.2
the morphological number
5.3
where *T*_{fϕ}=d*T*_{f}/d*ϕ*|_{ϕ=ϕi} and *G*_{T} is the imposed thermal gradient; *Γ*, which measures the effects of surface energy,
5.4
a kinetic undercooling parameter given by
5.5
and the kinetic anisotropy parameter defined by
5.6

When we substitute equation (5.1) into the system (4.3)–(4.7), the dimensionless governing equation becomes
5.7
where for convenience we have dropped the stars from the rescaled variables. The equation is subject to the dimensionless boundary conditions
5.8
5.9
where the nonlinear terms in *h*_{x} have been neglected.

We define the dimensionless quantity
5.10
which is proportional to the difference in the freezing-point depression at the perturbed interface *z*=*h*, and that at the unperturbed interface *z*=0. If we recall from equation (4.1) that the non-equilibrium temperature of the perturbed interface, *T*_{I}, is related to *T*_{f}(*ϕ*_{i}) by
5.11
then it is seen that after substituting equations (5.10) and (5.11) into equation (3.10), we obtain the following equation that must be satisfied at the perturbed interface:
5.12
The scaled basic state of the system (5.7)–(5.9) is given by
5.13

We impose small perturbations of the basic state as
5.14
5.15
Inserting equations (5.14) and (5.15) into the governing equations, linearizing in the perturbed quantities, Taylor expanding the boundary conditions about *z*=0 and seeking normal mode solutions of the form
5.16
where α is real and σ=σ_{R}+*i*σ_{I}, leads to the following characteristic equation for the growth rate σ
5.17
where . Since *α** is a complex number, we can explicitly compute the corresponding real and imaginary parts leading to the following coupled system of equations for the growth rate σ:
5.18
where
5.19

Assuming exchange of stabilities, the marginal stability curves are given by σ_{R}=0. From equation (5.19), this gives
5.20
One can analyse equation (5.20) explicitly in one simple case given by
Then *r*Pe−Re(*α**)<0 so that the condition σ_{R}=0 is equivalent to
5.21
and equation (5.21) is the equation of the marginal stability curve for

Consider the relation
from equation (5.18). In principle, this yields a smooth function such that σ_{R}=0 on the marginal stability curve (5.21). If we move slightly above the curve, say *M*=1/(1−*α*^{2}*Γ*)+*ϵ*, where *ϵ*<<1 is small, then we have |σ_{R}|<<1 and
5.22
again from equation (5.18). For a fixed *α*>0, because |σ_{R}| is much smaller than either Pe or *α* for *M*=1/(1 − *α*^{2}*Γ*) + *ϵ*, we have . This immediately shows that the right-hand side of equation (5.22) is positive and therefore σ_{R}>0. A similar analysis shows that when we move below the marginal stability curve, say when *M*=1/(1−*α*^{2}*Γ*)−*ϵ*, then the right-hand side of equation (5.22) is necessarily negative and consequently σ_{R}<0. Since for fixed , σ_{R} in equation (5.22) is a smooth function of *M* and *α* and σ_{R}=0 on the marginal stability curve (5.21), we deduce that σ_{R} has fixed sign on either side of the marginal stability curve and
5.23

The region above the curve corresponds to the region of instability, whereas we have the region of stability below the curve. We can use this analysis near the point (*α*,*M*)=(0,1). It follows straightaway from equation (5.23) that constitutional supercooling (*M*>1) leads to instability. That is, σ_{R}>0 for *M*>1 and σ_{R}<0 for *M*<1. Further, one can use equation (5.21) to study the effect of surface energy on the stability of the planar interface through the parameter *Γ*. We note that from equation (5.21), *M* diverges as *α*^{2}→1/*Γ*^{−}. As *Γ* increases, the marginal stability curves translate vertically upwards in the (*α*,*M*)-plane and the region of stability beneath the marginal stability curve grows in size, i.e. the planar interface is stable for a larger range of parameter values (*α*,*M*).

The parameter has a stabilizing effect, as can be seen from equation (5.20). One can rewrite equation (5.20) as
and since and for fixed *α*≠0, one can see that translates the marginal stability curves vertically upwards compared with the curve (5.21) with . Therefore, the range of stability grows in size as is increased.

The marginal stability curves for bentonite with are computed numerically and shown in figure 6. The parameter values shown there lead to the dimensionless groups Pe=0.5, *Γ*=1×10^{−5} and , with ranging from 0 to 0.5. Instability is predicted to occur when *M*>1, in agreement with experimental data (figure 7). Interestingly, as increases above about 0.45, a second short-wavelength mode develops. For this mode, which could be analysed in terms of cubic normal forms (Wettlaufer 1992), σ_{I}≠0 and hence the perturbations travel along the interface. This result could potentially explain the intriguing *spiral* and *polygonal* ice lenses observed in freezing bentonite (Taber 1929; Peppin *et al.* 2007). Indeed, the combination of the onset of unstable modes translating parallel to the interface with the local space–time variation in the compression of the colloidal suspension could lead to such patterns. Experimental measurements of for clays and colloidal suspensions may enable quantitative comparisons to be made of critical freezing velocities for the onset of such travelling waves.

### (b) Evolution equation for the interface

The structure of the linear stability diagrams suggests that the nonlinear stability of a long-wavelength mode should be assessed, and hence we examine equations (5.7)–(5.9) near the instability threshold of *M*=1 and *α*=0 by setting
5.24
where *ϵ* is a small positive number and *s* is a sign parameter, i.e. *s*=−1 (subcritical) and *s*=+1 (supercritical). We rescale the governing equations (5.7)–(5.9) and (5.12) by introducing the variables (Sivashinsky 1983; Young *et al.* 1987)
5.25
where the dimensionless variables are defined in equation (5.1). In terms of the rescaled variables (5.25), the governing system for *ϕ* becomes
5.26

We seek solutions as a power series in *ϵ* as
5.27
where *ϕ*_{0}=−1/*r*Pe e^{−rPe Z} and *H*=0 is the planar interface solution. Inserting equation (5.27) into equation (5.26) and solving the system up to O(*ϵ*^{3}), we obtain the following explicit expressions for *ϕ*_{1} and *ϕ*_{2}
5.28
and the following evolution equation for the interface shape at order *ϵ*:
5.29
where *ω*=*r*^{2}Pe^{2} *T*_{fϕϕ}/*T*_{fϕ} is a measure of the nonlinearity of the freezing temperature curve. The manipulations are similar to those in Young *et al.* (1987) and we omit the details here for brevity.

Equation (5.29) has the operator structure of the long-wavelength equation derived by Young *et al*. (1987) but with additional parameters related to the physics of colloidal suspensions. In particular the finite consolidation layer leads to a Peclet number that does not appear in their semi-infinite system, and the parameter *ω* is a consequence of the highly nonlinear liquidus curve.

We can study the linear stability of the planar interface by taking the ansatz , where *a* is real, σ=σ_{R}+iσ_{I} is the growth rate, *δ*≪1 is a positive constant and we neglect nonlinear terms of order *δ*^{2}. We substitute this ansatz into equation (5.29) and obtain
5.30
Equating imaginary terms gives
and equating real terms leads to the following explicit expression for the growth rate σ_{R} in terms of the sign parameter *s*, Pe, *r* and *Γ*:
5.31
For *s*=−1, σ_{R}<0 for all such travelling-wave disturbances and for *s*=+1, σ_{R}>0 for
5.32
i.e. the planar solution is linearly unstable with respect to perturbations with wavenumbers in the range (5.32) in the supercritical case. We note that this is consistent with the linear stability analysis where the marginal stability curve (5.21) has an asymptote at *a*^{2}=1/*Γ* and the planar interface is stable for all wavenumbers *a*^{2}>1/*Γ*.

Wollkind & Segel (1970) and Caroli *et al*. (1982) have studied the nonlinear stability of a similar set of equations in the context of alloy solidification and found that as the segregation coefficient *k*_{ s} approaches zero, the instabilities are subcritical, clearly owing to the enhancement of constitutional supercooling, which in our system is driven by the efficiency of premelting-induced rejection of particles. Abrupt transitions to instability were observed in experiments on solidification of bentonite (Peppin *et al.* 2008), indicating that the instability may indeed be subcritical. Owing to experimental uncertainties, however, it was not possible to determine if *M* was less than 1 at the transition. Future more careful experiments are needed to determine the sign of *s* at the onset of instability.

## 6. Discussion

We have studied the morphological instability of an ice–colloid interface, with vanishing segregation coefficient, including the effects of interfacial curvature, kinetic undercooling and kinetic anisotropy. Essential aspects of our work are associated with the nonlinear freezing temperature curve. We explore the regime *T*>*T*_{B} where *T*_{B} is the critical breakthrough temperature away from the shrinkage limit *ϕ*→*ϕ*_{p} in equation (2.3) so that the diffusivity *D* can be treated as a constant. The linear stability analysis of the unidirectional solidification problem reveals the growth rate of travelling-wave instabilities in terms of the Peclet number Pe, the morphological number *M*, the surface energy parameter *Γ*, the kinetic undercooling parameter and the kinetic anisotropy parameter . This allows us to explicitly demarcate the linearly stable and unstable regions in the (*α*,*M*)-plane, where *α* is the wavenumber of the travelling-wave disturbance. The weakly nonlinear analysis leads to an evolution equation for the interface shape in terms of Pe, *r*, *Γ*, and *ω*, where *ω* is a measure of the nonlinearity of the freezing temperature curve, or colloidal liquidus.

It is useful to summarize some implications of these results for the tendency of segregated ice to form during freezing of colloidal suspensions and soils (‘frost susceptibility’). In soils, the formation of segregated ice is often undesirable owing to deleterious effects on roads and buildings (frost heave and/or thaw consolidation). On the other hand, segregated ice, or other solid, is desirable in industrial colloidal suspensions, as it enables the freeze fabrication of novel microaligned porous materials. Thus, it is necessary to establish conditions under which ice segregation will occur. A main conclusion is that the onset of constitutional supercooling is an important factor determining the tendency of segregated ice to appear. Other crucial factors, such as pre-melting, have been elucidated before. Here, we summarize what we believe to be four essential criteria for the frost susceptibility of colloidal suspensions and hence of many natural soils.

(i) The ice must pre-melt against the particle surface (Dash

*et al.*2006). This allows for the segregation coefficient*k*_{ s}to be zero at slow freezing velocities, hence enabling the ice to push particles rather than engulf them.(ii) The confining pressure on the particle matrix must be such that

*ϕ*_{0}<*ϕ*_{B}. If*ϕ*_{0}≥*ϕ*_{B}, then no ice will form in the system until the temperature drops below*T*_{B}, at which point the ice will simply invade the pore space and the morphological instability considered here will not occur.(iii) The freezing velocity must be less than the critical engulfment velocity

*V*_{c}. Otherwise, the particles are trapped by the ice owing to rapid freezing (Uhlmann*et al.*1964; Rempel & Worster 1999; Wettlaufer & Worster 2006).(iv) The freezing conditions must be such that constitutional supercooling occurs ahead of the ice interface.

Small temperature gradients and high colloid-content soils (leading to large gradients in particle fraction ahead of the interface) both enhance the tendency for constitutional supercooling to occur. If all conditions are met except for the supercooling criterion (iv), one may expect to observe a stable layer of ice growing at the cooled boundary of the soil, and segregated ice will not invade the soil interior. Condition (iv) should be treated as a bound on freezing conditions, in that it is possible that the ice interface may become unstable even before constitutional supercooling has occurred (subcritical instability). In contrast to this, the volumetric expansion caused by the density difference between ice and suspension drives a flow towards the positive *z* direction, expanding the colloid boundary layer and diminishing the tendency for constitutional supercooling to occur. The effects of latent heat and differing thermal conductivities between ice and particles, while not important in the system considered here, may modify the instability condition in other systems, and should be accounted for to ensure system-specific quantitative predictions. Many future questions present themselves, ranging from (i) solely mathematical questions regarding the quantitative and qualitative nonlinear evolution of such systems, with their possibly system-specific compressibility of the colloidal matrix, (ii) the complete study of the morphological instability near the shrinkage limit *ϕ*→*ϕ*_{p} in equation (2.3) where the diffusivity changes dramatically, (iii) the extension of the model to other particularly environmentally relevant clays and silts, and (iv) engineering patterned composite materials be they lamellar, polygonal or otherwise orthotropic.

## Acknowledgements

This research was supported by the King Abdullah University of Science and Technology (KAUST), Award no. KUK-C1-013-04, by the US National Science Foundation (0PP0440841) and by the Department of Energy (DE-FG02-05ER15741). J.S.W. is grateful for support from the Wenner-Gren Foundation, the Royal Institute of Technology and NORDITA, all in Stockholm.

## Footnotes

- Received July 23, 2009.
- Accepted September 3, 2009.

- © 2009 The Royal Society