A model for entrainment of sediment into sea ice by aggregation between frazil-ice crystals and sediment grains

Abstract

Avertical numerical model has been developed that simulates tank experiments of sediment entrainment into sea ice. Physical processes considered were: turbulent vertical diffusion of heat, salt, sediment, frazil ice and their aggregates; differential growth of frazil-ice crystals; secondary nucleation of crystals; and aggregation between sediment and ice. The model approximated the real size distribution of frazil ice and sediment using five classes of each. Frazil crystals (25 μm to 1.5 cm) were modelled as discs with a constant thickness of their diameter. Each class had a constant rise velocity based on the density of ice and drag forces. Sediment grains (1–600 μm) were modelled as constant density spheres, with corresponding sinking velocities. The vertical diffusion was set constant for experiments based on calculated turbulent rms velocities and dissipation rates from current data. The balance between the rise/sinking velocities and the constant vertical diffusion is an important feature of the model. The efficiency of the modeled entrainment process was estimated through α, an aggregation factor. Values for α are in the range 〈0.0003, 0.1〉, but average values are often close to 0.01. Entrainment increases with increasing sediment concentration and turbulence of the water, and heat flux to the air.

1. Introduction

Sediment-laden sea ice was observed during the first journey across the Arctic Ocean (Nansen, 1906), and later expeditions have confirmed that such ice can be found in all parts of the ocean (Barnes and others, 1982; Nürnberg and others, 1994; Meese and others, 1997). Sediment in sea ice is often patchy in nature; concentrations are usually in the range 5–500 mg L⁻¹ and vary both horizontally and vertically at the cm scale. The entrainment processes appear to be governed by episodic events (Eicken and others, 2000).

Frazil-ice crystals are thin dendritic crystals that form in turbulent supercooled water (Martin, 1981). Granular ice is congealed frazil ice, and sea ice with incorporated sediment is often granular ice. Such sediment-laden ice is assumed to be a result of sediment entrainment during formation of frazil ice in open water, described as suspension freezing (Reimnitz and others, 1992). When frazil ice stays in suspension due to turbulent diffusion it may “scavenge” suspended sediment (Osterkamp and Gosink, 1984), i.e. collide and aggregate, and bring the sediment along to form the sediment-laden granular ice at the surface. If frazil-ice particles collide and aggregate with coarse objects at the bottom, or a lot of sediment in suspension, the frazil ice will sink and form anchor ice. Such anchor ice may rise with its aggregated sediment, or become buoyant enough to bring coarse material to the surface, if the ice grows in volume.

Experiments conducted in small tanks with a short duration have shown high levels of entrainment for processes that can be termed “suspension freezing” (Kempema and others, 1993; Reimnitz and others, 1993; Ackermann and others, 1994). The experiments had high levels of turbulence and heat fluxes, but these were not quantified. Recently, some experiments were completed in a larger tank at a scale more likely to represent natural conditions in time and space (Smedsrud, 1998, 2001; Haas and others, 1999).

›Different approaches have been applied to model the entrainment numerically. Two different vertical models have been developed, both showing the temporal nature of the process (Eidsvik, 1998; Sherwood, 2000). The efficiency of the aggregation was estimated using a simple box model based on classic aggregation theory, and the first set of laboratory experiments (Smedsrud, 1998). Frazil-ice growth and concentration in these models follow the simplifications in Omstedt (1985).

In this paper the model of Sherwood (2000) is developed further by including several important new processes. These include a size distribution for frazil ice and sediment, differential growth of crystals, and secondary nucleation. Aggregation between frazil ice and sediment is modelled using different size classes. Results from Smedsrud (2001) are used to tune the model parameters, and a size-dependent aggregation factor is calculated. Finally, the model sensitivity to its key parameters is tested, and predictions are made for different basic forcings to the model.

2. Model Development

The frazil and sediment model (Frasemo) represents time- dependent vertical profiles of vertical mixing, temperature, salinity and concentrations of sediment and ice. The model was developed by Sherwood (2000) and mainly used to calculate co-concentrations of sediment and ice for a given site in the Kara Sea, Russia, as a result of surface wind forcing and low air temperatures.

The model is here set to represent the 1 m deep tank in experiments A–D, described in Smedsrud (2001). Hereafter these experiments are referred to just by their letters (see table 1 in Smedsrud (2001) for a summary of experimental parameters). No attempt is made to model the horizontal flow, and the aim is to model the vertical diffusion and concentrations satisfactorily.

Table 1. Sediment size classes and sinking velocities

2.1. Diffusion and turbulence

The vertical diffusion has implications for the vertical gradients of heat, salinity, sediment and ice. Equations for the approximate conservation of these are:

(1)

(2)

(3)

(4)

Here T _w is the water temperature, S _w is the salinity, C _s is the total volume concentration of sediments, C _i is the total volume concentration of ice, and K is the vertical eddy diffusivity. G _T is the source of heat from the freezing crystals, G _s is the source of salt, and G _I is the source of frazil ice. w _i is the rise velocity of the crystals, and w _s is the sinking velocity of the grains. A constant uniform vertical eddy viscosity, K, is set through the depth. K was estimated for A–D using the mean rms fluctuation velocity, q, and the turbulent rate of dissipation, ε, as stated in Smedsrud (2001).

Equation (1) assumes that sediment C _s is at the same temperature as the local water, while the frazil ice C _i is assumed to be at the freezing point. C _s(j) in Equation (3) consists of (j = 5) size classes (Table 1), where the grain diameter (2r _s) and the sinking velocity (w _s(j)) are given. C _i(k) in Equation (4) consists of (k = 5) size classes, where a diameter d _i(k) and a rise velocity w _i(k) are required (Table 2).

Table 2. Frazil-ice size classes, C _i(k), and calculated rise velocities based on equations (1–3) in Gosink and Osterkamp (1983)

The vertical diffusion of T _w, S _w, C _i and C _s is modelled with a fully implicit diffusion scheme (Patankar, 1980). The routine is run separately for the different size classes of sediment and ice, as well as for T _w and S _w.

2.2. Model initialization and forcing

Each model run is started with both the measured T _w (close to the freezing point) and the measured S _w, vertically homogeneous.

The initial concentration of C _s is not straightforward, because the measured mass concentration at 0.5 m depth of, say, 18 mg L⁻¹ contains an unknown size distribution. The distribution is a balance between the upward diffusion and the sinking velocities of the variously sized grains. To solve this, a volume concentration, comparable to the total added mass of sediment, is used for the initialization, with the observed size distribution. The model is then run for 1 hour with diffusion, but no cooling, to allow the largest grains to sink, and avoid any aggregation with frazil ice before the steady-state balance is reached for the sediments. It is then controlled so that the total C _s matches the total measured mass concentration.

A very light snowfall is specified (1.0 × 10⁻⁸ ms⁻¹ ∼0.9 mm per 24 hours) to simulate the “snow” falling from the ceiling. This snow was observed to consist of crystals >2–3 mm in diameter, so it is specified as equally distributed over the three largest size classes. The initial concentration of ice crystals is zero for all classes, and increases after the first hour when the “snow” starts falling.

Air–sea heat flux is calculated by the model at the upper boundary, using standard bulk sensible-heat-flux parameterizations (Gill, 1982). Because the circulating water spent ∼65% of the time under the insulating ice cover, the calculated heat-flux coefficient c _H = 2.73 × 10⁻³ (Smedsrud, 2001) is reduced to the standard “open-sea” value of c _H = 1.1 × 10⁻³ (Simonsen and Haugan, 1996). This value is kept constant, and the observed air temperatures of the tank, from the corresponding experiments A–D, are used during each run.

Numerical solutions are obtained using finite differences on a staggered vertical grid with constant spacing of 5 cm and a time-step of 1 s. Apart from the specified flux of heat and snow at the surface, insulating boundaries are used for the other parameters at the top, and for all parameters at the bottom.

2.3. Sediment

For the sediment in suspension, a balance is reached between the upward diffusion and the downward advection of sediment grains, depending on their size and corresponding sinking velocity.

The size classes used in the model are shown in Table 1. All sizes are given as a grain diameter. These five classes approximate the real size distribution in A–C as shown in Smedsrud (2001, fig. 9). The sediment used in D consisted of clay and silts only, with a median diameter of 2.5 μm. The initial sediment volume was in this case set to 50% of C _s(1) and 50% of C _s(2). The sinking velocities are calculated from the Stokes settling velocities (Neilsen, 1992), using T _w = −2.0°C, S _w = 35 psu (practical salinity units) and a sediment density of 2650 kg m⁻³.

Assumptions about the form of the grains are necessary for modelling the aggregation, but are not needed to model sediment concentrations other than in these calculations. The “sphere assumption” is applied for the grain form, although many different shapes can be found. The sediment is assumed to be non-cohesive. That is, it does not form sediment–sediment aggregates.

2.4. Frazil ice

Frazil ice is the key factor in the entrainment process. It is also a very dynamic and complex feature to model. The volume concentration of the frazil ice increases steadily over time due to the heat flux, but also the size distribution changes as time passes. The frazil-ice size range in the model is similar to the range in other numerical studies and experiments (Gosink and Osterkamp, 1983; Hammar and Shen, 1995; Svensson and Omstedt, 1998). The frazil-ice size distribution is based on video images from D, and the constructed size groups are shown in Table 2.

The rise velocities are calculated using a constant crystal thickness of as suggested by measurements (Gosink and Osterkamp, 1983). A steady-state situation is reached in 0.3 s or shorter for the given size classes, so the rise velocities are treated as constants in the main model. The calculated rise velocities are given in Table 2. For the “Small” crystal class the given value is obtained by linear interpolation between zero and the “Fine”class.

As Frasemo uses w _i and K to calculate vertical gradients of C _i with Equation (4), no assumption needs to be made about the crystal form. However, when collision between crystals, and aggregation between frazil-ice and sediment grains, are modelled, the crystals are assumed to be spheres. When this is done, an equivalent radius r _ie is used, instead of the given d _i. r _ie is calculated so that a disc of diameter d _i and thickness t _i has the same volume as the sphere of radius r _ie. Values for r _ie are also given in Table 2.

2.4.1. Differential growth

The heat flux from the ice crystal of class (k) to the surrounding water is described by

(5)

Here N_u is a Nusselt number describing the ratio between the actual(turbulent) heat flux and the heat conduction. N_u may vary with the flow conditions, and has earlier been set on the order of 1 (Svensson and Omstedt, 1994), or higher when N_u is set dependent on the turbulent dissipation rate (ε) and the Kolmogorov length scale (λ) (Hammar and Shen, 1995). The thermal conductivity of sea water is set to K _w = 0.564 W m⁻¹ °C⁻¹ (Caldwell, 1974).

Equation (5) is comparable with Equation (3) in Svensson and Omstedt (1994), transferred from number of crystals per volume of water to volume of crystals per volume of water. The radius is chosen as the characteristic length scale following Hammar and Shen (1995).

The models of Svensson and Omstedt (1994) and Hammar and Shen (1995) were validated by fresh-water laboratory experiments like that of Carstens (1966). Svensson and Omstedt (1994) also used direct observations of particle size and total number of crystals. It seems that there must be a difference in the residual supercooling between fresh and salt water, because Carstens (1966) reports only “a few thousand” after the first few minutes, while A–D show “a couple of hundred” during 24 hours.

The salt-water experiment of Tsang and Hanley (1985) may be compared to A–D, but there seeding took place at different levels of supercooling. The residual supercooling was not mentioned explicitly, but the reported salinity (S _w = 29–30 psu) and residual T _w gives a residual supercooling of 0.08–0.2°C after the first frazil formation took place. The longest experiment lasted about 15 min, and had a residual supercooling of 0.141°C (S _w = 30 psu), comparable with A–D. Omstedt (1985) examined these experiments and found that the data could be represented using N_u = 4.0, d _i = 1 mm and .

The models mentioned here describe the situation from a state of supercooling, through a seeding process, until “equilibrium” is reached. While the seeding gives an accurate starting point for the ice-growth process in the calculation, it does not represent most natural conditions where some snow will fall into the water before, as well as when, the water reaches its freezing point. Here an attempt is made to model the situation from a start with a few crystals at the surface, through a maximum supercooling, and for a long period of time during the subsequent “equilibrium” as well.

The crystal is assumed only to grow at the edges, so the active freezing area is πd _i(k)t _i(k) in Equation (5). , and the heat flux from a crystal is therefore independent of the thickness. The total heating of the surrounding water from the growth of each crystal class in Equation (1) can then be calculated as

(6)

Here ρ _w is the density of water, ρ _i is the density of ice, C _pw = 3989 J kg⁻¹ °C⁻¹, the heat capacity of water, and is the total ice concentration in the gridpoint. The total heating in a gridpoint in Equation (1) is then . This heat flux results directly in freezing calculated with

(7)

ΔC _i(k) is calculated for each size class (k = 1, 4), and represents the growth of a number of ice crystals. Notice that the “Large” crystals (k = 5) are not permitted to grow, because they have already reached their maximum size. Because of the constant mean diameter in the size classes (d _i(k)) this growth has to be transformed to a certain volume (or, alternatively, number) of crystals being transferred to the next size class (Hammar and Shen, 1995):

(8)

Here v _i(k) is the volume of the ice crystals in that specific class, and Δv _i(k) = v _i(k + 1) − v _i(k). With this formulation, the “Large” class increases in volume because of growth in the “Coarse” class, while the “Smalf”crystals (k = 1) always have a pure loss to the “Fine” class due to the growth. Total growth of frazil ice is then ,and the growth results in a salt flux calculated as Gs = S _w G _I in Equation (2).

2.4.2. Secondary nucleation

Secondary nucleation is the term used for production ofnew small crystals by removal of nuclei from the surface of parent crystals. The main processes thought to occur are collision between crystals and a resulting small piece of ice (collision breeding), and detachment of surface irregularities by fluid shear (Daly, 1984). The simplified approach of Svensson and Omstedt (1994) is followed here. A crystal in relative movement to the fluid will sweep a volume ΔV _i during a time interval Δt:

(9)

where

(10)

incorporating both the rise velocity and the turbulence intensity. r _ie, the equivalent radius, has to be used here, because the crystal can twist and turn in all directions, and no way has been found of modelling a disc in a turbulent flow. The increase in volume for the “Small” size class from collision between all the different size classes is then calculated as:

(11)

is the average number of all the different ice crystals in the grid volume. ΔC_i (k = 1) in Equation (11) is always positive, and there is a corresponding loss of (exactly the same) volume for the other classes for each part of the summation in Equation (11). is a calibration constant, and a maximum value has been set for each model run, limiting the secondary nucleation process. Values found are generally much lower than the value of 4.0 × 10⁶ given in Svensson and Omstedt (1994), and a further discussion of n_i is given in section 4.2.

2.5. Aggregation

All substances suspended in water have the potential to collide and aggregate with others. Suspended frazil ice tends to form flocs (crystal aggregates) in fresh-water experiments, but this was not reported in Hanley and Tsang (1984) or A–D. Frazil will also collide and aggregate with sediment to some extent, and this occurred during experiments A–D.

The term aggregation is preferred as the general mechanism, instead of the terms coagulation or flocculation, which should refer to the process that makes two substances aggregate (Shamlou, 1993). Coagulation is the preferred term when two substances are not in physical contact but separated by a thin liquid film, kept together by London–Van der Waals forces. Flocculation is when there is physical contact between the two particles: ice-crystal to ice-crystal aggregation is thought to be primarily of this kind.

Aggregation here also includes what has been termed “scavenging” of sediment by frazil ice in suspension, or “filtration” of sediment in the surface slush (Osterkamp and Gosink, 1984), as well as “mechanical trapping”, or whatever process makes the sediment and frazil stick together when they collide.

The basic relationship for the number of collisions per unit time, J _is, between frazil ice and sediment is given by Smoluchowski (1917):

(12)

Here, frazil ice has the radius r _ie and sediment r _s, and their number concentrations are n _i and n _s, respectively.

Ω is a diffusive coefficient, dependent on the flow. For conditions where there is no velocity, Ω is the Brownian motion (molecular diffusion). For turbulent conditions, and if the colliding particles are smaller than about 10λ (Shamlou, 1993),

(13)

λ = (v ³/∊)^1/4 ∼ 350 μm, using the turbulent dissipation rate ∊ from A–C, and a constant kinematic viscosity v = 1.8 × 10⁻⁶ , and only the largest crystals approach 10λ. In Equation (13), r is the size of the eddies.

By setting the size of the eddies, r, equal to (r _ie + r _s), the “particle diameter”, an expression for the aggregation frequency, can be obtained. The increase in number of new aggregates n_x can then be expressed as:

(14)

Here t _T isTaylor’s time-scale, representing turbulent strength, t _T = (15v/∊)^1/2 (Tennekes and Lumley, 1994). The aggregation factor, α, describes the statistical chance of aggregation between the frazil and the sediment, and incorporates a colliding efficiency. This colliding efficiency describes the statistical chance of a collision between the two, and is assumed to be constant since collisions involve soft and non-uniform aggregates. Only a fraction of the collisions will lead to an aggregate, and values for α are expected to be less than unity, with zero as the lower boundary. The aggregation factor α is similar, but not identical, to the collection efficiency E (Osterkamp and Gosink, 1984), and also associated with the “frazil stickiness”“discussed by Kempema and others (1993).

The theoretical framework has been tested on mono-sized polystyrene and spherical latex particles, and fits well to experimental data (Higashitani and others, 1983; Gierczycki and Shamlou, 1996). It is expected that every substance has its own aggregation factor α (and also colliding efficiency), and α should therefore be looked upon more as an empirical constant.

Here two different types of aggregates are classified, with their own size distribution. Ice that has one or more sediment grains aggregated to the crystal will be termed hybrid ice C _is(k). C _is has the same sizes as C _i(k). Likewise, all the sediment that has aggregated to ice will be termed hybrid sediment C _si(j), with the same size as C _s(j).

It is now assumed that the basic collision Equation (12) describes the situation even if some of the particles are not floating alone. The sediment is assumed non-cohesive, and frazil flocs are not very abundant.

Looking at the loss of sediment in the “clay”“ class (j = 1), Equation (14) becomes;

(15)

This is identical to equation (3) in Smedsrud (1998), adding the five size classes. The transfer to volumes is achieved by multiplying with the volume of the respective (spherical) particles, . When a clay grain collides with one of the ice crystals, or one of the hybrid ice crystals that has already aggregated to a sediment grain, aggregation may occur. The sediment volume then has a statistical chance, α, of being transferred to C _si. The loss and gain are exactly the same volume, so ΔC _si (j = 1) = −ΔC _S(j = 1). Hereafter IRS refers to measurements of ice-rafted sediment, i.e. sediment that remained on drained frazil crystals or was sampled within solid ice in A–D, while C _si refers to model results.

In the same way, the equation for loss in an ice class (e.g. k = 5), becomes:

(16)

which describes exactly the same processes. C _is and C _si are always increasing in volume.

There are no such terms as C _is [C _s + C _si] as one might expect, because even if a hybrid ice crystal aggregates with a sediment grain, the volume concentration of C _is does not increase. The ice crystal is already “tagged” as hybrid. Likewise, a hybrid sediment already on an ice crystal does not change its type if it aggregates with yet another ice crystal, i. e. there are no terms like C _si[C _i + C _is] either.

The aggregation factor a is now the only parameter not measured in A–D, and the model is used to estimate it. From the box-model approach in Smedsrud (1998) α ∼0.025, but it might vary between the different size classes.

The results from the box model (Smedsrud, 1998, fig. 2) showed that Equations (15) and (16) give rise to exponentiallike solutions of the increasing C _si. This is caused by the product of the different volume concentrations in Equations (15) and (16). This fits qualitatively with the IRS data presented in Smedsrud (2001).

C _is and C _si are modelled separately by their volumes, but in reality they have formed several combinations of hybrid particles, with an unknown combination of the 10 size classes. These hybrids will now have a different rise or sinking velocity, depending on how much of each type is in a specific gridpoint. The average density of the hybrid particles at any level is calculated using:

(17)

Here ρ _i = 920 kg m⁻³, the density of pure ice, and ρ _s = 2650 kg m⁻³, the density of the sediment. ρ _h is used to calculate the rise or sinking velocity by linear interpolation. For the hybrid ice C _is,

(18)

and for the hybrid sediment C _si,

(19)

is the time-averaged water density over the length of the experiment. This makes the range of w _is and w _si equal to the range of w _i and w _s. The hybrid aggregates will tend to sink if the volume of C _si in a gridpoint is larger than 6% of C _is. Then anchor ice will be produced, i. e. frazil ice kept close to the bottom by sediment or larger pieces of material. Anchor ice may grow due to adherence of additional crystals, or from heat transfer to the supercooled water (Martin, 1981). This may happen with the hybrid ice and sediment that has sunk in the model as well, and this ice/sediment may rise again if enough frazil aggregates or grows so that the hybrid ice volume increases significantly. From Equation (17) the theoretical maximum concentration of IRS in neutrally buoyant surface slush using is found to be 159 g L⁻¹. This is much higher than IRS observations from A–D, or from the field, which are usually <1g L⁻¹ (Nürnberg and others, 1994).

3. Model Tuning

The first hour of each run has no snow and no cooling (T _a = −2.0°C), so the sediment is allowed to reach a steady-state vertical concentration. Then snow is added at the given rate for 1 hour, equally distributed over the three largest crystal sizes. After 2 hours the air temperature is decreased linearly over 1 hour to the observed mean of the given experiment (Ta = 〈−10.43, −17.17〉°C).

3.1. Initial sediment concentrations

The steady concentration of sediment in the model is a balance between the vertical diffusion K and the sinking velocities w _s. The size distribution is based on observations and does not change between experiments A–C, but the mass of added sediment is changed in the model to match the measured initial concentration. Hereafter the term suspended particulate matter (SPM) is used when referring to measurements, and C _s when referring to model results. As B will be used to calculate aggregation factors in section 4, the initial C _s of that experiment is discussed here as well.

C _s is fairly homogeneous in the upper 0.5 m, and the value at 0.5 m matches the observed SPM of 18 mg L⁻¹ in the beginning of B (Fig. 1). As the total volume of sediment in the tank remains constant, the steady-state situation will remain until the aggregation starts to be effective, and a significant volume of C _s is transferred to the hybrid class C _si. This decreases the sinking velocity, and depending on the volume of the hybrid ice Cis, some sediment will be diffused upwards by the turbulence from the situation shown in Figure 1. The added SPM corresponds to a homogeneous vertical concentration of 40 mg L⁻¹, or a total dry weight of 5.16 kg. The total added mass of sediment was 7.8 kg, including about 35% water.

Fig. 1. Modelled vertical volume concentration of sediment, C _s, as cooling starts in experiment B.

The SPM used in D consisted of C _s (1) and C _s (2) only. Due to their small w _s and the stronger vertical diffusion, initial C _s was homogeneous with depth, at 12.0 mg L⁻¹, equally divided between the two classes.

3.2. Supercooling

Experiment D is used to tune Frasemo to the observed supercooling and frazil size distribution. The supercooling drives the formation of frazil ice as described by Equation (5), and thereby also influences the frazil size distribution. Experiments A–C have salinity and temperature records, and are used as a control of the modelled supercooling. The given equations and forcing leave two parameters to tune, the Nusselt number N_u in Equation (5), and the in Equation (11).

The choice of N_u has minimal impact on the total C _i, and the measured supercooling becomes insignificant at 24 hours (∼0.002°C). Observations show that T _w is ∼0.02 below T _f 2–5 hours into experiments, and that T _w increases to ∼0.01 below T _f after 6–20 hours. This is achieved using N_u ∼1.

This is somewhat surprising since it indicates that the turbulence does not increase the crystal growth, and it is close to a situation with conduction only. A constant Nusselt number implies that the crystals are large in comparison to the turbulent dissipation length scale. This is discussed further in section 4. Avalue of N_u = 1.5 provides results that fit A–D quite well.

Figure 2 shows the time evolution of T _w and S _w, and the corresponding freezing point T _f (same data as shown in figure 1 of Smedsrud, 1998). The modelled salinity differs < 0.05 psu from the observations (not shown here), which makes T _f differ < 0.003°C. The calculated T _w follows observations to within 0.01°C as seen in Figure 2. The decrease of T _w to a minimum temperature, and then a return to an “equilibrium” temperature slightly less than T _f is a well-documented feature and is found in all frazil-ice experiments (Carstens, 1966; Daly, 1984).

Fig. 2. Temperature T _w (solid line) salinity S _w (dot-dashed line) and the corresponding freezing point T _f (dashed line) at 0.5 m depth as calculated by the model for experiment D. The measured T _w (jagged line) is also plotted for comparison.

For the 24 hour experiments A–C, model results are qualitatively the same. Supercooling of about 0.02–0.04°C persists for up to 10 hours, and then slowly decreases to < 0.01°C at 24 hours.

3.3. Frazil-ice volumes

Volumes of frazil ice in suspension were measured during all experiments A–D, but most measurements were taken during C and D, so they will be used here. Hereafter the term C _i is used for modelled values of frazil ice, and “frazil ice” is used for measured values. C _i at different depths is dependent on two parameters, the rise velocity, w _i(k), and the (constant) vertical turbulent diffusion, K.

w _i(k) depends on size as well as thickness, and a mean crystal-thickness to -diameter ratio of given by Gosink and Osterkamp (1983) is used. In A–D the thickness was not measured directly, only qualitatively estimated as t _i < 1 mm, for the crystals with di ∼20 mm. Values for K are estimated using the rms fluctuation velocity q, and the turbulent dissipation rate ∊ (equation (5) in Smedsrud, 2001). In C K = 5.9 × 10⁻³, using the mean q = 10.3 cm s⁻¹, and .

In D K = 111.4 × 10⁻³, using the mean q = 6.8 cm s⁻¹, and . The vertical profiles of C _i(k) in D show that the two smallest classes have concentrations of <1 mg L⁻¹ (Fig. 3). C _i (3) and C _i (4) are close to vertically homogeneous, and the major part of the volume is found as C _i (5). The total calculated volume is close to the measured range of frazil ice 3.0−4.3gL⁻¹ between 0.1 and 1 m depth during D (Smedsrud, 1998, fig. 1). At 6 hours C _i has a surface concentration of 5 g L⁻¹,decreasing to 3.7 g L⁻¹ at 1.0 m depth.

Fig. 3. Modelled vertical profile of frazil ice in suspension, C _i, 5 hours into experiment D.

In C, more of C _i is close to the surface due to the lower vertical diffusion. Surface concentrations reach 23 g L⁻¹ in 10 hours, and 58 g L⁻¹ in 20 hours. Frazil-ice measurements indicate a linear increase at 0.5 m depth to around 1g L⁻¹ at 10 hours. This is also calculated by the model. Frazil-ice values from 0.25 and 0.75 m depth are also close to 1g L⁻¹, but C _i has a larger gradient, giving values from 0.7 g L⁻¹ at 0.75 m to 7 g L⁻¹ at 0.25 m.

3.4. Frazil-ice size

The relative concentrations of frazil ice in the surface were counted from video images filmed during D. Countings resulted in a close to constant distribution in time, where ∼60% of the ice belonged to the Large class, ∼37% to the Coarse class and ∼3% to the Medium class.Volumes of the two smallest classes are insignificant, both in the model and in the observations. This is consistent with the thermodynamic stability of frazil crystals, whereby the crystals have a maximum diameter given by their thickness t _i, the level of supercooling and a number concentration (Forest, 1986). For a supercooling of 0.005°C and a number concentration of 10⁷m³, the maximum diameter of 2 cm corresponds to a thickness of ∼47 μm. This process limits the growth of frazil-ice crystals over long time-scales.

At 1.25 hours into D, frazil-ice sizes were 52.3% for the Large class (k = 5), 43.5% for the Coarse class (k = 4), and 4.1% for the Medium class (k = 3).With time, the volume of the Large class increased most efficiently, and at 6 hours the relative concentration was 63.7%.

The size distribution predicted by the model generally differs <5% from the observations. Towards the end of the experiment the model predicts too high an increase for C _i(5), and corresponding lower values for C _i(4) and C _i(3). The relative concentration of C _i(5) at 5 hours at the surface in Figure 3 is 69.4%. For A–C, C _i(5) dominates, with ∼85% of the volume at the surface, while C _i(4) holds about 14% of the volume. Less than 1% is left for C _i(3), and again the two smallest sizes are insignificant.

3.5. Size-dependent aggregation

Observations of IRS show high variation, and the focus here is therefore to model the average values of the aggregation process; the variation and sensitivity will be discussed in the next section. Samples of SPM and IRS were taken in all experiments A–D, but the size distribution of the IRS was only measured in B. B is therefore used to calculate the size-dependent aggregation factors, and A, C and D are used to discuss the model sensitivity.

The increased chance of collision between larger particles is accounted for by the different radii in Equation (15). This means that the chance of aggregation for each sediment size class reflected by α is dependent on more than size, for instance shape of the grains.

Because of the large variations of the IRS concentrations, α(j) is adjusted to give the average IRS of 100 mg L⁻¹ with the measured relative concentrations. Some of the measured IRS concentrations are >100 mg L⁻¹. These are probably local maxima from the horizontal variations. The measured SPM of 7.6 mg L⁻¹ at the end of B provides a control for the mass of sediment entrained into the surface slush and “lost” from the water.

No grains larger than ∼150 μm were found in the IRS, so C _si (5) = 0. The four remaining C _si size classes show similar vertical gradients (Fig. 4). This is because they are aggregated to ice and have a hybrid rise/sinking velocity w _si depending on C _is and C _si. In A–D C _is was usually 1–100 g L⁻¹, much higher than C _si, in the range 1–100 mg L⁻¹. This results in ρ _h close to the ice density, calculated by Equation (17), and corresponding rise velocities calculated by Equation (19). w _si in B is 13.0 mm s⁻¹ at the surface, and decreases to 6.3 mm s⁻¹ at the bottom after 23 hours.

Fig. 4. Vertical profiles of sediment aggregated to frazil ice at 23 hours in experiment B.

The relative concentrations of the different grain-sizes at the surface are close to the observed values. The major mass is C _si(4) (66.8%). There is 17.5% of C _si(3), and 18.3% of C _si(2) and only 0.4% of C _si(1). The differences between the calculated size distribution and the measured size distribution of the surface IRS (Smedsrud, 2001, fig. 17), are < 1 mg L⁻¹.

As α(j) is treated as an empirical constant in Equation (15), any concentration of the different classes can be reproduced, up to the level where all SPM has transferred to IRS. The estimated values are α(1, 2, 3, 4) = (0.2, 5.0, 0.85, 16.5) × 10⁻³. α(5) is set to zero to reproduce the observations of no grains larger than ∼150 μm in the IRS.

Values for α(j) are found within a fairly compact range. It is therefore tempting to include a constant a case. An a of 0.01 gives a surface C _si of 103 mg L⁻¹ and the same vertical gradients as in Figure 4. However, the C _si size distribution shows a significant trend towards smaller grains, and the values for the surface are C _si (1, 2, 3, 4, 5) = (2.5, 17.1, 35.7, 43.4, 1.3)%. Notice that only 1.3% of C _s(5) ends up as IRS, despite the large size of these coarse and medium sand grains. This indicates that collisions between C _s(5) and ice occur quite rarely. This must be due to the low value of C _s(5) close to the surface (Fig. 1) where C _i has the largest volume, as well as the small C _i volumes close to the bottom of the tank where C _s(5) has high concentrations. Vertical distribution of modelled C _i in B and C is very similar, so surface C _i are roughly 10 times higher than bottom values. For the three smallest sediment classes, almost nothing remains as C _s(1–3) in the constant α case, and they are mostly transferred to hybrid sediments.

4. Model Sensitivity

Much of the forcing and initialization in the model is robust and based on measurements where the physical process is well understood, and the variations are known to be within a given range (T _w, S _w, T _a, q and C _s). Other settings (∊, K and C _i) are partly based on measurements, where values are within the observed range, or measurements were qualitative only, or showed high variations. A few settings (N_u , and α) are tuned constants, values chosen to make the model predictions comparable to observations. The sensitivity to the less robust parameters is discussed in this section.

4.1. Turbulence

The vertical turbulent diffusion, K, controls the vertical gradients of C _s and C _i as described by Equations (3) and (4). If K is very low, sediments are found close to the bottom, ice stays at the surface, and no aggregation can occur. Turbulence also affects the collision frequency, through t _T in Equation (15), and the secondary nucleation of frazil crystals through U _r in Equation (10).

K is estimated with a common parameterization assuming a homogeneous vertical density (Mellor and Yamada, 1982), using q and ∊ (Smedsrud, 2001). Measurements of q are robust, but the 95% confidence level of the ∊ estimate covers an order of magnitude (Smedsrud, 2001). The overall range of K calculated from the current data of A–D was 〈1.7,111.4〉×10⁻³. There is spatial and time variation of the turbulent parameters in the tank. For instance, the increasing ice concentration at the surface affects the flow, through both wave dampening by the slush, and increase in viscosity.

A density interface also develops at the lower side of the slush, while the bed load of sediment has a minor effect on the density. The assumption of a homogeneous water column, and thereby the constant K with depth also, is supported by the fact that it is possible to reproduce measurements.

C _i shown in Figure 3 is not very sensitive to K. Increasing K by 100% to 222.8 × 10⁻³ decreases the vertical gradient, and C _i becomes almost vertically homogeneous. Decreasing K to 55.7 × 10⁻³ increases the gradient, and now a surface concentration of 4.3 g L is reached in 5 hours with 2.5 g L⁻¹ at the bottom.

4.2. Frazil ice

The volume and size of the snow are qualitative observations, and have to be considered. Completely removing the snow causes the supercooling to reach 0.45°C in 6 hours, and still no frazil ice is produced. Reducing the snow from 0.9 to 0.09 mm d⁻¹, produces a maximum supercooling of 0.12°C, much higher than observations. A ten-fold increase of the snow to 9 mm d⁻¹ makes the maximum supercooling 0.01°C, which is less than observed. This volume of snow starts to be significant compared to C _i; in 6 hours the snow reaches 11.5 g L⁻¹ if it is distributed over the upper 0.2 m.

If the size of the snow is changed to include the smallest ice classes, then the supercooling is absent at 0.5 m, and only a very small supercooling is present in the upper 0.2 m. The C _i size distribution is altered quite drastically to ∼45% of C _i(2) and C _i(3), and vertical concentrations become homogeneous due to the small rise velocity of these ice classes.

The model sensitivity to volume and size of snow thus seems to be quite large. The measured time development of T _w could not be reproduced using small sizes of snow, even when tuning the two other frazil coefficients N_u and . The volume of C _i does not change significantly with the variations of snow, except for the added volumes themselves.

The delicate balance between the snow size and the supercooling indicates the tuning of the other two parameters controlling the ice growth, N_u and . N_u controls the heat conduction from each crystal, and the production of small crystals, which grow faster than the larger ones. As for the snow, N_u and does not alter the volume of C _i significantly, only the size distribution and the supercooling.

Size- and turbulence-dependent N_u ’s are not included in the model, but may be calculated using equations in Hammar and Shen (1995). For ∊ = 3.7 × 10 from D, the Kolmogorov length scale becomes A = 1120 μm. This results in N_u (1, 2, 3, 4, 5) = (239.7, 24.5, 4.1, 2.0, 1.3). These N_u would further increase the more efficient growth of the smaller ice crystals and make it more difficult to reproduce the measured development of T _w. To validate the theoretical values of N_u properly, in situ size distribution over time would be appropriate, and this has not been accomplished to date.

Increasing N_u in Equation (5) to 3.0 reduces T _w to −2.005°C in Figure 2. Further increase of N_u to 6.0 reduces T _w to −1.995 °C. Neither size distribution nor vertical gradients are altered significantly by the change in N_u . Setting N_u = 1.0, the lower limit, reduces T _w to −2.04°C, a supercooling of 0.02°C too much.

The value of in Equation (11) represents the average value for the number of crystals in a gridpoint. Number concentrations of crystals are quite high, with C _i(5) = 10g L⁻¹, , or 123.0 × 10⁶ m⁻³ The process of secondary nucleation is not well understood, and many of the new crystals that are created by collisions between the larger ones (with highest size, inertia and rise velocities) might be smaller than the critical size at the present level of supercooling (Daly, 1994). These new, very small crystals will then dissolve, and not start to grow. The ones that are larger than the critical radius will start to grow. As the level of supercooling decreases, the critical radius increases, and thus there is less chance of producing new nuclei that are large enough. This limits the process, at the same time as increases, and provides a physical explanation of why an upper limit on is necessary.

The maximum value of that gave satisfactory results was 1.0 × 10³ m⁻³. No significant changes appear with a 10-fold decrease or increase of this value. Setting = 1.0 × 10⁶, of the same order as proposed by Svensson and Omstedt (1994), leads to quite a large difference. The minimum T _w reaches only −2.0°C, and returns to T _f before 2 hours is past, and follows T _f closely. The size distribution changes significantly, comparable to including small ice particles in the snow, with ∼45% in C _i(3) and C _i(4).

The surface frazil size distribution is based on counting the individual frazil crystals on a monitor during D, and is qualitatively observed in the other experiments. While the range is robust, the individual concentrations in the size groups must be considered as estimates, but when both the size distribution and the T _w development are altered, this indicates consistency in the measurements, and demands a tight forcing of the model.

Calculation of w _i is done using one set of experiments (Gosink and Osterkamp, 1983). Alternative equations exist, which generally give higher rise velocities. The velocity range 0.01–80 mm s⁻¹ is suggested by Daly (1984) for the size range given in Table 2. By setting , the thickest crystals in the range of Gosink and Osterkamp (1983), w _i(5) = 73.9 mms⁻¹, close to the suggested value from Daly (1984). This gives 3.2 g L⁻¹ of frazil at the bottom, 4.0 g L⁻¹ at 0.5 m depth, and 5.3 g L⁻¹ at the surface, again compared to Figure 3.

4.3. Aggregation

Volumes and size of C _s are considered robust, which provides a control for the value of K, because the total C _s at a given depth is known, and the model also “predicts” this volume concentration. The size distribution of C _si can only be considered an estimate, but the range is robust.

The calculation of α is one of the main results in this work. Smedsrud (1998) showed a to be most sensitive to the dissipation rate, ∊, and less sensitive to particle sizes. As a spectrum of sizes is now used, and the ranges are considered robust, sensitivity studies of size are not needed. The sensitivity of the constant α is discussed, as a size-dependent α can always be calculated, and the calculated C _si did not deviate too much from the measured IRS values using a constant α.

As the aggregation factor is empirical, the sensitivity of a should be tested against the other parameters of the model. The difference in IRS concentrations between A–D should also be discussed in this context. To do this, the entrainment factor, X _E = IRS/SPM, is utilized, as discussed in Smedsrud (2001).

Measured IRS, or calculated C _si, is compared with the initial SPM, or C _s, in each experiment. In this way the different experiments, A–D, can be easily compared. X _E has values of <1 when the surface C _si is smaller than the initial Cs at 0.5 m. With time the IRS concentration increased in the tank experiments (Fig. 5). The calculated surface C_si in B at 23 hours, 100 mg L⁻¹, results in X _E = 5.6.

Fig. 5. Normalized concentrations of sediment with time from experiment A–D (X_e = IRS/SPM). The different X_e ’s from the model are shown as (dashed line), andfrazilslush and granular IRS refers to whether the IRS was sampled in frazil ice, surface slush or surface granular ice.

For a given time evolution with associated concentrations and size of C _i and C _s,the relationship α/t _T from Equation (15) has to be constant to give a certain volume of C _si in the surface. This indicates the connection between a and the strength of the turbulence. For experiment B, α/t _T = 0.037. Increasing and decreasing e by one order of magnitude results in α = 〈0.003, 0.032〉, which corresponds approximately to the observed range of ∊ in A–C (Smedsrud, 2001, fig. 2). Setting α = 0.03 gives X _E = 9.1, and setting α = 0.002 gives X _E = 3.9, both after 23 hours. As seen in Figure 5, both these values are within the measured variation.

The maximum IRS concentration (200 mg L⁻¹ in B, X _E = 11.1) indicates an empirical maximum for the aggregation process in the tank experiments. Using α = 0.1 gives C _si = 202 mg L⁻¹, and the time evolution of X _E is the upper curve in Figure 5. Now all the sediment in classes C _s(1, 2, 3) is found only as C _si, and the bed load of the two remaining classes has only 20 mg L⁻¹. The size distribution improves from the other constant a runs, but only 0.5 mg L⁻¹ of C _s is left at 0.5 m, much lower than the observed 7.7 mg L⁻¹. The abrupt change in X _E at 4 hours occurs because there is no more Cs in the near-surface layer; it has already transformed to C _si. The empirical minimum value of α is found to be 0.0003 (Fig. 5), so the overall range is α = 〈0.0003, 0.1〉 based on A–D.

Experiments A and C have been modelled in the same way as B. The same parameters are used, except the measured initial (homogeneous) values of T _w and S _w, the measured average T _a, the measured average q, and the measured initial SPM (C _s). A gives an X _E of 3.35 after 23 hours using α = 0.01, and is also plotted in Figure 5. This is equal to an IRS concentration of 42 mg L⁻¹, and the measured range was 13–44 mg L⁻¹. In experiment C an X _E of 2.8 is predicted (43 mg L⁻¹), higher than the observed range of 6–33 mg L ¹.

Experiment D is different because of the high K, the smaller size of the sediment, and the short duration of 6 hours. At the end of the experiment the impellers were switched off, and all frazil crystals in suspension rose to the surface with their aggregated sediment. The IRS sampling was done after this. At 6 hours the total calculated C _is over the 1 m depth is added and divided over the two upper grid-points. Using α = 0.06, a surface C _si of 12.6 mg L⁻¹ is reached, the average of the five measurements. This is consistent with Figure 5, indicating that X _E ∼1 after 6 hours.

While the model cannot reproduce the exact concentration of IRS in each experiment, it does reproduce the main features, such as the highest IRS levels in B, and lower levels in A and C. The higher IRS values in B are a product of the high initial SPM and the high heat flux creating more frazil ice.

Why then does C have lower surface IRS than A when both the initial SPM and the heat flux are larger? The explanation is that there is more C _si in C, but less at the surface due to the larger K.

Waves are incorporated into the model in the way that the measured higher turbulent rms velocities increase K. In this way waves can increase turbulent diffusion of frazil ice downwards, and sediment upwards. But there is no mechanism that increases aggregation or entrainment by waves.

Aggregation seems to occur mainly in the upper half of the 1 m deep water column in A–C due to the relatively high frazil concentrations there. Towards the bottom, less frazil ice is found, as well as the larger sediment grains. Even with a constant a for all classes of C _s (j) almost no C _si(5) was created.

With the larger K and smaller sediment in D, the aggregation process takes place efficiently at all depths. Under such conditions the “box-model approach” (Smedsrud, 1998) seems appropriate.

5. Model Predictions

To further increase the understanding of the aggregation process, the model predictions are tested against three of the central and robust parameters. These are the constant air temperature, T _a, driving the formation of frazil ice, the initial sediment concentration C _s (1, 2, 3, 4, 5) and the turbulent rms velocity, q. B is used to test these predictions, and the corresponding size-dependent α’s.

The mean T _a during B was −17.17°C. Using a T _a of −25.0°C increases the total heat flux from 90 W m⁻² to 136 W m⁻². The maximum supercooling reaches 0.07°C instead of 0.05°C. The surface slush increases from 207 g L⁻¹ at the surface at 23 hours to 298 g L⁻¹. The surface C _si increases from 100 mg L⁻¹ to 118 mg L⁻¹. The size distribution of the C _si does not change qualitatively.

With a T _a of −10.0°C, a mean heat flux of 46 W m² leads to a supercooling of 0.03°C, and 112 g L⁻¹ of ice in the surface. The C _si in the surface decreases to 68 mg L⁻¹ at 23 hours.

Doubling the initial C _s gives an initial concentration of 40 mg L⁻¹ at 0.5 m. The response in C _si is almost linear, giving 204 mg L⁻¹ at 23 hours at the surface. This is also true for lowering the initial C _s to 9 mg L⁻¹ (half the initial volume), which gives a C _si of 52 mg L⁻¹ at 23 hours at the surface. The size distribution of C _si does not change significantly.

Changing the rms fluctuation velocity q changes only the value of K. It does not alter the efficiency of the aggregation process, only the vertical distributions. Increasing q from 7.4 cm s⁻¹ to 14 cm s⁻¹ leads to a 10 times increase in K to 21.2 × 10⁻¹. The initial C _s becomes 38 mg L⁻¹. The increased vertical diffusion reduces the surface C _i to 40 g L⁻¹ at 23 hours, but as the volume does not change, more C _i is in suspension. The surface C _si decreases too, reaching only 43 mg L⁻¹ , but the total sum of C _si increases. If all C _si rose to the surface, it would now be 32.2 g m⁻²,compared to 18.4 g m⁻² in the normal case. Using q = 5.0 cm s⁻¹ (K = 0.34 × 10⁻³) makes the vertical sum of C _si 3.7 g m⁻².

The surface size distribution of C _si changes significantly with q. The high q makes C _si consist of larger grains, while the low q leads to smaller grains in the C _si. This is because high vertical diffusion tends to distribute larger grains more uniformly through the water column.

The model predictions of C _si are close to linear for increased heat fluxes (or T_a with constant wind), and initial concentrations of C _s. The turbulence in the model depends on both q and e, and the associated balance between K and the rise/sinking velocities is quite delicate.

The effect of changing T_a, C_s and q shows that the aggregation process is strongly dependent upon these parameters. They are all needed to estimate the IRS content of newly grown ice in natural settings. In the future, Frasemo will be used to estimate the IRS content of new sea ice on the shallow shelf in the Kara Sea.

6. Conclusion

A new version of a vertical model for frazil-ice and sediment concentrations is developed. New processes included are dynamic growth of frazil ice and sediment-size-based aggregation. Frasemo satisfactorily simulates the tank experiments A–D in Smedsrud (2001). The model is forced with constant (observed) air temperature, wind, turbulence and snow flux. It then reproduces measured time series of water temperature and salinity, and vertical concentrations of sediment and frazil ice.

Frazil crystals are modelled to grow according to their size, and secondary nucleation is set dependent upon the number of crystals at a gridpoint. A counted size distribution can be reproduced using a Nusselt number of 1.5, close to the lower limit, and a low value for the maximum number of ice crystals in a gridpoint.

The vertical diffusion is set constant through experiments A–D based on calculated turbulent rms velocities and dissipation rates from current data. The balance between the rise/sinking velocities and the constant vertical diffusion is somewhat delicate, but values within the observed ranges are used, and found to reproduce measured gradients.

The efficiency of the modelled entrainment process is estimated through a, an aggregation factor.Values for α are in the range 0.0003–0.1, but average values are close to 0.01. This indicates that even with constant forcing, the entrainment process has quite high variations in a tank. These variations cannot be described in detail by the developed model, but results satisfactorily explain the differences in entrainment level between laboratory experiments. Entrainment increases with increasing turbulence and waves, and increases close to linearly with sediment concentration of the water. Lower air temperatures lead to larger heat fluxes from the water to the air, and corresponding higher concentrations of frazil ice, which also increases entrainment.

The model sensitivity to the different parameters is tested, and the average values from experiments A–D should provide a basis for using the developed model in natural case-studies.