Introduction
Sea-ice growth from sea water is a unidirectional freezing process. Under quiescent conditions freezing is normally directed downwards, with a sea-ice layer growing in thickness. The ice–water interface during this freezing process is known to be cellular, with concentrated sea water sandwiched between elongated vertical plates. This structure involves the entrapment of concentrated brine between the plates and makes sea-ice saline. First notes of this microstructure date two centuries back, describing the ice as spongy and porous (Scoresby, Reference Scoresby1815), with a vertically striated structure (Walker, Reference Walker1859) and lamellar (Ruedorff, Reference Ruedorff1861). It was also noticed that the vertical orientation of the ice plates or lamellae, in the direction of freezing, corresponds to the horizontal orientation of the c-axis, normal to the ice plates (Hamberg, Reference Hamberg1895; Drygalski, Reference Drygalski1897).
An early documentation of this structure, a drawing after a tinfoil replica from the bottom of sea ice, is reproduced in Figure 1. The figure shows two basic microstructure length scales emerging from horizontal cross-sections – the plate spacing or brine layer spacing and the grain size. Within each grain, a few centimetres in dimension, one observes elongated sub-grains that are parallel within each grain and spaced by less than a millimetre. So far the most frequently used technique to determine the plate spacing is via optical micro-graphs of thin sections (Weeks and Ackley, Reference Weeks, Ackley and Untersteiner1986; Shokr and Sinha, Reference Shokr and Sinha2015). More recently X-ray micro-tomographic imaging has been used (Maus and others, Reference Maus2015). Figure 1 also shows a slice from a 3-D image obtained by this technique.
Left: Drawing after a tinfoil replica from the bottom of sea ice (Drygalski, Reference Drygalski1897) showing the crystal structure of grains and sub-grains. Right: Slice from a 3-D X-ray micro-tomographic image near the bottom of ice grown in a tank study (S. Maus, unpublished), revealing the same features.
The first quantitative descriptions of the plate spacing, a 0, ranged from 0.2 to 0.5 mm for thin ice (Fukutomi and others, Reference Fukutomi, Saito and Kudo1952; Anderson and Weeks, Reference Anderson and Weeks1958; Weeks and Hamilton, Reference Weeks and Hamilton1962) and 0.5 to 1 mm for thick Arctic sea ice (Assur, Reference Assur1958; Schwarzacher, Reference Schwarzacher1959). The largest plate spacings in sea ice were 1.3–1.5 mm and have been documented by Cherepanov (Reference Cherepanov1964) for a 10–12 m thick ice island.
A search for the predictability of the plate spacing was once motivated by two aspects – its relevance for the ice salinity and its strength. Anderson and Weeks (Reference Anderson and Weeks1958) were the first to present a microstructure-based strength model, where the plate spacing was a constant. Assur and Weeks (Reference Assur and Weeks W1963) extended this study by introducing a relationship between plate spacing and growth velocity, to construct a model of sea-ice strength evolution based on growth conditions. However, the relationship was only based on limited field data. Since then there has been little interest in this problem, presumably because large scale sea-ice properties have been in the focus, and there is little knowledge about their dependence on small-scale sea-ice properties. With the ongoing development of more advanced sea-ice models the role of microstructure scales like the plate spacing is likely to receive more interest again (Weiss and Dansereau, Reference Weiss and Dansereau2017). Regarding the salinity, an early approach to relate it to the plate spacing, via structural entrapment of brine in layers between ice plates, was proposed by Tsurikov (Reference Tsurikov1965) and considered extremely interesting by Weeks and Lofgren (Reference Weeks and Lofgren1967). However, the ideas by Tsurikov (Reference Tsurikov1965) and Weeks and Lofgren (Reference Weeks and Lofgren1967) did not merge into a consistent theory of plate spacing and salinity. This is likely related to the large discrepancy between field and laboratory observations that followed. Lofgren and Weeks (Reference Lofgren and Weeks1969) performed an extensive set of laboratory experiments of freezing NaCl solutions with salt concentrations ranging from 1 to 100 g kg−1. The results differed in several ways from other observations. The obtained concentration dependence was non-linear and rather different from the one obtained experimentally by Rohatgi and Adams (Reference Rohatgi and Adams1967) for higher growth velocities. Most importantly for sea-ice problems, the plate spacing observations from Lofgren and Weeks (Reference Lofgren and Weeks1969) were a factor of 2–3 smaller than in the field (Weeks and Hamilton, Reference Weeks and Hamilton1962; Assur and Weeks, Reference Assur and Weeks W1963). The authors attributed the considerable complexity and deviation from field data to uncontrolled convection processes in their laboratory set-up.
Nakawo and Sinha (Reference Nakawo and Sinha1984) found a rather different inverse correlation between growth velocity and plate spacing, and were also able to show that plate spacing strongly correlates with salinity in the ice. The authors also suggested, based on limited data, a dependence between crystal orientation and plate spacing. In a later study, the results were basically confirmed (Sinha and Zhan, Reference Sinha and Zhan1996), yet the growth velocity range covered was still limited to $0.5 \ndash 2\, {\rm cm}\, {\rm d}^{-1}$
. The authors proposed that a much larger set of ice cores would be needed to resolve the problem of how plate spacing and salinity depend on growth conditions.
In the latest study on the topic the author had compiled available observations of plate spacings from other studies and proposed a model for the plate spacing (Maus, Reference Maus2007a, Reference Maus2007c). The model has two main ingredients: (i) it is based on morphological stability theory (MST) developed by Mullins and Sekerka (Reference Mullins and Sekerka1964) and (ii) contains a parameterisation of convection due to brine release. With only a few parameters adjustable within plausible ranges it is able to predict the observations over 5 orders of magnitude in the growth velocity. It was further shown by the author that a model for the plate spacing allows for prediction of the salinity of sea ice (Maus, Reference Maus and Leppäranta2008). Compared to other convection-based approaches that do not consider microstructure to predict salinity summarised by Worster and Rees Jones (Reference Worster and Rees Jones2015), such a micro-structural model depends on less adjustable parameters. Strength and salinity are likely not the only physical sea properties influenced by the plate spacing. As first discussed by Anderson and Weeks (Reference Anderson and Weeks1958), the plate spacing is the starting point of an evolving sea-ice microstructure, and a concise prediction of this point is essential to understand the role of sea ice in the environment (Weeks and Ackley, Reference Weeks, Ackley and Untersteiner1986; Weeks, Reference Weeks2010; Shokr and Sinha, Reference Shokr and Sinha2015).
In the current study I present a simplified model that appears to be more consistent with the observations. To the data compilation of plate spacings from Maus (Reference Maus2007a, Reference Maus2007c) the results of some more studies are added (Fukutomi and others, Reference Fukutomi, Saito and Kudo1952; Paige, Reference Paige1966; Jeffries and others, Reference Jeffries, Weeks, Shaw and Morris1993; Sinha and Zhan, Reference Sinha and Zhan1996). I also suggest a different presentation of the results from Lofgren and Weeks (Reference Lofgren and Weeks1969), rather using the temperature gradient to relate plate spacings to growth velocity. This removes most of the inconsistency between the dataset from Lofgren and Weeks (Reference Lofgren and Weeks1969) and all other observations. The data compilation and model results are then analysed to identify processes that lead to the variability in the plate spacing, and to derive a relationship between growth velocity and plate spacing.
Observations of plate spacing
Table 1 summarises studies where both plate spacings and growth velocities were documented, or may be deduced from observations and images. Some of these studies were already mentioned in the Introduction. Most of these studies reveal an increase of plate spacing with decreasing growth velocity. In practice growth velocity often decreases with thickness which leads to an increase of plate spacing with depth. Such a relationship has been reported in other studies (Tabata and Ono, Reference Tabata and Ono1962; Paige, Reference Paige1966; Lange, Reference Lange1988). In a few studies no such depth dependence was found (Gow and Weeks, Reference Gow and Weeks1977; Jeffries and others, Reference Jeffries, Weeks, Shaw and Morris1993). However, as growth velocities in the latter studies were not reported, this behaviour might be related to some unusual growth-thickness relationship (e.g. persistent temperature decrease after freeze-up) and/or to the general scatter of observations.
Studies of plate spacing versus growth velocity
In all studies except for those of Rohatgi and Adams (Reference Rohatgi and Adams1967) and Milosevic-Kvajic and others (Reference Milosevic-Kvajic, Kvajic and Brajovic1973), denoted with $\Uparrow$
, the ice growth was downwards, in the direction of gravity. The columns from left to right give the study reference, location, typical ice thickness, growth rate and water salinity (if reported), followed by the observed range of plate spacings and the number of observations for which growth velocities were reported or could be deduced from ice thickness observations. The bracketed values indicate the total number of observations, including those for which ice growth velocities are unavailable or, in the study from Lofgren and Weeks (Reference Lofgren and Weeks1969), those where solution salinities were rather different from sea water.
Table 1 lists the studies of the plate spacing in chronological order: Fukutomi and others (Reference Fukutomi, Saito and Kudo1952) reported plate spacings for a 1–2 cm thick initial ice crust. They also reported the thermal gradients in the ice, which here were converted to growth velocities. Weeks and Hamilton (Reference Weeks and Hamilton1962) reported five sets of plate spacing and growth velocity for natural sea ice at Barrow, Alaska and six data points for laboratory grown NaCl ice. Paige (Reference Paige1966) reports the vertical variation of plate spacings for 3 m thick land-fast ice for three stations in his Figure 15 with as much as 31 data points. However, the thickness evolution for the stations in Figure 3 of that paper only allows to estimate the growth rate for the two bottom-most observations at station 2. Rohatgi and Adams (Reference Rohatgi and Adams1967) reported plate spacings in ice growing upwards from different aqueous solutions in the laboratory. The four data points for a 0.25 N NaCl solution were picked from their Figure 16 (not using data at much higher and lower concentration than sea water). Growth rates were estimated by assuming linear temperature gradients between a −70°C cold surface and the freezing temperature. Lofgren and Weeks (Reference Lofgren and Weeks1969) measured plate spacings for ice grown at different NaCl concentrations ($1 \ndash 100\, {\rm g}\, {\rm kg}^{-1}$
NaCl) in a cylindrical tank of 14 cm diameter. These authors estimated growth velocities by fitting a 4th order polynomial to the thickness data. Here we selected the four runs (at $15 \ndash 30\, {\rm g}\, {\rm kg}^{-1}$
concentration NaCl) that most closely reflect natural growth conditions with 23 data points in total. Milosevic-Kvajic and others (Reference Milosevic-Kvajic, Kvajic and Brajovic1973) report plate spacing for ice grown upwards from 1% NaCl aqueous solutions. As these authors did not explicitly report growth velocities, we use the temperature gradient to estimate the latter. Nakawo and Sinha (Reference Nakawo and Sinha1984) performed a detailed study of plate spacing observations in 1.1 m thick sea ice in Eclipse Sound, including crystal orientation to be discussed below. The data of growth velocity and plate spacing analysed here were obtained from their Figure 10. From Gow and others (Reference Gow, Arcone and McGrew1987) two values for an outdoor tank experiment (CRRELEX-84) are estimated from Figure 6, with growth velocities from Figure 3 (see also Arcone and others, Reference Arcone, Gow and McGrew1986); one value for CRRELEX-85 was obtained from Figure 16 with velocity from Figure 13. Jeffries and others (Reference Jeffries, Weeks, Shaw and Morris1993) observed plate spacings in ice cores from different places in McMurdo Sound. They distinguish between regimes of strict congelation (most studies reported here) and interstitial congelation ice that formed between platelets (The term platelets is used for crystals that either settle from the water column or grow at the ice–water interface (e.g., Dempsey and others, Reference Dempsey2010), not to be confused with the plate spacing of lamellae, the current study's subject.). As these authors found little difference between plate spacings formed under simple congelation and those measured between platelets (their Fig. 15 and Table 4), we do not distinguish between these growth modes. All plate spacings from Figure 14, were averaged for 10–15 cm spaced levels between 1.4 and 2.25 m depths, for which ice growth velocities were obtained from Figure 3 of that paper. Nghiem and others (Reference Nghiem1997) report two plate spacing – growth rate data points from a tank experiment and one from a natural lead in their Figure 6 (showing both a 0 and V). Cole and others (Reference Cole, Dempsey, Shapiro, Petrenko, Dempsey and Rajapakse1995) show the growth history for ice from Elson Lagoon, Barrow, with nine data points of plate spacing (Table 3) and growth velocity (Fig. 3). Sinha and Zhan (Reference Sinha and Zhan1996) report the vertical distribution of plate spacings in 2 m thick sea ice from Resolute Bay, together with growth velocity estimates based on a polynomial fit to ice thickness data. They compare these with their earlier study (Nakawo and Sinha, Reference Nakawo and Sinha1984). Here we plot their laboratory measurements from Figure 5.
Crystal structure observations from the 10 m thick ice island SP-6, on which the drifting station North Pole 6 was installed, were presented by Cherepanov (Reference Cherepanov1964). The vertical extent of the parallel fibres was very large – some lengths over 2.5 m were measured, indicating slow uniform growth of this ice. The authors also concluded that a considerable part of the ice island, likely 6–8 m, must have melted from its top. The origin of this ice island is unclear, yet the strong horizontal c-axis orientation would be consistent with very thick land-fast ice that existed over many years. The plate spacing for SP-6 was reported as 1.3 mm at 0.5 m depth, 1.4 mm at 2.5 m and 1.5 mm at 8 m depth. We make the following approach to associate these horizons with certain growth rates. We add to all depths the estimated 7 m ablation to obtain the thickness 7.5, 9.5 and 15 m, when these ice levels were at the ice–water interface. For an estimate of the ice growth rate we consider modelled and observed vertical temperature distributions (Maykut and Untersteiner, Reference Maykut and Untersteiner1971; Legen'kov and others, Reference Legen'kov, Chuigui and Yeremin1974). Based on these we assume that considerable ice growth is prevalent over 8 months of the year. Next, the result of a 1-D sea-ice model from Maykut and Untersteiner (Reference Maykut and Untersteiner1971) is considered. While models have been advanced since them, these authors performed simulations with forcing from the 1950s and 1960s, and for absence of oceanic heat flux, which we assume to allow the large thickness of ice island SP-6. Simulations for zero oceanic heat flux give 0.39 m a−1 bottom accretion (and surface melting) when the 5.7 m equilibrium ice thickness is reached. For the thicker ice from SP-6 this equilibrium growth rate is scaled by the thickness ratios 5.7/7.5 and 5.7/9.5 and 5.7/15, and multiplied with 12/8 to account for growth restricted to the coldest 8 months. Such a scaling is not exact, but accounts for the dominating effect of increased thickness. One then arrives at typical growth velocities of 0.12, 0.09 and 0.06 cm d−1 for the 0.5, 2.5 and 8 m horizons. Due to the assumptions made, model uncertainties and the unknown origin of ice island SP-6, these values are the most uncertain, but probably valid within a factor of 2.
The largest a 0 to date of 5 mm has been documented for the bottom of an Antarctic Ice Shelf (Zotikov and others, Reference Zotikov, Zagorodnov and Raikovsky1980). The authors obtained two consistent estimates of the growth velocity. A value of 2 cm a−1 was based on the conjecture of growth of a distinct bottom layer since a drilling event with a flame jet, while a slightly smaller number of 1.6 cm a−1 was estimated from the temperature gradient.
In Figure 2 all observations are shown. It is seen that, while observations agree reasonably well for certain growth velocities, the data from Lofgren and Weeks (Reference Lofgren and Weeks1969) deviate from all other datasets. I believe that, most likely, in their laboratory experiments there were strong lateral heat fluxes as well as thermal convection present, that led to low growth velocities, even when the temperature difference between the ice surface and the tank water was large. However, if one uses the temperature gradient in the ice to compute a ‘virtual growth velocity’ these data match observations at similar growth velocities much better. In this study we will thus only show these virtual growth velocity data points from Lofgren and Weeks (Reference Lofgren and Weeks1969) when compared to model predictions. Note that also for the laboratory dataset from Milosevic-Kvajic and others (Reference Milosevic-Kvajic, Kvajic and Brajovic1973) growth velocities are estimated from the temperature gradient. For the latter study growth velocities were also obtained from their Figure 4, which resulted in slightly (10%) different values (not shown).
Plate spacing observations from various sources summarised in Table 1. (a) The three coloured curves are empirical fits by Assur and Weeks (Reference Assur and Weeks W1963), Lofgren and Weeks (Reference Lofgren and Weeks1969) and Nakawo and Sinha (Reference Nakawo and Sinha1984) based on the data points with the same colour. The grey curve from Maus (Reference Maus2007a) is based on a morphological instability model prediction. In addition to the blue open circles showing the fitted growth velocities from Lofgren and Weeks (Reference Lofgren and Weeks1969) these data are also presented by estimating the growth velocities based on the temperature gradient (blue dots, see text for the justification).
Several investigators fitted their observations empirically and the figure summarises these fits: Assur and Weeks (Reference Assur and Weeks W1963) assumed a relationship a 0 ~ V −1/2 arguing with findings from the solidification of metals (magenta curve). Nakawo and Sinha (Reference Nakawo and Sinha1984) also considered research from the metallurgical literature and proposed that the solution at low velocities from Bolling and Tiller (Reference Bolling and Tiller1960) would yield a 0 ~ V −1 (red curve). In the double-logarithmic plot these power laws appear as straight lines with constant slopes −1/2 and −1. Lofgren and Weeks (Reference Lofgren and Weeks1969), in an extensive laboratory study, were not able to find a constant power law exponent, and proposed some other relationship found by trial and error, yet without a physical basis. While these three empirical fits are consistent with the respective limited datasets, they disagree strongly over most of the parameter space and none of them is able to predict the plate spacings variation over the whole growth velocity regime $0.005 \ndash 500\, {\rm cm}\, {\rm d}^{-1}$
. The prediction by Maus (Reference Maus2007c) is the only solution that is based on a physically based model (grey curve). In this model regimes of solute diffusion and convection are distinguished, predicting a behaviour that is reasonably consistent over the whole range of growth velocities.
Morphological stability theory
During the 1960s the problem of cellular pattern formation during solidification received increasing interest in the field of metallurgy. Rutter and Chalmers (Reference Rutter and Chalmers1953) formulated qualitatively the theory of Constitutional Supercooling based on the following considerations: solute, rejected at the freezing interface, lowers the freezing point and a solutal boundary layer forms by diffusion. Due to much faster thermal than solutal diffusion, the thermal boundary layer has a much larger extension than the solute boundary layer and the liquid will be constitutionally supercooled near the freezing interface. A small perturbation of the freezing interface moving into this supercooled liquid will have the tendency to grow spontaneously. This is the basic mechanism of cellular instability. The wavelength of the most rapidly growing perturbations is expected to increase with the solutal boundary layer width (and thus ∝ D). As more rapid growth allows less time for lateral diffusion normal to the solidification direction, this leads to smaller cell sizes. Rutter and Chalmers (Reference Rutter and Chalmers1953) supported all their qualitative arguments by observations. Tiller and others (Reference Tiller, Jackson, Rutter and Chalmers1953) postulated an appropriate quantitative solution of the problem. They proposed that the necessary condition for cellular growth is
wherein V is the growth velocity, D the coefficient of solute diffusion and k the solute distribution coefficient at the freezing interface, defined as C s/C int, the ratio of salinity in the ice to that at the interface in the liquid. C ∞ is the liquid salinity outside the solutal boundary layer, distant from the interface, and m the variation of freezing point with salinity (∂T f/∂S b). Hence for a constant m the freezing temperature T f is the same as −m C ∞. In the steady state these properties determine the interfacial gradient mG c in the freezing temperature, that has to exceed the actual temperature gradient G b = ∂T b/∂z in the liquid for constitutional supercooling to occur.
Although a considerable theoretical step forward, Eqn (1) was not capable of predicting the cell size of instabilities, and the first attempts to do so were not fully convincing (Tiller, Reference Tiller1963). However, a decade later Mullins and Sekerka (Reference Mullins and Sekerka1964) came up with an elegant solution to this problem, known today as Morphological Stability Theory (MST). The original analysis from Mullins and Sekerka (Reference Mullins and Sekerka1964) has later been presented with slight modifications (Sekerka, Reference Sekerka1968; Langer, Reference Langer1980; Coriell and McFadden, Reference Coriell and McFadden1993). The following description refers to the derivation by Coriell and others (Reference Coriell, MacFadden and Sekerka1985) and Coriell and McFadden (Reference Coriell and McFadden1993).
Linear stability of a planar interface
Consider the z-directional solidification of a saline solution at constant velocity V in a (x-y-z) coordinate system. Assume the limit of large ratios of thermal to solutal diffusivity (κ b/D ≫ 1). In the steady state the diffusion equations to be fulfilled for the temperature T b and T i in the (liquid) brine and (solid) ice are
and
The solute concentration C b in the liquid obeys
In the solid solute diffusion is neglected and C i
is given by C int at the interface in the liquid and the interfacial solute distribution coefficient k. The other boundary conditions at the interface are
for heat and
for solute, where V is the local solidification velocity (in x–y–z direction) and n the unit vector normal to the interface. L v is the volumetric latent heat of fusion, K b and K i are the thermal conductivities of the brine and ice and m is a linearised local liquidus slope. The boundary condition for C ∞ in the liquid far away from the interface is
Equilibrium at the interface implies
with melting temperature T f of the pure liquid (water). Γ is the Gibbs–Thomson parameter
where γsl is the solid–liquid surface free energy per unit area. Γ has unit Kelvins × metres and multiplying it with the local mean curvature ${\cal K}$
of the interface gives the freezing point depression due to surface tension. For an interface given by z = h(x, y) the linearised mean curvature is given as
All transport coefficients are assumed constant for a certain concentration range C int to C ∞. They are computed at C ∞ throughout the current study. Fluid flow due to the change in density upon solidification is neglected. The assumption of a thermal steady state is justified for large ratios of thermal to solutal diffusivity as for the ice–brine system (κ i/D ≈ 1.6 × 103 and κ b/D ≈ 2 × 102 near 0°C).
In a linear stability analysis approach the temperature and concentration fields are now written as the sum of an unperturbed part only depending on z and a perturbed part exp(Σt + i(ωx x + ωy y)). This analysis leads to the following equation for the freezing interface (defined as z = 0 in the steady state) for the perturbed state:
Therein Σ describes the time-dependent behaviour of an infinitesimal perturbation and ωx and ωy its horizontal wave number vectors. As shown in several similar treatments (Mullins and Sekerka, Reference Mullins and Sekerka1964; Coriell and others, Reference Coriell, MacFadden and Sekerka1985; Coriell and McFadden, Reference Coriell and McFadden1993) this analysis leads to the following dispersion relation
The parameters therein are
a horizontal wave number
while
is an effective temperature gradient at the interface, based on the temperature gradients G b and G i in the liquid and ice, respectively. At equilibrium freezing the solute gradient G c at the interface is given by
Equation (14) is valid in thermal steady state. It solves for those wavelengths λ = 2π/ω with Σ > 0 where the interface is unstable to infinitesimal perturbations, while it is stable for wavelengths with Σ < 0. To understand the influence of different terms in Eqn (14) it is useful, with help of the z-component of Eqn (7), to express the effective temperature gradient G eff in the form
with n k = K i/K b. This emphasises the important role of latent heat L v in the problem.
To proceed to a solution of Eqn (14) one defines two other non-dimensional parameters. These are the absolute stability parameter ${\cal A}$
from Mullins and Sekerka (Reference Mullins and Sekerka1964)
that contains the surface energy as an important property, and the ratio
which characterises the degree of constitutional supercooling. As discussed by Mullins and Sekerka (Reference Mullins and Sekerka1964), absolute stability implies ${\cal A} \leq 1$
and corresponds to the state when surface energy completely dominates the solute effect of constitutional supercooling. The parameter ${\cal G}$
is the ratio of the effective temperature gradient at the freezing interface, G eff, and the freezing temperature gradient due to the solute gradient, m G c. It needs to be <1 for constitutional supercooling to occur.
Sekerka (Reference Sekerka1965) showed that the stability of a planar freezing interface depends on the parameters ${\cal A} \leq 1$
and ${\cal G}$
, as well as the solute distribution coefficient k. For an explicit validation of this dependence he derived the following characteristic equation for the solution of the stability problem
where the variable r is defined as r 4 = 1 + (2 Dω/V)2 and determined by the one real root greater than unity of
while the wavelength (making use of λ = 2πω) at the onset of instability is given as
The condition for the onset of instability may then be written (Sekerka, Reference Sekerka1965) as
It expresses the modification of constitutional supercooling ${\cal G}$
for the case of interface instability due to surface tension. In the absence of surface tension effects, for ${\cal A}\ll 1$
and ${\cal S} \approx 1$
, the condition ${\cal G}< 1$
is sufficient to render the interface unstable. The larger the effect of surface tension, the smaller the necessary ${\cal S}$
for instability, and the higher the corresponding constitutional supercooling, given by smaller ${\cal G}$
according to Eqn (25).
Constitutional supercooling revised
A first principle result of MST is a re-interpretation of constitutional supercooling. Writing condition (25) for instability in the form
one can simplify it for two cases. If the liquid temperature gradient G b dominates over latent heat, stability condition becomes
The stability function ${\cal S}( {\cal A},\; \, k)$
is mostly in the range 0.9–1 under natural ice growth conditions, and hence the influence of surface energy on the selected wavelength is weak (Maus, Reference Maus2007a). Neglecting a slight deviation of ${\cal S}$
from unity, this equation may be compared to the constitutional supercooling condition Eqn (1) from Tiller and others (Reference Tiller, Jackson, Rutter and Chalmers1953). For K i = K b, equal thermal conductivities in ice and brine, these equations are equal. For the ice–brine system however, with K i/K b ≈ 4, the interface will be less stable, as G eff is smaller than G b, and a factor ≈5/2 less constitutional supercooling is needed for instabilities to grow.
For the second case it is assumed that 2G bK b is small compared to L vV – corresponding to zero oceanic heat flux. Then Eqn (26) simplifies to
where the growth velocity is no longer directly present except through its influence on the surface tension stability function ${\cal S}$
. Hence, the interface stability becomes independent of V when the temperature gradient in the liquid is small, i.e. the case of low oceanic heat flux. In summary, constitutional supercooling implies the following qualitative aspects:
• A positive liquid temperature gradient G b stabilises the interface which is the classical constitutional supercooling.
• Latent heat L v stabilises the interface as it has to be removed before freezing can proceed.
• The larger the conductivities K i and K b, the more rapidly the latent heat is removed, the less stable the interface.
• A faster removal of solute due to a larger solute diffusivity D stabilises the interface.
• Constitutional supercooling increases with concentration C ∞ and liquidus slope m, which both are destabilising.
• Surface tension stabilises the interface. This effect is only relevant at small wavelengths (high growth velocities) and implied in the decrease of ${\cal S}< 1$
.
Principal results – onset of instability
Before proceeding to a solution of the above equations, some principal results are illustrated (To obtain quantitative results we use the thermodynamic properties of NaCl solutions given in the appendix.). Equation (28) may be used in the limit of low growth velocities, when ${\cal S} \approx 1$
, to estimate the solution salinity at which a planar interface will become unstable. Assuming a bulk interface value k ≈ 10−3 for dilute NaCl solutions (Tiller, Reference Tiller1963), one finds that cellular instabilities should occur at solution salinities C ∞ above ≈0.0013 g kg−1.
Next, the behaviour of Eqn (14) is illustrated in Figure 3 for an ice growth velocity of 10−4 cm s−1 and with the liquid temperature gradient G b set to zero. Focus first on the solid curve obtained for a value of C ∞ = 0.0015 g kg−1, slightly above the stability limit. The figure shows three wavelengths that are characteristic for the solution. Moving from the left to larger wavelengths, the perturbation growth rate becomes positive, then reaches a maximum and becomes negative again at higher wavelength. These three wavelengths correspond to the surface-energy dominated lower bound ($\lambda _\Gamma$
), the wavelength where the growth rate is a maximum ($\lambda _{\max } = 3^{1/2}\lambda _\Gamma$
), and the longest unstable wavelength controlled by diffusion (λ D). Next, the solution is illustrated for a solute concentration very close to the onset of instability for a planar freezing interface ≈0.0013223 g kg−1. Under this condition the wavelengths $\lambda _\Gamma$
, λ max and λ D are very close to each other. At the critical concentration close to the onset of instability the Σ = 0 line is touched at a single wavelength, which is termed λ mi, the marginal stability wavelength. It has been noted by Weeks and Ackley (Reference Weeks, Ackley and Untersteiner1986) that, if such low salt concentration is sufficient to create constitutional supercooling for the onset of instabilities, it can be expected in any normal lake water. An interesting early observation relating to this topic has been made by Quincke (Reference Quincke1905a, Reference Quincke1905b) based on the analysis of ice frozen from dilute saline solutions. He observed that a remarkable ‘milkiness’ or ‘foam-like’ character of the ice was already observed at a solution concentration of ≈0.001 g kg−1 NaCl.
Characteristic wavelengths at the onset of instability. The solid curve presents the result slightly above the critical concentration of instability, the dashed curve gives it very close to this value. At the transition the three wavelengths $\lambda _\Gamma$
, λ max and λ D are equal to the marginal stability wavelength λ mi.
Approximate wavelength solutions
In general, Eqn (22) needs to be solved numerically, yet one can derive limiting cases (Coriell and others, Reference Coriell, MacFadden and Sekerka1985). As cell spacings during solidification are found to decrease with V −b, yet with an exponent b < 1 (e.g. Tiller, Reference Tiller1963), one finds ωD/V ≫ 1 at high growth velocities far from the onset of instability. This leads to (Coriell and others, Reference Coriell, MacFadden and Sekerka1985)
If one further sets ${\cal G} \ll 1$
, the case of high supercooling, one obtains by inserting ${\cal A}$
from (Eqn (20)) the relationship λ Γ ~ (V/D)−1/2, which is an often assumed power law (e.g. Assur and Weeks, Reference Assur and Weeks W1963). The corresponding most rapidly growing wavelength is λ max ≈ 31/2λ Γ (Coriell and others, Reference Coriell, MacFadden and Sekerka1985).
An approximation for the onset of instability is obtained by assuming ${\cal A}^{1/6} \ll 1$
. Equation (22) may then be written as
while the critical wavelength is now given by
Hence, at the onset of instability the wavelength is expected to scale as (V/D) −2/3.
Modifications of MST for sea ice
The classical MST theory as outlined so far describes a planar freezing interface. Based on a steady state, it cannot describe how and under what conditions the interface evolves into a system of plates or cells with deep grooves between them. However, one can make the assumption that, after changing from a planar to cellular state, the freezing interface continues to operate at the onset of instability, as proposed by Maus (Reference Maus2007a). The new cellular interface then needs to be described in terms of effective thermal conductivities and temperature gradients at the interface. For lower growth velocities and growth in the direction of gravity (most natural sea ice) also a modification of the purely diffusive equations due to solutal convection is needed.
Marginal stability – maximum k
When the interface becomes unstable the interface solute distribution coefficient k is no longer a fixed parameter. To first order, neglecting the ice–brine density difference, k now reflects the brine volume fraction or porosity at the interface. As the interface concentration and solute gradient decrease with k one wants to know the maximum k at which the interface can stay marginally unstable. Focusing on small ${\cal A}$
far from absolute stability the solution can be approximated by Eqn (30). We further define two new parameters
and
to rewrite Eqn (30) in the form
which is a fifth-order polynomial to be solved for k. The basic result of this solution can be illustrated in the limit of ${\cal G} = 1$
. Then Eqn (33) becomes $k = ( 1-{\cal G}_0) ^{-1}$
. Assuming now zero temperature gradient in the liquid (no oceanic heat flux) one obtains
In this approximation k only depends on material properties and the solute concentration C ∞. Inserting the properties of an 35 g kg−1 NaCl solution at its freezing point one obtains k ≈ 0.9693. With ΔC = C int − C ∞ one has
for the interface salinity increase over the far field value. ΔC is to first-order independent of growth velocity or salinity (Note that, as L v, D, K i and K b and m are all temperature and salinity dependent, ΔC will also be so.). For C ∞ = 35 g kg−1 one obtains ΔC = 1.1 g kg−1, and a constitutional supercooling of mΔC = −0.069 K.
Near-interface temperature gradient
Morphological stability depends essentially on temperature and solute gradients at the interface (Fig. 4). To illustrate how these are changed for a cellular interface, consider a planar interface for the case G b = 0. The ice growth velocity is then given as V = G iK i/L v. When the growth is dendritic or cellular, this relationship is changed. Now heat is not only conducted away from the interface through the solid ice but also through the liquid brine. This changes the relationship between the growth velocity V and the bulk temperature gradient G i at the dendritic interface. G i in the tips is smaller to the degree to which the heat flow through the liquid enhances the growth V of ice dendrites. For a given growth velocity the temperature gradient G i at the interface may thus be less than for a planar ice interface. To account for this K i at the interface may be replaced by an effective thermal conductivity K eff that obeys
To keep the problem solvable one needs to parameterise K eff in terms of k. Assume a geometry where the edges of knife-like plates are rounded and have a certain tip radius, and associate k with the average solid ice volume fraction ϕ tp of these tips. The salinity in this thin interface layer can then be expressed as
where q = ρ b/ρ i is the ratio of densities of brine and pure ice. Combining this equation with C i = k C int, the solid fraction ϕ tp in the cell tips becomes
Let ϕ tp be the average ϕ s integrated from the interface (ϕ s = 0) to the value ϕ s = ϕ rt at the root of cell tips, one radius away from the interface. Then the latter is given as
Apparently this scaling only makes sense while ϕ tp < 2/π, as at higher solid fractions the assumed geometry is not valid and would imply ϕ rt > 1. However, in this case the present concept can still be applied by setting ϕ rt = 1. The geometric details will not play a large role in this situation, as the effect of heat conduction through the liquid is small anyway.
Sketch of the freezing interface with plate spacing a 0 and the solutal boundary layer ahead. On the right hand the temperature gradients are indicated, with the weaker gradient in the ice near the interface. Within the solutal boundary layer the concentration is changing from C int to C ∞. The corresponding freezing temperature is lower than the actual temperature which implies that this layer is constitutionally supercooled (CS). Solid fraction ϕ tp and C int within the tip regime relate to k that sets this supercooling. The vertical level one radius up, termed root of the tip, has solid fraction ϕ rt. When the boundary layer D/V reaches its critical thickness for the onset of convection, the corresponding k nu and ϕ nu are defined slightly longer up in the ice. The boundary layer than receives solute also from a regime between the cells. Note that the solid fraction profiles above ϕ rt are only tentative.
The effective thermal conductivity K eff for the tip regime can now be written as
where w K is the fraction of the heat flow through the liquid that contributes to the growth of dendrites. In this study we consider three simple bounds for w K. At the lower bound w K = 0 all heat flowing through the intracellular liquid is drawn from the infinite water reservoir. Then no liquid conduction effect on G i is present, K eff = K i, and the temperature gradient is a maximum. For w K = 1 all heat conducted contributes to enhanced dendritic growth, implying a minimum temperature gradient. A third solution is constructed from simple geometric arguments based on circular cell tips. The heat drawn vertically upwards through the brine (at ϕ rt) must come from a volume between two cell tips, bounded laterally by the tip surfaces, and vertically by the ϕ = 0 interface. Assuming that the heat is drawn equally through these boundaries one has
which will be taken as the standard solution to be used in our calculations.
Solutal convection
Solutal convection is incorporated into the present MST approach based on earlier study on the convective instability of an infinite fluid layer. For the freezing of sea water this has been addressed theoretically and experimentally (Foster, Reference Foster1968, Reference Foster1969a). The critical parameter to evaluate the onset of convection is the solutal Rayleigh number
based on the concentration difference ΔC = (C int − C ∞). β, ν and g are the haline contraction coefficient, kinematic viscosity and gravity acceleration. H is the layer thickness that below will be identified with D/V. When Ra s exceeds a critical value, convection sets in and the solute flux increases above its diffusive value F 0 = DΔC/H. For this enhancement we use the well-known relation
for the Nusselt number Nu. The salt flux F s then becomes
It is independent of the layer thickness, which is the classical 4/3-flux law (e.g. Turner, Reference Turner1973; Grossmann and Lohse, Reference Grossmann and Lohse2000). The problem of cooling or density increase from one boundary differs from the classical two-sided Rayleigh–Benard problem, and in that c nu is a factor of 24/3 larger for the one-sided case, as can be seen by inserting half thickness and density difference in Eqn (43). We use the range c nu = 0.15 ± 0.02 (corresponding to the range 0.052–0.068 for the two sided problem) found in experimental studies (Turner, Reference Turner1973; Katsaros and others, Reference Katsaros, Liu, Businger and Tillman1977; Selman and Tobias, Reference Selman and Tobias1978; Goldstein and others, Reference Goldstein, Chang and See1990) to account for a plausible range in c nu.
The increased solute transport at the interface implies an effective k nu based on the Nusselt number that is
This k nu may be imagined to be related to a solid fraction increase further up in the ice, above the root of the cell tips. It is then in turn expected to affect ϕ rt and thus the heat conduction and temperature gradient given by Eqn (41). To account for this effect we use Eqn (39) and compute
We then define a modified solid fraction ϕ rt′ that replaces ϕ rt in Eqn (41) as
With this formulation the effective thermal conductivity near the cell tips, K eff(k), is now a function of the solid fraction ϕ rt(k) near the tip and the solid fraction ϕ nu(k, Nu) further up in the ice. While ϕ rt in the tip regime is given by k and the minimum constitutional supercooling, ϕ nu is a higher solid fraction further up in the ice and associated with solute transport due to convection. The convective solute transport takes place in a thin horizontal boundary layer situated ahead of the interface. As the interface is cellular this layer will be deformed upwards between the cells. The downward solute transport implies lateral cell thickening and increase in the horizontally averaged solid fraction that decreases the effective thermal conductivity K eff(k) at the interface, Eqn (41). This effect is here approximated through an effective solid fraction ϕ′rt. While this approach ignores the details of lateral cell growth, it keeps the problem solvable and represents the basic feedback: the higher solid fraction leads to lower K eff, which in turn implies a higher temperature gradient G i near the interface, Eqn (37). One effect of convection is thus, via increased solid fraction, to let K eff approach K i. However, the main effect of convection is to increase the solute gradient G c and the solute transport from the ice interface into the sea water.
The physical process behind Eqn (45) is the build-up of a diffusive boundary layer, until it convects and solute is transported away in plumes, followed by a new build-up. The critical time for this build-up is related to a critical Ra s at the onset of convection. Foster (Reference Foster1968) obtained such a relationship numerically. Defining the diffusional distance as H = (D t c)1/2, he obtained an equation for the critical time t c at the onset of convection
Inserting H = (D t c) 1/2 in Eqn (43) one obtains a corresponding critical Ra sc = 593/2 ≈ 453. Results from other studies with different boundary conditions and Schmidt numbers (Sc = D/ν) cover the range 100 < Ra sc < 500 (Mahler and Schlechter, Reference Mahler and Schlechter1970; Onat and Grigull, Reference Onat and Grigull1970; Gresho and Sani, Reference Gresho and Sani1971). To be consistent with Eqn (44) we use $Ra_{sc} = c_{nu}^{-1/3}$
, which ensures that Nu exceeds 1 above Ra sc. This gives Ra sc ≈ 296 for c nu = 0.15 and 204 < Ra sc < 455 for the assumed range in c nu.
It is ΔC = (C int − C ∞) that enters Eqns (43) and (49). ΔC in turn is obtained from MST assuming marginal stability and minimum supercooling, given by Eqn (36) yet modified for surface tension effect (recall the role of ${\cal A}$
). The computation of critical time t c as well as Rayleigh number Ra thus becomes possible through this boundary condition of minimum constitutional supercooling. For the present steady state solution the critical time will rather represent the characteristic time scale for intermittent build-up and breakdown of the boundary layer (Foster, Reference Foster1968). However, it is well known that D/V is the appropriate length scale of a diffusive boundary layer in solidification problems (e.g. Trivedi and Kurz, Reference Trivedi and Kurz1994) (This may be seen by differentiating the diffusion length (2Dt)1/2 with respect to t to obtain V = (D/2t)1/2, which leads to D/V = (2Dt)1/2.). Inserting D/V in Eqn (43) one can compute, for a given critical Ra s value, the critical growth velocity V c at which convection sets in as
This gives V c ≈ 15 cm d−1 for Ra sc = 296. Based on the 4/3-law parameterisation, for growth velocities below V c the salt flux will be given by F s = DΔC/(D/Vcr) = ΔC V c and only depend on ΔC, not on the growth velocity V. Given the process of intermittence, a concise incorporation of the convection scaling into MST would be rather complex and is not attempted here. However, one can make an approach based on the approximate solution (31) for the wavelength at the onset of instability. Assuming G eff ≈ m G c for marginal stability it may be written in the form
This equation contains three length scales on the right hand side: a thermal length (ΔT/G eff), the solutal diffusion length D/V, and the capillary length Γ/ΔT. Assuming that after the onset of solutal convection the solutal length scale is given by D/V c, the main effect of convection will be a change from Eqn (31) to
with V c given by Ra sc. This implies that accounting for convection leads to a weaker dependence of λ mi on the growth velocity, and that the expected scaling law has the form λ mi ~ V −1/3.
Results
In summary, there are three modifications of MST that we propose for sea-water freezing:
(i) It is conjectured that the cellular interface is in an operating state of marginal stability, where the constitutional supercooling at the interface is at a minimum. This corresponds to maximum k obtained by solving Eqn (34) with an approximation by Eqn (35).
(ii) The temperature gradient at the interface should be modified for high solid fractions. A simple geometric model for the cell tips is assumed and the effect is implemented in dependence of k through Eqn (41).
(iii) Convection changes the solute flux and the thickness of the solutal boundary layer. This happens below a critical velocity V c ≈ 15 cm d−1. Also convection is parameterised in terms of k, leading to a modified solution for the wavelength at the onset of instability (Eqn (52)).
In the following the model predictions will be shown for pure diffusion and with convective parameterisations included, typically resembling the cases of upward growth and downward growth (that is, natural sea ice).
Diffusive solution for fixed k
In Figure 5 we compare the results of MST wavelength predictions with observations of sea-ice plate spacings. The model results are shown as a stability balloon in growth velocity – wavelength space. Only diffusion is taken into account and a fixed k and C ∞ = 35 typical for sea water are used. Inside the balloon a planar interface is unstable, outside it is stable. The results are shown for three temperature gradients in the liquid, being expressed as heat fluxes (Q = G bK b), and respectively representing the cases of very small to moderate to very large oceanic heat fluxes during the sea-ice growth season. Note that here k was fixed at 0.15, as this value implies a concentration of C ∞/k ≈ 233 g kg−1. Lower values of k would imply a concentration exceeding the eutectic concentration and are therefore not considered. The wavelengths relevant for MST are indicated – the short-wave branch on the left related to surface free energy (λ Γ), the red curve giving the wavelength of maximum growth rate (λ max), the long-wave branch on the right controlled by diffusion (λ max). The wavelength at the onset of instability, (λ mi) is given approximately in the middle of the bottom of the balloon. Here all branches bifurcate. They merge again at the top of the balloon at the absolute stability limit indicating that at velocities larger than ≈ 102 cm s−1 the interface is stable.
Stability bounds based on MST for different heat flux from the liquid/ocean, shown together with plate spacing observations from sources in Table 1 and Figure 2. The extreme branches for $\lambda _\Gamma$
and λ D are indicated, the wavelength of maximum growth rate λ max as a red dotted curve. The onset of stability (for increasing growth velocity) takes place at λ mi at the bottom of the balloons where the stability branches bifurcate.
It is seen that the heat flux can affect the stability, yet that for these settings it would be of little relevance for sea-ice growth. Only the slowest growing ice falls out of the stability balloon for high heat flux. In general, the heat flux is seen to have no effect on the wavelength selection when the interface is unstable. It just moves the onset of stability to higher growth velocities. Another message from the figure is that the observations are one to two orders of magnitude away from the short and long-wave branches. Does this mean that the plate spacing cannot be predicted?
Marginal stability – the role of k
In Figure 6 the neutral stability balloon is plotted in a different way. Now a very small oceanic heat flux is fixed, yet k is treated as the variable parameter. The larger k, the better the agreement between observations and theory becomes. The essence of this figure is that it supports the marginal stability hypothesis of k being maximised, which means that constitutional supercooling is a minimum. Indeed when using k ≈ 0.96 close to the estimate obtained with Eqn (34) the stability balloon shrinks further and the agreement between theory and observations starts becoming reasonable.
Stability bounds based on MST for different values of k, the solute distribution coefficient at the ice–water interface. The outer stability branches for $\lambda _\Gamma$
and λ D (refer to Fig. 5) move towards each other and the whole stability balloon shrinks while k increases. The stage when the balloon shrinks to a single point corresponds to some maximum k and the marginal stability wavelength λ mi.
Diffusive solution – λ mi
Figure 6 has demonstrated the importance of marginal stability in the problem, corresponding to the maximum k and minimum ΔC. We now obtain the corresponding solution for λ mi in two ways. First we solve Eqn (14) numerically for k, finding the transition to positive perturbation growth rates by iteration. Second, we use the approximation Eqn (31) with k obtained via (34). Both solutions are calculated for the discussed range in temperature gradients at the interface, by employing Eqn (41) for three values of w K. The range of these purely diffusive solutions is shown in Figure 7. It is seen that the numerical prediction (dashed curve) and the approximation (red dashed-dotted curve) do not differ much, the difference being largest at high velocities.
Marginal stability wavelength λ mi predicted by MST with plate spacing observations from sources in Table 1 and Figure 2. The shading shows the range of the purely diffusive numerical model prediction based on the interfacial temperature gradients, where the upper bound corresponds to the minimum temperature gradient (w K = 1 in Eqn (41)). The dashed curve is the proposed standard solution (temperature gradient 0.74 weaker than the maximum). The red dash-dotted curve is the approximation (31).
Note that the dashed curves correspond to the standard case based on Eqn (37), and a temperature gradient 0.74 weaker than the maximum. The solution with w K = 1 corresponds to a temperature gradient at the interface being a factor of 5 smaller, to be discussed further below. Above a growth velocity of ≈20 cm d−1 there is, given the data scatter and range in predictions, reasonable agreement between the predicted λ mi and the observed plate spacing. Most observations fall between the predicted temperature gradient bounds. It seems that the observations are closer to the theoretical upper bound corresponding to the minimum temperature gradient. It is also evident that at lower growth velocities the diffusion theory predicts wavelengths that are far too large when compared to observed plate spacings.
Solution with convection
To account for convection we solve Eqn (34) for k and ΔC, evaluate the critical growth velocity V c for the onset of convection through Eqn (50), and compute the salt flux via Eqn (45). The critical wavelength λ mi is then obtained from the approximate Eqn (52). The plate spacing predicted in this way is shown in Fig. 8. The shading now indicates the bounds from the minimum and maximum temperature gradients corresponding to 0 < w k < 1 and from the assumed uncertainty in the parameterisation of the Nusselt number 0.13 < c nu < 0.17.
Marginal stability wavelength λ mi predicted by MST with the parameterisation of solutal convection in the approximate Eqn (31). The shading shows the range of solutions based on bounds for convection and the interfacial temperature gradient: The upper bound represents lower bounds for the temperature gradient (corresponding to w K = 1 in Eqn (41) and convection parameterisation c nu = 0.13, while the lower bound of the shading represents upper bounds of temperature gradient and convection strength (w K = 0 and c nu = 0.17). The dashed curve is for the standard settings (w K from Eqn (42) and c nu = 0.15). All data points as in Fig. (7).
Note that Fig. 8 shows, for clarity, only the results based on the approximate Eqn (52), not the full numerical solution shown in Fig. 7. However, for growth velocities below 10 cm d−1 these only differ by a few percent (compare the dashed black and red dash-dotted curve in Fig. 7), and the approximation is thus considered as sufficiently accurate for the regime of natural sea ice and marine ice growth.
The standard case, corresponding to our most likely temperature gradient reduction as well as c nu = 0.15, is shown as a dashed curve. At high growth velocity this dashed curve corresponds to the diffusive approximation shown as a red dashed curve in Figure 7. As mentioned, in this regime the difference between the approximation and the full numerical solution is largest. The observations of plate spacings fall slightly above the approximated upper limits of λ mi (given by the shading).
Figure 8 demonstrates, most importantly, that the MST predictions based on convection also agree with plate spacing observations at lower growth velocities. They match the observations down to the ice shelf data with a growth velocity near 2 cm a−1. Note that the different temperature gradient scenarios, given by the parameterisation of K eff near the interface (Eqn (41)) imply different k and ΔC, and thus different critical velocities for the onset of convection. While the value is V c ≈ 14.9 cm d−1 for the standard settings, the maximum and minimum temperature gradients imply 16.1 and 10.4 cm d−1. These transitions are reflected by a change in the slope in the log–log plot that corresponds to a shift from λ mi ~ V −2/3 to ~ V −1/3. This slope change behaviour is crucial for the agreement of predictions and observations.
The shaded bounds of the predictions in Figure 8 show an interesting behaviour, in that the range of predictions first increases once the critical velocity for convection onset is reached, and then decreases reaching a smaller range at low growth velocities. The explanation is that, once convection sets in, the predicted range is roughly created by the sum of uncertainties in the temperature gradient and the convective parameterisation. With increasing convection the solute transport at the interface increases, implying a lower k nu and an increase of the solid fraction ϕ nu at the root of cell tips (Eqn (47)). With larger near-tip solid fraction also the temperature gradient has to increase (for the same V) and all solutions approach the maximum temperature gradient condition w K = 0. The uncertainty that remains at low growth velocities is the uncertainty in the convection parameterisation, which has a weaker effect than the temperature gradient at the interface.
Discussion
The main purpose of the current study is the theoretical prediction of the plate spacing over the whole range of growth rates, based on MST from Mullins and Sekerka (Reference Mullins and Sekerka1964). MST is suitable for predicting the onset of instabilities on a planar freezing interface and its breakdown into a cellular interface. For a given freezing rate MST predicts that instabilities will fall between a short wave branch governed by surface energy and a long wave branch governed by diffusion, and that the dependence on growth velocity is λ D ~ V −1 for the latter branch and $\lambda _\Gamma \sim V^{-1/2}$
for the former. Wettlaufer (Reference Wettlaufer1992) presented theoretical calculations based on MST and showed that observed sea-ice plate spacings indeed fall between these bounds. However, due to the large range spanned by the bounds (see Fig. 5) he concluded that plate spacings are not predictable by linear theory and that the problem likely requires a more complex non-linear treatment.
The current study follows and extends the proposal by Maus (Reference Maus2007b). When the onset and margin of instability is considered, MST theory implies the dependence λ mi ~ V −2/3 (Coriell and others, Reference Coriell, MacFadden and Sekerka1985). The physical background for this scaling, as well as other length scales during freezing and directional solidification, has been summarised by Trivedi and Kurz (Reference Trivedi and Kurz1994). At the onset of instability, there are three length scales of equal importance for predicting the critical wavelength. Combining the notation from Trivedi and Kurz (Reference Trivedi and Kurz1994) with ours these are a thermal length, based on temperature gradient and supercooling near the cell tips (l T ~ ΔT/G eff), the solute diffusion length (l D ~ D/V) and the capillary length ($l_\Gamma \sim \Gamma /\Delta T$
). The critical wavelength at marginal stability scales as cubic root of the product of these three length scales which means $\lambda _{\rm mi} \sim ( l_{\rm T} l_{\rm D} l_\Gamma ) ^{1/3}$
. As ΔT cancels in the product $l_\Gamma l_{\rm T}$
, this leads to the proportionality λ mi ~ (ΓD/G eff V) 1/3. When the ice growth velocity mainly depends on the temperature gradient in the ice, G eff ~ V, and as Γ is velocity independent, this implies λ mi ~ V −2/3. Indeed, this is the dependence found in Figure 7.
When applying this marginal stability principle, one obtains a wavelength or plate spacing that follows from a state of minimum constitutional supercooling. This is achieved in MST by treating the interfacial solute distribution coefficient k as a free parameter of a cellular interface, not as a fixed value. For the standard simulations for sea water this leads to
The agreement of MST with observations of the plate spacing increases considerably when the marginal stability criterion is applied, as shown in Figure 7. However, there remains a mismatch between model and observations that becomes larger with decreasing growth velocity. Hence, another fundamental change in the theory is needed – the account of convection.
The influence of convection was formulated with a Rayleigh number-based framework assuming that convection is controlled by the interface salinity excess related to constitutional supercooling. Basically the approach formulates a coupling of the marginally stable k and ΔC to the 4/3-flux law for convection, Eqn (45). This allows to compute the critical growth velocity at the onset of convection (V c ≈ 15 cm d−1) and the enhancement of solute fluxes below this velocity. The fundamental change to the diffusive solution is now that, below V c for onset of convection, the solutal diffusion length, or boundary layer thickness, becomes constant (l D ≈ D/V c). From this it emerges a relationship where the critical wavelength is λ mi ~ V −1/3.
Validations of MST parameters
The present approach predicts the marginal stability wavelength by prescribing how an interface with an effective k interacts with solutal and thermal boundary layers and their respective salt and heat transport. The present predictions are consistent with the compiled dataset of plate spacings. While there is considerable scatter in the observations of a 0, in particular near a sea-ice growth velocity of 1 cm d−1, the range of plausible model parameterisations leads to a comparable range in predictions a 0. To validate and constrain the present model, it will now be discussed where the approximations made in the model are confirmed by other than plate spacing data (of, e.g. the temperature gradient), and where there is potential and need for observations. Also the observations of a 0 are discussed in more detail to answer the question if their variability in plate spacing can be related to physically justified changes in model parameterisations. The most relevant properties and processes in the theory are the solid fraction and ΔC at the interface, the near interface temperature gradient in the ice, and the parameterisation of solutal convection.
Interfacial k and ΔC
Direct observations of the concentration field near the freezing interface are difficult to obtain. The only study known to the author is by Terwilliger and Dizio (Reference Terwilliger and Dizio1970), who used micro-conductance probes to monitor the solutal boundary layer during upward freezing of NaCl solutions. From the observed concentration profiles the authors obtained k at several solution salinities from 0.29 to 5.8 g kg−1, and over the freezing velocity range 4 to 40 × 10−4 cm s−1. Solving the same equations for diffusion as in MST, the authors found k to depend on C ∞, while ΔC across the boundary layer was largely independent of the latter, in agreement with Eqns (35) and (36). The slightly lower k (and larger ΔC ≈ 2) compared to MST predictions can be explained by high liquid temperature gradients in the experiments (see Fig. 9.1 in Maus, Reference Maus2007a). Another set of observations suitable to validate the present approach has been published for the upward freezing of aqueous KCl solutions by Kirgintsev and Shavinskii (Reference Kirgintsev and Shavinskii1969). These authors performed experiments for a larger concentration range at a growth velocity of V = 4.7 × 10−3 cm s−1. k was evaluated based on sectioning completely frozen samples. Results from this indirect approach to evaluate k are more scattered, yet also agree fairly well with the present predictions (see Fig. 9.2 in Maus, Reference Maus2007a).
Looking at the concentration dependence of the plate spacing one may insert the approximation (35) into Eqn (31) to obtain
With K i + K b/(L vD) ≈ 15K −1 and mC ∞ being approximately the freezing point, one finds only a few percent change above water salinities of 10 g kg−1, yet some influence at dilute concentrations. This behaviour is consistent with the limited observations available (Lofgren and Weeks, Reference Lofgren and Weeks1969; Maus, Reference Maus2007c).
Solutal convection
Consider now the conditions for the onset of convection as derived from the approach above. Once the critical growth velocity of V c ≈ 15 cm d−1 is reached (Eqn (50)), the diffusive boundary layer starts to convect and emit plumes. Its thickness at this critical condition is approximately given as H c = D/V c ≈ 0.4 mm. The time to reach this critical thickness according to Eqn (49) is t c ≈ 4.5 min and may be associated with the frequency of intermittent convection. In addition to this time scale Foster (Reference Foster1968) also proposed a characteristic wavelength L c of the convection as
For our standard case the value obtained for L c is 2.7 mm. Are there observations available to support these numbers? Foster (Reference Foster1969a) performed sea-water freezing experiments and observed the wavelength of convection for different growth conditions by means of optical Schlieren imaging. He observed a convection cell spacing of $2\hbox {--}3$
mm, consistent with the theoretical prediction. Critical times were not well defined in his study, due to the difficulty to know exactly the onset of freezing in the presence of supercooling. In a more recent study by Middleton and others (Reference Middleton, Thomas, Wit and Tison2016) also Schlieren optical photography was used to observe convection patterns in a thin growth cell. In addition to streamers emerging from the ice's interior, the authors report on small convection cells emerging from the ice–water interface, having a wavelength of the order of 5 mm. The latter value is larger than predicted, yet the cell was very thin (3 mm) which may affect the convection behaviour. The observations in both studies are consistent with the theory proposed here. In both studies the convections cells remained at rather fixed positions.
No observations of the thin boundary layer H c mm at convection onset have been reported so far. A closer look at the supplementary material (Video 2) from Middleton and others (Reference Middleton, Thomas, Wit and Tison2016) indicates a bright thin interfacial layer that appears to be of the order of 0.5 mm. However, observations of this boundary layer remain a challenge. While Middleton and others (Reference Middleton, Thomas, Wit and Tison2016) reported that the layer associated with the influence of the plumes remained 1–2 cm in thickness during all experiments, Foster (Reference Foster1969a) traced the plumes vertically through the 25 cm deep tank.
With regards to the critical Rayleigh number for onset of convection a range of 100 < Ra sc < 500 has been proposed as most reliable based on experimental studies (Mahler and Schlechter, Reference Mahler and Schlechter1970; Onat and Grigull, Reference Onat and Grigull1970; Gresho and Sani, Reference Gresho and Sani1971). Changing Ra sc in this range will affect the results for L c, t c and H c by some 20–30%. The few observations presently available thus are insufficient to find an optimum parameterisation of convection.
With regards to the critical Ra c several authors have studies the coupled problem of the onset of convection and morphological instabilities at a planar interface. It was found that in a steady state approximation (as used in MST), convective motion sets in at rather low Rayleigh numbers when based on D/V (Hurle and others, Reference Hurle, Jakeman and Wheeler1982; Caroli and others, Reference Caroli, Caroli, Misbah and Roulet1985). For a solute distribution coefficient of k close to unity, convective instabilities will start at Ra sc ≈ 12. At first glance this seems to be a large discrepancy to values mentioned above and applied in this study. A plausible explanation has been given by Foster (Reference Foster T1969b), who observed that critical times theoretically predicted by the diffusive penetration depth tend to be too short by a factor of 4. This implies through Eqns (49) and (43) an eight times smaller critical Ra sc. As an explanation he proposed that the time disturbances need to develop is long compared to the time for growth of the boundary layer. Hence, while the first convective motion may be initiated at Ra sc ≈ 12, manifest boundary layer convection requires larger Ra, and the scaling laws for Nu noted above will remain largely unchanged. Note that these considerations have been taken by Foster (Reference Foster1968) into account when formulating the semi-empirical scaling Eqns (55) and (49).
We can finally note that the critical V c for Ra sc ≈ 12 would be ≈43 cm d−1. The only downward freezing experiments, with potential convection, at such high growth rates are those from Lofgren and Weeks (Reference Lofgren and Weeks1969). A look at Figure 8 indicates that existing plate spacing observations are insufficient to evaluate the location of the slope change. However the overall effect of an earlier (in terms of higher V) convective motion on the plate spacing is probably small, as it is the Nusselt number parameterisation that is more relevant. Also, sea ice will very seldom grow with that high growth velocities.
Near-interface temperature gradient
We have introduced K eff (Eqn (41)) as an enhanced modified thermal conductivity to account for the fact that dendrite or cell growth at the interface can be enhanced by heat flow through the liquid between plates. For the same growth velocity this implies a smaller temperature gradient at the interface. The formulation gives K eff > K i and notably differs from the application to estimate an effective bulk thermal conductivity, based on brine and ice phases. The latter is always smaller than the ice conductivity K i. In the present approach the growth velocity is prescribed and the focus is on the modification of the temperature gradient near the interface.
The three cases considered are (i) no growth enhancement at all, (ii) a plausible geometrical model (Eqn (37)) and (iii) the minimum temperature gradient where all heat drawn through the liquid brine advances the growth of ice cells. For the present standard model setting with C ∞ = 35 and k ≈ 0.97 the corresponding values of K i/K eff are 1, 0.74 and 0.20. Hence, while a plausible reduction of the interfacial temperature gradient should be 0.74, a maximum reduction of 0.2 is possible.
Detailed observations of the temperature gradient near the sea–ice interface have not been presented so far. However, when looking at temperature profiles documented for young ice growth, the near interface temperature gradients are generally reduced (Cox and Weeks, Reference Cox and Weeks1975; Niedrauer and Martin, Reference Niedrauer and Martin1979; Wettlaufer and others, Reference Wettlaufer, Worster and Huppert1997). With respect to the bulk ice temperature gradient this reduction is typically in the range 0.6–to 0.8. Hence, the proposed reduction in the near-interface temperature gradient is consistent with observations.
When convective solute transport comes into play, the temperature gradient bounds are changed. In the present formulation we use two solid fraction estimates. ϕ rt is controlled by diffusion in a tiny layer representing the tip geometry, while ϕ nu relates to enhanced solute flux by convection and the lateral freezing of cells. Increased solute transport implies a lower k nu and an increase of ϕ nu at the root of cell tips (Eqn (47)). This in turn increases the near-interface temperature gradient (for the same V). Finally, at low enough V and large ϕ nu, the maximum temperature gradient condition is approached. In this regime the predicted range for the plate spacings is related to the uncertainty in convective solute transport.
How realistic is the minimum temperature gradient approximation? We consider two settings where it might be rather relevant. The first relates to ice growth in finite sample containers. Then the heat flow drawn through the brine fraction is taken from a finite volume that can become supercooled. This supercooling applies to the liquid as a whole, it is not constitutional. It will enhance the growth of ice cells and dendrites into the lab container, yet not affect the wavelength selection in the solutal boundary layer. A look at the diffusive growth regime in Figure 8a shows that most observations above 20 cm d−1 are close to the upper limit in the plate spacing, the minimum temperature gradient limit. All these data are obtained in laboratory studies, which are most prone to this effect. The second proposed mechanism relates to crystal orientation. In the standard case heat conduction aligned with brine layers and ice columns heat conduction takes place in parallel. The ice is not advancing due to heat transport in the brine. However, this heat transport pattern will be changed when the dendrites are inclined. Now the heat transport near the cells tips approaches a serial conduction model which, for high liquid fractions gives similar increase in K eff/K i as for our minimum temperature gradient bound.
For an exact solution of the heat and solute flow near the cell tips one would need to know the tip geometry, and employ relationships between the local supercooling and the tip radius, while the present approach is using global values for interface. Observations indicate a complex 3-D pattern for fresh ice dendrite tips (Furukawa and Shimada, Reference Furukawa and Shimada1993), and similar observations would be needed for sea ice. While models to predict plate spacings based on tip radius geometry have been suggested for other systems, these seem not to perform too well for sea ice (Maus, Reference Maus2007c).
Processes affecting the plate spacing
To evaluate the processes versus observations we now focus on the growth velocity regime around 0.5–2 cm d−1, for which the data and model predictions are shown in Figure 9. This appears to be the regime with largest variability in observed late spacings. Due to the discussion in the last paragraph we focus on two candidate processes to affect the plate spacing – the role of crystal orientation and choice of model for the solute transport and Nusselt number. These are discussed based on Figure 9, where observed plate spacings are shown along with model predictions. The light grey shading highlights the range due to the maximum and minimum interface temperature gradient, while the darker shading the effect of convection parameterisation.
Model predictions from MST with solutal convection implemented. Data points for the plate spacing are shown with the same symbols as in Figure 8. The light shading shows the range of the predictions related to bounds for interfacial temperature gradient, while the dark shading refers to the effect of the 4/3 flux law convection parameterisation. The green curves emerge from assuming convection-based Nusselt number relationships Nu = 0.313 Ra 0.29 and Nu = 0.28 Ra 0.309 after Pandey and others (Reference Pandey, Verma and Mishra2014) and Silano and others (Reference Silano, Sreenivasan and Verzicco2010) respectively. Measurements of samples with predominating c-axis orientation along the coastline are marked by blue ×; bottom samples are denoted with red dots.
Disregarding two extreme values from Nakawo and Sinha (Reference Nakawo and Sinha1984), one finds an overall variability in the plate spacing of $\pm 20\%$
. Statistics of plate spacing measurements have been provided by Weeks and Hamilton (Reference Weeks and Hamilton1962) and Nakawo and Sinha (Reference Nakawo and Sinha1984). The standard deviation of the plate spacing obtained from thin sections was typically $20\hbox {--}30\%$
of the mean values. This suggests that much of the variability in Figure 9 may be related to the fact that always a certain range of plate spacing is observed.
Crystal alignment
Crystal alignment has been documented with plate spacing observations by Nakawo and Sinha (Reference Nakawo and Sinha1984) and Sinha and Zhan (Reference Sinha and Zhan1996). In the latter study the whole ice core showed an alignment of the c-axis with the direction of the coast, while in the former only bottom samples were affected. In Figure 9a all data points with crystal alignment are highlighted with blue crosses. Overall, the plate spacings for samples with crystal alignment appear to lie $\approx \! 10\%$
above the average of the data. Nakawo and Sinha (Reference Nakawo and Sinha1984) investigated the dependence of plate spacing on crystal alignment in some more detail, comparing the a 0 values for crystals within a thin section (their Fig. 8). For five of seven thin sections analysed they found that plate spacings of crystals aligned with the coastline were of 5–20% larger than those with little alignment. Jeffries and others (Reference Jeffries, Weeks, Shaw and Morris1993) investigated a large amount of thin sections and found slightly larger plate spacings for samples with lower standard deviation in the crystal orientation. Although not being statistically significant, also these results indicated an increase of the plate spacing with crystal alignment.
Crystal alignment along the coastline or current direction has been discussed in a couple of studies (Weeks and Gow, Reference Weeks and Gow1978; Langhorne and Robinson, Reference Langhorne and Robinson1986). A possible explanation for larger plate spacings can be based on theoretical and experimental studies (Langhorne and Robinson, Reference Langhorne and Robinson1986; Forth and Wheeler, Reference Forth and Wheeler1992; Huang and others, Reference Huang, Geng and Zhou1993). These authors showed that the interaction of a shear flow with the solute rejected from growing crystals implies a redistribution of solute that creates the largest supercooling on the upstream side of crystals/plate. As a consequence the growing crystals turn upstream. Such inclinement would imply a more effective vertical heat conduction from brine to the ice cell tips, and yield a reduced temperature gradient, corresponding to higher plate spacing predictions. As the same time also fluid flow would enhance the heat transfer between supercooled brine and the ice crystals, leading to the same effect. Last but not least, the fluid flow will also decrease ΔC in the tip layer which implies an increase in the plate spacing due to a decrease in the salinity gradient. Models and observational studies for other systems than sea water support both the effect of fluid flow (Forth and Wheeler, Reference Forth and Wheeler1992; Huang and others, Reference Huang, Geng and Zhou1993) and due to the temperature field (Xing and others, Reference Xing2015).
Presence of platelet ice
Many of the plate spacing observations from Jeffries and others (Reference Jeffries, Weeks, Shaw and Morris1993) for McMurdo Sound are for congelation ice that grows in the interstitials of platelets. It has been shown that this implies higher ice growth rates (Smith and others, Reference Smith, Langhorne, Frew, Vennell and Haskell2012). One may suspect that this happened for the ice studied by Jeffries and others (Reference Jeffries, Weeks, Shaw and Morris1993). Below 1.7 m thickness all cores contain platelets (Fig. 14 of that paper) and the ice growth seems to speed up slightly (Fig. 3 of that paper). While the overall data did not indicate a significant difference between plate spacings of congelation ice and interstitial congelation ice (Table 4 of that paper giving 0.71 ± 0.08 and 0.69 ±0.10 mm respectively), Figure 9 indicates the interesting aspect that plate spacings were constant or slightly increasing with the estimated growth velocity. This observation is consistent with the present predictions when considering another interesting observation from Smith and others (Reference Smith, Langhorne, Frew, Vennell and Haskell2012). These authors found not only faster growth in the presence of platelet ice, but also a larger ratio of growth velocity and temperature gradient in the ice (their Fig. 11). In terms of the present theory this would correspond to larger plate spacings at a given growth velocity. A shift to the upper bound of predictions in Fig. 9 is indeed what is observed. It is then worth a note that also the two data points from Paige (Reference Paige1966) stem from field work in McMurdo Sound and might have been influenced by the presence of platelet ice. They also appear near the upper limit of plate spacing predictions in Figure 9. While theoretical expectations and available data regarding the presence of platelet ice appear consistent, an alternative explanation for larger plate spacing could be the suppression of solutal convection during freeze-front propagation into a sub-ice platelet layer. More research on this aspect is needed.
Modelling of solutal convection
The dark shading in Figure 9 shows the bounds 0.13 < c nu < 0.17 in the 4/3-flux law for solutal convection. It is seen that, variation of the convective parametrisation in this range only affects the plate spacing prediction by $\pm \! 5\%$
. However, two aspects of the turbulent convection parameterisation are further considered as uncertain. While the Nusselt number scaling based on Eqn (44) and the 4/3-flux law, Eqn (45) is often used, its general applicability is in question (e.g. Kelley, Reference Kelley1990). For example, with Ra s based on D/V, our model predicts Rayleigh numbers <107 when considering ice growth rates V > 0.5 cm d−1. For such a regime a different flux law may be more reasonable. In addition, one might have to consider a Prandtl or Schmidt number dependence (Grossmann and Lohse, Reference Grossmann and Lohse2000).
To show the potential influence of such a modification we adopted slightly different flux laws obtained in two recent studies of the Nusselt–Rayleigh number relationship (Silano and others, Reference Silano, Sreenivasan and Verzicco2010; Pandey and others, Reference Pandey, Verma and Mishra2014). Both scaling laws have in common that they predict a larger Nusselt number at moderate Rayleigh numbers, and thereby lead to a prediction of smaller plate spacings, as shown by the dashed and dash-dotted curves in Figure 9.
Last but not least, there might be an influence of convection emerging from the interior of the ice on the boundary layer convection. In the experiments shown by Middleton and others (Reference Middleton, Thomas, Wit and Tison2016) both processes are evident but operating on different length and time scales. While theoretical modelling based on mushy layer theory may be applied to predict the onset of both convection modes (Worster and Wettlaufer, Reference Worster and Wettlaufer1997), more theoretical research based on real microstructure is needed.
Dependence on ageing
The plate spacing is a microstructure scale that is set at the freezing interface. Based on the data of very thick ice and its consistency with the predictions, it seems likely that the plate spacing under certain conditions is retained in old ice. However, it will be less easy to measure when the ice has become less saline and plate boundaries are less well defined. One may suspect that slow diffusion processes may first remove the smallest sub-grains, and that this aspect may introduce a bias towards overestimating the plate spacing for older ice. However, Shokr and Sinha (Reference Shokr and Sinha2015) have presented observations of the structure for different ice age and imaging conditions, showing that many of the original sub-grain boundaries are retained, while discrimination may require different techniques. Based on the limited number of studies it seems likely that (i) during the growth season plate spacings are well retained, (ii) for ice that has warmed considerably, even above its freezing temperature, much (perhaps not all) of the signature will be lost and (iii) the transition between these regimes will depend on thermal history.
In Figure 9 the data points from the bottom sections of ice cores are marked by red dots. At first glance these appear slightly lower than average. However, when considering all the points marked as blue to be influenced by crystal alignment, then the bottom observations appear as less exceptional for the remaining data points. One finally should note that bottom growth conditions often relate to decreasing ice growth at the end of the season. Considering that plate spacings need a certain distance or time to adjust to the new growth velocity, bottom values are likely to be related to somewhat higher than actual growth velocities. With regards to these processes, and the observational uncertainties, there is no evidence of significantly larger plate spacings at the bottom, at least for the present data sources considered.
Summary and conclusions
The plate spacing is probably the most fundamental microstructure scale of sea ice and has been measured by many investigators. However, the dependence of the plate spacing on growth velocity has remained an open question. A few investigators have attempted to derive this dependence based on experimental data, yet have come to rather different results (Assur and Weeks, Reference Assur and Weeks W1963; Lofgren and Weeks, Reference Lofgren and Weeks1969; Nakawo and Sinha, Reference Nakawo and Sinha1984; Maus, Reference Maus2007b). To address this problem in detail, this study presented compiled observations from published studies, and compared these to an approach to predict the plate spacing based on MST.
As summarised in Figure 10 the theory is suitable to predict the plate spacing a 0 of sea ice and marine ice over 5 orders of magnitude in the growth velocity. The model is based on a continuum approach in which all processes and parameters at a cellular interface are formulated in terms of the variable k, an effective solute distribution coefficient. This allows a broad investigation of the model's response to changing k and evaluation of its influence through parameterisations of heat and convective solute transport. The transport bounds, set by plausible parameter ranges, then yield lower and upper bounds for a 0. Based on the model results, largely supported by observations, one can distinguish three regimes:
(i) A high velocity regime (V > 20 cm d−1) where a 0 ∝ V −2/3. Solute diffusion is the dominant process that controls the plate spacing. Observations fall close to the upper bound of the present model. This is consistent with the fact that the data in this regime stem from laboratory studies and freezing of aqueous NaCl saline solutions in finite containers.
(ii) An intermediate velocity regime (0.1 < V < 20 cm d−1) where a 0 ∝ V −1/3 is approached. This is the growth regime most relevant for natural sea ice. Diffusion and solutal convection at moderate Rayleigh numbers control the plate spacing. The latter is expected to vary between upper and lower model bounds related to convection, fluid flow, crystal alignment and some unaccounted microscopic details of the ice–water interface.
(iii) A slow growth regime (V < 0.1 cm d−1) where a 0 ∝ V −1/3. This regime is relevant for marine ice shelves. The plate spacing is controlled by diffusion and solutal convection at high Rayleigh numbers. High solid fractions at the interface imply more constrained model predictions than in (ii). However, so far only a few observations exist.
Marginal stability wavelength λ mi predicted by MST, demonstrating the difference between a purely diffusive solution and one with solutal convection implemented in the approximation Eqn (31). For the legend of all data points see Figures 7 and 8. The shading shows the range of the predictions based on bounds for convection and the interfacial temperature gradient. Red curve: Numerical prediction for diffusion; blue curve: approximation for convection. The green lines, clearly matching the simulations, are least square fits to the diffusion model above V = 15 cm d−1 and to the convection model below this value, assuming slopes of −2/3 for diffusion and −1/3 for convection.
The model results can be presented in the form of log–log least square fits of a 0 versus V. A fit to the diffusive and convective model solutions gave exponents only slightly different from −2/3 and −1/3 respectively. Fixing the exponents to these values, the following best fits to the standard model were obtained:
where a 0 is in mm and V in cm d−1. Note that these approximations match at a growth velocity of 18.1 cm d−1, which is slightly different from the value deduced for V c at the onset of convection.
The parameterisations of heat and solute transport, leading to the present plate spacing predictions, are based on a simplified geometry of a cellular interface. More studies on the interface morphology, and associated heat and solute transport are needed for fine-tuning the model parameters, in particular at low growth velocities. As the standard deviations in the measurements of plate spacing were typically 20–30%, well controlled laboratory experiments followed by 3-D imaging of the interface, are one path to improvement. It is also likely that unpublished and published field datasets are already available for revision and reanalysis to be compared with the present theory. In particular for the low growth rates of ice shelves there is a need for more observations. If the present theory can be confirmed, this would open new possibilities to trace growth rates based on microstructure, both for marine ice shelves and sea ice.
Another path is detailed numerical simulations of crystal growth into sea water. A potential method to make progress may be the phase-field method (e.g. Steinbach, Reference Steinbach2013) for which some proof of concept studies for sea ice exist (Berti and others, Reference Berti, Fabrizio and Grandi2013; Moravetz and others, Reference Moravetz, Thoms and Kutschan2017). Progress may also achieved by computationally less expensive mesoscopic crystal growth models that require certain assumptions about the crystal growth seeds and morphology (e.g. Wongpan and others, Reference Wongpan, Langhorne, Dempsey, Hahn-Woernle and Sun2015). While such studies will likely lead to improvements in predictions, the present model is a fair first step towards fundamental modelling of sea-ice microstructure and its physical properties. An important application, first proposed by Tsurikov (Reference Tsurikov1965), is to use the plate spacing in a structural model of sea-ice salinity entrapment and evolution. The author showed that a simplified approach led to promising results for the prediction of ice salinity (Maus, Reference Maus and Leppäranta2008). Concise salinity prediction will depend not only on the plate spacing, but also on internal convection and microstructure evolution within the brine layers between the lamellae (e.g. Assur and Weeks, Reference Assur and Weeks W1963), suggesting the need for 3-D modelling and observations.
Acknowledgments
This project was funded through the Research Council of Norway (RCN) programme PETROMAKS2 grant 243812 (Microscale Interaction of Oil with Sea Ice for Detection and Environmental Risk Management in Sustainable Operations, MOSIDEO). I like to thank three anonymous reviewers for their input to improve this paper.
Appendix
Model properties of aqueous NaCl solution at a water salinity of S b = 35 g/kg used in this study adopted from (Maus, Reference Maus2007a):
Freezing point, T f = −2.105°C
Liquidus slope, m = dT f/dS b = −0.0619 K ‰−1
Pure ice density, ρ i = 916.97 kg m−3
Brine density, ρ b = 1027.1 kg m−3
Thermal conductivity of brine, K b = 0.555 W m−1 K−1
Thermal conductivity of ice, K i = 2.167 W m−1 K−1
Diffusion coefficient NaCl in water, D = 6.16 × 10−10 m2 s−1
Kinematic viscosity of brine, ν = 1.93 × 10−6 m2 s−1
Latent heat of fusion, L v = 3.031 × 108 J m3
Ice–water interfacial free energy γsl = 0.0299 J m2
Gibbs–Thomson parameter, Γ = (T f + 273.15)γsl/L v = 2.674 × 10−8 K m
Open access
