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 Reference SmedsrudSmedsrud (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 Reference Gosink and OsterkampGosink and Osterkamp (1983)

The vertical diffusion of T _w, S _w, C _i and C _s is modelled with a fully implicit diffusion scheme (Reference PatankarPatankar, 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 (Reference GillGill, 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⁻³ (Reference SmedsrudSmedsrud, 2001) is reduced to the standard “open-sea” value of c _H = 1.1 × 10⁻³ (Reference Simonsen and HauganSimonsen 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 Reference SmedsrudSmedsrud (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 (Reference NeilsenNeilsen, 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 (Reference Gosink and OsterkampGosink and Osterkamp, 1983; Reference Hammar and ShenHammar and Shen, 1995; Reference Svensson and OmstedtSvensson 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 (Reference Gosink and OsterkampGosink 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 (Reference Svensson and OmstedtSvensson and Omstedt, 1994), or higher when N_u is set dependent on the turbulent dissipation rate (ε) and the Kolmogorov length scale (λ) (Reference Hammar and ShenHammar and Shen, 1995). The thermal conductivity of sea water is set to K _w = 0.564 W m⁻¹ °C⁻¹ (Reference CaldwellCaldwell, 1974).

Equation (5) is comparable with Equation (3) in Reference Svensson and OmstedtSvensson 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 Reference Hammar and ShenHammar and Shen (1995).

The models of Reference Svensson and OmstedtSvensson and Omstedt (1994) and Reference Hammar and ShenHammar and Shen (1995) were validated by fresh-water laboratory experiments like that of Reference CarstensCarstens (1966). Reference Svensson and OmstedtSvensson 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 Reference CarstensCarstens (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 Reference Tsang and HanleyTsang 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. Reference OmstedtOmstedt (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 (Reference Hammar and ShenHammar 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 (Reference DalyDaly, 1984). The simplified approach of Reference Svensson and OmstedtSvensson 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 Reference Svensson and OmstedtSvensson 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 Reference Hanley and TsangHanley 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 (Reference ShamlouShamlou, 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 (Reference Osterkamp, Gosink, Barnes, Scholl and ReimnitzOsterkamp 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 Reference SmoluchowskiSmoluchowski (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λ (Reference ShamlouShamlou, 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 (Reference Tennekes and LumleyTennekes 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 (Reference Osterkamp, Gosink, Barnes, Scholl and ReimnitzOsterkamp and Gosink, 1984), and also associated with the “frazil stickiness”“discussed by Reference Kempema, Reimnitz, Clayton and PayneKempema and others (1993).

The theoretical framework has been tested on mono-sized polystyrene and spherical latex particles, and fits well to experimental data (Reference Higashitani, Yamanchi, Matsuno and HosokawaHigashitani and others, 1983; Reference Gierczycki and ShamlouGierczycki 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 Reference SmedsrudSmedsrud (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 Reference SmedsrudSmedsrud (1998) α ∼0.025, but it might vary between the different size classes.

The results from the box model (Reference SmedsrudSmedsrud, 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 Reference SmedsrudSmedsrud (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 (Reference MartinMartin, 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⁻¹ (Reference NürnbergNürnberg and others, 1994).

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 Reference SmedsrudSmedsrud, 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 (Reference CarstensCarstens, 1966; Reference DalyDaly, 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 Reference Gosink and OsterkampGosink 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 Reference SmedsrudSmedsrud, 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 (Reference SmedsrudSmedsrud, 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 (Reference Forest, Lunardini, Wang, Ayorinde and SodhiForest, 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 (Reference SmedsrudSmedsrud, 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.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 (Reference Mellor and YamadaMellor and Yamada, 1982), using q and ∊ (Reference SmedsrudSmedsrud, 2001). Measurements of q are robust, but the 95% confidence level of the ∊ estimate covers an order of magnitude (Reference SmedsrudSmedsrud, 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 Reference Hammar and ShenHammar 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 Reference Svensson and OmstedtSvensson 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 (Reference Gosink and OsterkampGosink and Osterkamp, 1983). Alternative equations exist, which generally give higher rise velocities. The velocity range 0.01–80 mm s⁻¹ is suggested by Reference DalyDaly (1984) for the size range given in Table 2. By setting , the thickest crystals in the range of Reference Gosink and OsterkampGosink and Osterkamp (1983), w _i(5) = 73.9 mms⁻¹, close to the suggested value from Reference DalyDaly (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. Reference SmedsrudSmedsrud (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 Reference SmedsrudSmedsrud (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 (Reference SmedsrudSmedsrud, 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” (Reference SmedsrudSmedsrud, 1998) seems appropriate.