Simulation of supercooling and size distribution in frazil ice dynamics

cold regions science and technology ELSEVIER

Cold Regions Science and Technology 22 (1994) 221-233

Simulation of supercooling and size distribution in frazil ice dynamics Urban Svensson a, Anders Omstedt b aComputer-aided Fluid Engineering, Krokvdgen 5, S-60210 Norrk6ping, Sweden bSwedish Meteorological and Hydrological Institute, S-60176 NorrkiJping, Sweden (Received 1 March 1993; accepted after revision 14 September 1993)

Abstract

The objective of the work presented is to formulate a mathematical description of frazil ice dynamics. The formulation is to be in balance with the current knowledge of the physical processes, for example secondary nucleation. As the knowledge of some of these processes is fragmentary, this means that a conceptually simple formulation is sought. A number of processes are known to influence the supercooling rate and the frazil ice production. The present formulation accounts for the following processes: initial seeding, secondary nucleation, gravitational removal, growth due to cooling of water volume and flocculation/break up. Equations are formulated for these processes considering a resolution in time and radius of particles but not in space (well-mixed jar). The equations are solved using a simple explicit numerical scheme. Preliminary results indicate that the model can be calibrated to describe the experimental results reported in the literature. It is mainly the supercooling curves that are used for comparison but some information about the crystal size distribution is also considered. It is to be noted that the model is calibrated to fit the experiments, due to the lack of detailed mathematical description of some of the physical processes. Sensitivity analysis is also used in order to establish that the model behaves according to experimental findings and expectations. The main conclusion of the study is that a fairly simple mathematical model can be formulated and calibrated, which fits the experimental data reported in the literature hitherto. It is further concluded that a resolution in radial space gives additional insight into the dynamics of the process. The evolution of the size distribution and its sensitivity to seeding and dissipation rate has been predicted with results that look physically plausible.

1. Introduction

The initial ice formation in turbulent water bodies is due to frazil ice, which is the term used for the small ice crystals formed in supercooled water. Frazil ice formation is common in rivers, large lakes and in the oceans. The world's largest frazil ice producer is in fact the Antarctic Ocean, where the main sea ice formation is due to frazil ice production. In rivers, frazil ice can cause serious problems at intakes to, for example, hydro-

power plants. Similar problems have also been encountered at fresh water intakes in large lakes. The initial ice formation in turbulent waters begins when the water is supercooled and seeded with ice crystals from the atmosphere, see Fig. 1. As the ice formation proceeds, the suspended ice crystals increase in size and number, stick to the bottom or rise towards the surface, flocculate and form clusters with high porosity. Due to the action of waves the frazil slush freezes together and forms pancake ice. In the ocean, salt is rejected

0165-232X/94/$07.00 © 1994 Elsevier Science B.V. All fights reserved SSDI 0165-232X (93)E0030-M

222

U. Svensson, A. O m s t e d t I Cold Regions Science a n d Technology 22 (1994) 2 2 1 - 2 3 3

'b

Q

Net heat toss •

"o

¢~] Net heat toss

Net heat toss

° • ~ °oe o

o

o

o

z~

~ •° o

Turbutence

Turbutence

c~b ~ o 8:' ° o

Sott-

~ rejection LJ

o

Turbutence

7 i l l A \ \ \ \ ~ i i l i i ~ \ \ \ \ ' i l f l A \

Supercooting, seeding

Frazit ice formation, growth, break-up, multiplication, snltreJection

Floccul.otion, surfoce ice and anchor ice formation

~e

d

~

Net heat loss

~

Supercooting Max

Ice drift

s,upercool.ing I ~ecover

Time

o o

Turbutence

Icecovergrowth,freeze-up

Cooing

iceformation

Time temperature evotution

Fig. 1. Sketch of initial ice formation in turbulent water bodies, together with a typical time history of water cooling during ice formation.

from the ice crystals and forms dense water which may sink into deeper layers. During ice formation, the crystal density and size distribution vary considerably (Daly, 1991 ). During seeding, the typical dimension of the crystals is 10 -5 to 10 -4 m. Once frazil ice starts to form, the typical size is ranging from I 0-5 to 10-2 m. At the larger limit, the crystals become morphologically unstable, and ice fragments are shed from the crystal, generating new and smaller crystals. The typical dimension of frazil flocs is 10 - 3 t o 10 - 1 m, and due to the action of waves they become circular. A general review on frazil ice formation can be found in Daly (1984). There is a large number of models on mixed layer dynamics available, only few have, however, considered frazil ice dynamics• Omstedt and Svensson (1984) developed a frazil ice

model for a turbulent Ekman layer. The basic idea was that the frazil ice can be treated by a boundary layer approach, in which buoyancy effects due to ice crystals are included. The model was verified against laboratory data by Omstedt (1985a) and extended to the upper layers of the ocean by Omstedt (1985b). The model by Omstedt and Svensson (1984) has also been applied to rivers by Nyberg (1986). In these models only one mean size of the frazil ice crystals were considered. The reason was that multiplication of ice crystals was regarded as the dominant reason for ice growth. From comparison between modelled and measured data the supercooling rates and frazil ice growth were realistically simulated applying a given mean frazil ice crystal shape and a constant Nusselt number for the heat transfer calculations

U. Svensson, A. Omstedt I Cold Regions Science and Technology 22 (1994) 221-233

(Omstedt, 1985a ). As ice growth is partly due to crystal growth and partly due to multiplication, an assumption that only considers multiplication disregards the crystal growth from small seed crystals to flocs of frazil ice. The growth of the very first small crystals influences, however, the temperature recover from maximum supercooling. In the modelling efforts, the frazil ice dynamics has generally been treated in a quite rough way. Basically only the temperature response due to frazil ice formation has been considered. The time/temperature evolution has been the main concern. A notable exception is the work by Mercier (1984); he developed a model which considers most aspects of frazil ice dynamics mentioned above and also reported good agreement with data from laboratory experiments. Also the model by H a m m a r and Shen (1991 ) considered the size distribution of frazil ice by solving a number of volumetric concentration equations. These equations were solved together with the hydrodynamical equations, as applied to a river channel. The purpose of the present paper is to present a new approach to the modelling of supercooling and crystal size distribution. The model is fairly simple but accounts for the basic physical processes of the frazil ice dynamics. The outline of the paper is as follows: in the section to follow the model is presented together with a discussion of the physical processes that are under consideration. The next section gives the results based upon comparison between model simulations and measured data. Finally, a discussion about frazil ice dynamics and some conclusions are given in Section 4.

2. Conceptual/mathematical

223

models with vertical resolution, e.g. the model by Omstedt and Svensson (1984).

2.1. Physical processes considered The frazil ice particles are classified into Ndiscrete radius intervals, within which all particles are assumed to be of equal radius. The number of particles in each group, ni, is assumed to be a function of, see Fig. 2: - Initial seeding. It will be assumed that a number of particles, ninit(ri), is present in each radius interval at the time when supercooling starts. -Crystal growth. Due to the supercooling the crystals will grow. The heat transfer rate is a function of the turbulence characteristics and the dimensions of the crystal, often expressed as a Nusselt number. - S e c o n d a r y nucleation. This is the process by which smaller crystals are produced from collisions between larger ones. It will be assumed that this process is a source for the smallest radius interval and a sink for the rest. - Flocculation/break up. It is assumed that flocs can form by, for example, the sintering mechanism. Large flocs can also be expected to break up under the action of, for example, shear. It will be assumed that the net effect of the two processes is formation of larger aggregates. - Gravitational removal. The particles are given a rise velocity and hence have the capability to rise to the surface. If a steady state situation is to be considered, it is necessary to consider the gravitational removal, as the integral balance is between the surface heat loss and ice production. For a steady state to prevail it is further

model

In this section the basic assumptions and ideas forming the simulation model will be presented and put in mathematical terms. One of the more basic assumptions is that the model has zero space dimensionality, i.e. a well-mixed box is assumed. This is for convenience, and it is anticipated that the model will later be modified to fit

rf

ri ~

rN

J,occu[./ Breakup Crystal growth ~Sec. Nucl.

Fig. 2. Schematic figure of physical processes considered.

224

U. Svensson, A. Omstedt / Cold Regions Science and Technology 22 (1994) 221-233

necessary to assume that the gravitationally removed ice does not accumulate at the surface, as this would reduce the area for surface cooling. 2.2. The number continuity equation

The above processes can, formally, be represented in a crystal number continuity equation of the following kind: dn,

dt (I
N ~" E oljnj -olin ` j=2 (i=1) (2
-- flini

~Pi- lni- 1

-t-

(1 < i < N - 1) (2
-

lint

+

Fi_lni_l

-

yin, (1
(1)

(I
where a j, pit, 7~ and F; are coefficients [ s- 1] giving the strength of the process considered. The factor ~ is the ratio between volumes of particles of two neighbouring radius intervals. Eq. (1) gives the evolution of the particle size distribution from an initial stage, given the values of the coefficients. The left hand side of Eq. ( 1 ) gives the change in the number of crystals in radius interval i: This change is due to the processes on the right hand side of Eq. (1). The first term, secondary nucleation, gives a source for the smallest radius interval, i = 1, and a corresponding sink for other intervals. The two terms in flocculation/break up and crystal growth are due to the discretization concept used; crystals are entering from a smaller radius and leave the present interval to enter the next higher interval. The gravitational removal is a sink term for all radius intervals.

For discs it has been found that the heat transfer, and hence the growth rate, is larger at the edges than at the faces of the disc. In this study it will be assumed that the disc-shaped crystals grow only at the edges. The heat transfer per unit area can then be estimated as: Nuk

q=-~--(rice-rwater)

[ W i n -2 ]

(2)

where N u is a Nusselt number, k thermal conductivity, d a characteristic length (particle thickness), Ticeice temperature and Twaterwater temperature. The characteristic length will be chosen as the thickness of the disc; the thickness does not change as it is assumed that the disc grows only at the edges. As d<< r/, where r/is the dissipative length scale of the turbulence, it can be assumed that the Nusselt number is of the order of one. It will be assumed that Nu = 1.0. The ice volume produced per unit time by crystal growth for the radius interval i can now be written as: OVi

k =NUll( Tic~ - TwatCr)Ai¢¢ni/ (Pic¢L ) [m 3 s -1 ]

(3)

where Aide is the active freezing area per crystal, Pic~ density of ice and L the latent heat of ice. The active freezing area is 2rtr,d. By noting the difference in particle volumes between the crystals in the present radius interval and those in the next larger one, it is straight forward to calculate how many crystals are to be moved to the higher interval. The coefficient F~, see Eq. (1), is thus given by:

ov,

/', =-~-/Avol,

[s-ll

(4)

where Avo~,is the difference in crystal volumes mentioned above.

2.3. Crystal growth

If frazil ice - - introduced by, for example, seeding - - is present in a turbulent water volume, crystals will start to grow when the temperature drops below the temperature of freezing.

2.4. Secondary nucleation

Secondary nucleation is the term used for the production of small crystals by removal of nuclei

U. Svensson, A. Omstedt / Cold Regions Science and Technology 22 (1994) 221-233 IO0

225

I Fully Turbulent Range

A 10l -

oE

-

Intermediate Range

16 2 Stokes

,o,

i i~r t03

10"z

j

I Ill

I

t tl[

t

J ~t[

I0" I0O Veloeity (era/s)

J

L

I0'

1

i0z

Fig. 3. Terminal rise velocityof frazil disks. Basic figure from Daly (1984). lihood for collisions in the expression. A crystal in relative movement to the fluid will sweep a volume, di, during a time interval dt:

Number/(unif volume and radius)

1013

°°o

oo

o Oo

Ai = Urnr~ dt

o

o

(5)

where Ur is the relative velocity. A collision frequency for the n~ crystals in radius interval i can then be formulated as:

1011

f~o.~nAin,/dt

[s -1]

(6)

I0Ic

109

where n is the average number of crystals per unit volume. The relative velocity is assumed to be due to turbulent fluctuations and gravitational rise:

~o

108

t : r ~ , / t : , + Urise 2 lO-S

105

10-4

10-~

lOZ t-(m)

Fig. 4. Predicted size distribution spectra. Straight line gives the - 4 slope.

from the surface of parent crystals. This mechanism is poorly understood (Daly, 1991 ), and it is presently not clear if "detachment of surface irregularities by fluid shear" or "breeding by collisions" is the dominating process. Considering this uncertainty, the simplest possible formulae will be sought in the present analysis. It does, however, seem to be relevant to include the like-

(7)

where Ut=(1/15)w2(~/v)W2d (E is turbulent dissipation rate, v kinematic viscosity and d crystal diameter) is a turbulent velocity scale according to Delichatsios and Probstein (1975) and Ur~se the gravitational rise velocity (to be discussed below). Eq. (7) is correct when one of the two velocity scales approaches zero; there are, however, other possibilities to add the two scales, which may be more adequate when they are of the same magnitude. Eq. (6) is thus the expression to be used for the secondary nucleation. The proportionality - -

226

U. Svensson, A. Omstedt / CoM Regions Science and Technology 22 (1994) 221-233

b)

Q)

TIME (rain) 0.00

TIME ( r a i n )

0.00

2

0

4

6

0

2

4

6

8

-0.02

°0

I-

-0.02

fi.

I-

~

\.//

-0.04

m Q.

tu

-0.04

I-

x . d e,

-0.06

-0.06 C)

TIME ( r a i n ) 0

0.00

<

-0.02

,

2 ,

,

4 ,

,

6 ,

,

8 ,

tO ,

•

a. w I--

",*.,..%,,/

-oo4

Fig. 5. Comparisons with the experimental results (marked by II) from Michel (1963), top left figure, and Carstens (1966), top right and bottom figures. Solid lines are predictions according to Mercier (1984) and dashed lines according to the present model. In (a) the following values were used: Q= 1200 W/m 3, e=50× 10.4 m2 s -3, ninit=0.8 X 104. In (b): Q= 1400 W / m 3, e= 13X 10 -4 m 2 S-3, ninit~--1.1X 104. In (c): Q= 550 W/m 3, ~=4X 10 -4 m 2 s -3, ninit= 1.5 X 104. and hence the calibration m factor will be introduced by setting an u p p e r limit on n, thus: aj = nA ini/dt

[S - ] ]

(8)

lence and hence a reduction o f the collision frequency. It thus seems to be logical to limit n by the calibration factor nmax, as this will put a restriction on the collision frequency. 2.5. F l o c c u l a t i o n / b r e a k up

n=max~

lnj,nmax)

[m-3]

(9)

As already m e n t i o n e d , the process o f secondary nucleation is presently not understood, and it is hence not surprising that a calibration par a m e t e r is required. T h e r e is, however, a physical argument why it is i n t r o d u c e d the way it is. If the c o n c e n t r a t i o n o f frazil ice is increasing, one m a y expect gravitational d a m p i n g o f the turbu-

U n f o r t u n a t e l y also the processes o f flocculation and break up o f frazil ice crystals are poorly u n d e r s t o o d (Daly, 1991). By the sintering m e c h a n i s m it is, however, believed that crystals can form flocs. Sintering results from the tendency o f crystals to minimize their surface energy. In the present investigation it will be assumed that a transport to larger scales is the net effect o f flocculation and break up. It will further be

8

U. Svensson, A. Omstedt I ColdRegions Science and Technology22 (1994)221-233

227

that the drag coefficient will be influenced by the turbulence (Daly, 1984). During one time step, dt, gravity is assumed to remove all particles located in the volume A Un. ,edt, where A is the water surface area. The removal per unit time and volume is thus:

time (s) 16.00 1200 1000

yi=Urise/D

800

[S - 1 ]

(12)

where D is the depth of the well-mixed box. Ice removed by gravity can not enter the volume again and is assumed not to reduce the free surface.

600

t~0200I

tO0

I

~o ~

I

I

I

I

I

~00 5oo 600 700

~

:~

Q tW/

m3

2. 7. Overall heat balance )

Fig. 6. Relation between surface heat loss and time to maximum supercooling in the experiment (solid line) by Hanley and Michel ( 1977 ). Points give predictions with the present model with ninit= 12 000.

The water temperature in the well-mixed box can be calculated from the overall heat balance. N--1

d ( p C p T w a t e r ) = - Q n + ~ Qi

(13)

j=l

assumed that the process is more effective for larger crystals. A linear relation is assumed:

fli=anocri/rl

[s -1 ]

(10)

where anoc is a calibration parameter. A value on fli of, for example, 10 -2 means that one percent of the crystals in interval i will move to interval i + 1, per second.

2.6. Gravitational removal Daly (1984) reviews the aspects on the rise velocity of frazil ice discs. His conclusions are summarized in Fig. 3, where also a simplified equation is shown (solid thick line). This equation reads: Un,e = 30r] 2

[m/s]

(11)

where ri is the radius of the disc. In the present study this equation will be used, as it is in fair agreement with Daly's result and also in agreement with Nybrant (1945), Ur= 0.04 for r i = 2 . 5 X 1 0 -3, and Gosink and Osterkamp (1983), Ur=0.02 for r i = l . 2 5 X 1 0 -3. Eq. (11) may need to be modified when we have a mean flow and developed turbulence, as it is known

where p is density for water, C, specific heat for water, Qn net heat loss per unit volume at the surface and Qi release of heat due to freezing for radius interval i. If a steady state is reached, a balance between the surface heat loss and the ice production is required. In the frazil ice regime, a balance is then required between the production and the gravitational removal. Secondary nucleation and flocculation/break up will only redistribute the ice volume, i.e. modify the size distribution.

2.8. Numerical solution The set of equations was solved by a simple explicit (Euler) numerical scheme. Twenty radius intervals were found to give an adequate resolution in radial space and a time step of 0.5 seconds was found to produce a time step-independent solution for most of the cases to be presented.

3. Results

The mathematical model introduced has two calibration parameters, nmax and anoc. These will

228

U. Svensson, A. Omstedt I Cold Regions Science and Technology 22 (1994) 221-233

25O I

a) 00-

500

Time (s)

!

-01 - Number (m3)

c) I

109

Temp.

("C)

p

10 B

b)

Diam. (m)

2x~-

107

~o6 105

__j/

10a7

, I

/f"N._.j /

103

/

10z

>

0

2~0

500 Time (s)

250

dO Time(s)

Fig. 7. Experimental (dashed lines) and calculated (solid lines) time variation of supercooling (a), mean diameter (b) and total number of crystals (c). Experimental results by Daly (pers. commun., 1992).

first be estimated, then some predictions based on the model will be presented.

3.1. Calibration The calibration parameter ano¢ can be estimated by an analysis of the size distribution spectra. According to Mercier (1984, p. 42 ) the particle size distribution can often be described by an equation of the form:

g(l) =AI -b

(14)

where g(l) is the number density distribution (number of particles per unit fluid volume per unit particle length), l the particle length scale and A and b constants. Mercier (1984) finds that b = 4 gives the correct slope of the distribution. In Fig. 4 a calculation (Qn= 1000 W / m 3, = 10- 3 m 2 s - 3, ninit = 10 4) is shown together with the straight line giving the - 4 slope. The parameter anoc was in this calculation given a value of 10 -4 s - t , which, as can be seen, gives the correct size distribution. It should be men-

229

U. Svensson, A. Omstedt / C o l d Regions Science and Technology 22 (1994) 221-233

!

0

,o

•

¢

e. I ,I

Time = 150 s

1" "0

j

Time = 300 s

Fig. 8, A pictorial view of number and size distribution at four times.

tioned that the size distribution is sensitive to the value on anon, and it is not possible to get a correct slope without the inclusion of the flocculation/break up process, i.e. with ano~equal to zero. Next we consider the parameter nmax. Some classical laboratory experiments were considered for this purpose, namely Michel ( 1963 ) and Carstens (1966). These have also been analyzed by Mercier, and his specification of the dissipation rate in the experiments will be used. For the

experiments considered there is actually another unknown quantity, namely the initial seeding. In order to avoid the specification of both the amount, duration and size distribution, which all are unknown, we chose to specify a uniform number of crystals, ni,it, in each radius interval at time equal to zero. One can not expect a constant value on ninit for different experiments, and the calibration task is thus to find a common value on nm~x, allowing a variation on/'/init. Re-

U. Svensson, A. Omstedt / Cold Regions Science and Technology 22 (1994) 221-233

230

,,log (N i } loc (N i) s

4

2

" I0 -6

10.5

I0~

I0 -3' "

10~(rrl)

: i0-o

i

i ttlt,,]

i

lO-b

Time = lOO s

f ~IH.I

'

10 ./*

I I.,,,,I

I

r *,.,,.I

t

,

10::r fin)

10 I:J

Time = 200 s

,log{Ni}

tog (N i]

7

6

$

4

J

10-6

10-5

10-~

10-3

lO-2tlm)

T o n e "=- 4 0 0 s

i

104

. ,,.,.I

r

105

I Irlllll

i

10-~

i ''HII~

103

IIIIHI

tO:Pr (rn)

Time = 600 s

Fig. 9. A s e n s i t i v i t y a n a l y s i s ; ( ) r e f e r e n c e c o n d i t i o n s , ( - [] - [] - ) s e e d i n g rate i n c r e a s e d 5 t i m e s , ( - O - © - ) d i s s i p a t i o n rate i n c r e a s e d 5 t i m e s . Size d i s t r i b u t i o n s at 4 t i m e s .

suits of such a calibration are shown in Fig. 5. A value of 4 X 106 on nm~x and the values of ninit

3.2. P r e d i c t i o n s

given in the figure caption were found to give a close agreement with the experimental data.

The values on afloc and nrnax will now be considered as fixed, and the model will be used for predictions.

U. Svensson, A. Omstedt / Cold Regions Science and Technology 22 (1994) 221-233

231

Volum fluxes (rn31s)

5xlC~

-5×1ff -

,

10~

,

,,:,,,I

,

10-5

,

,,,,,,I

,

10-4

,

,,,,,,f

,

,

,J

10-3

Fig. 10. P r e d i c t e d steady state b a l a n c e o f v o l u m e fluxes at different r a d i u s intervals. g r a v i t a t i o n a l r e m o v a l , ( - X - X - ) crystal growth.

The first case considered is the experiment by Hanley and Michel (1977). Their measurements of the relation between time to maximum supercooling and surface heat loss are shown in Fig. 6 together with the predicted curve. The agreement found gives support to the model even if the value on ninit was used as a calibration parameter. Next we consider the time development of the mean radius and the number of crystals. Daly (pers. commun., 1992) is currently undertaking experiments that will demonstrate how the mean diameter and total number of crystals develop in time. Preliminary results from this experiment are shown in Fig. 7. A first attempt to simulate this experiment is also shown in Fig. 7. To somewhat conform with the seeding procedure in the experiment, the initial seeding was made a function of crystal size, i.e. ninit=nlrN/ri with n l = 100. The surface heat flux was put to 1222 W / m 3 and e = 2 0 × 10 -4 m 2 / s 3. As the experimental data are preliminary, only qualitative comparisons ought to be made. Nevertheless, the maximum value in mean diameter and the rapid

,,,f

10-~

,

,

,,,,,,I

_

10~ rim) ) F l o c c u l a t i o n / b r e a k up, ( - O - O - )

increase in number of crystals after the maximum supercooling are common features. A pictorial view of the simulation is given in Fig. 8. The number and size distributions are given as snapshots at four selected times. At 600 seconds conditions are close to steady state. In order to get a clear illustration the number of particles had to be normalized. The normalizing factor was chosen as 104 when the total number of particles was of the order of 105 and 106 when the number of particles was of the order of 109. A linear variation was used in between. Using these normalizing factors, a clear illustration of the size distribution at different times could be obtained, as seen in Fig. 8. It should perhaps also be mentioned that the crystals were distributed randomly in space, as the model has zero space dimensionality. A sensitivity analysis of the simulation was also carried out; the seeding rate was increased by a factor of 5 in the first test, and in the second test the dissipation rate was increased by a factor of 5. The effect on the time development of the size distribution is shown in Fig. 9. Starting with the

232

U. Svensson, A. Omstedt / Cold Regions Science and Technology 22 (1994) 221-233

increased seeding rate one can note that the number of particles will remain large up to the point when a steady state is approached, i.e. at 600 s. It is also clear that secondary nucleation already at 100 seconds has produced small particles, as the increase in this range is of the order of 10. For larger particles the increase is typically about 5 times. The increase in dissipation rate is expected to increase the collision frequency and hence the production of small particles by secondary nucleation. In Fig. 9 it is seen that this process has started already after 200 seconds; it increases the number of particles by two orders of magnitude. The number density for large particles is, as expected, less affected. At 600 s, which is close to steady state, the profile is still affected and can thus also be expected to influence conditions at steady state. The temperature curve is also affected by the seeding and dissipation rate. This will not be discussed here, but one may mention that the maximum supercooling was reduced by about 50% for the two cases discussed. The final case to be discussed deals with the steady state balance of Eq. ( 1 ). A typical situation with Q,= 1000 W / m 3 and ~=20X 10 - 4 m 2 / s 3 was predicted for that purpose. Calculations were continued till steady state conditions prevailed. In Fig. l0 the balance of volume fluxes [m3/s] is shown for different radius intervals. Gravity is always a sink term, as expected, while flocculation/break up redistributes from smaller to larger scales.

4. Discussion and conclusions The mathematical modelling of frazil ice dynamics includes many intricate physical processes, of which some are not well known. There is therefore a need to increase our knowledge of these processes, both through new laboratory and field experiments and through mathematical models which can generalize the observational findings. Up to now the modelling of the ice crystal size distribution during freezing in water bodies has not been extensively dealt with. The present work is one of the first attempts to model ice crystal size distributions. The mathematical

formulation is fairly simple, but several aspects of the dynamics are included. Based upon only two calibration coefficients, one limiting the secondary nucleation and one tuning the flocculation process, the model reproduces some important experimental findings. These are: time/temperature simulations during supercooling and ice formation, which gives reasonable estimates of maximum supercooling and time of recover, realistic recover times for different heat losses, crystal size distribution spectra with - 4 slopes, and realistic simulation of ice crystal growth and size distributions. The modelling approach is most promising, and together with more experimental data on the crystal size distribution it can easily be further developed. It is suggested that the present model should be used in combination with laboratory experiments. By using the numerical model as a synthesis of experimental findings and theoretical considerations, it is expected that improved understanding of frazil ice dynamics can be achieved.

Acknowledgements Special thanks are given to Steven Daly for making unpublished laboratory data available. This work has been financed by the Swedish Meteorological and Hydrological Institute and the Swedish Natural Science Research Council under contract G-GU 9151-304. Key subroutines of the computer code were developed and coded by Hans-Olof Kuylenstierna. The drawing by Ann-Margreth Hoist and the printing by Vera Kuylenstierna are also gratefully acknowledged.

References Carstens, T., 1966. Experiments with supercooling and ice formation in flowing water. Geofys. Publ. Norway, 26 (9): 3-18. Daly, S., 1984. Frazil ice dynamics. CRREL Monograph 841. CRREL, Hanover, NH, 46 pp. Daly, S., 1991. Frazil ice. In: K.C. Cheng and N. Seki (Editors), Freezing and Melting Heat Transfer and Engineering. Hemisphere Washington, DC, Ch. 16, pp. 523-544.

U. Svensson, A. Ornstedt / Cold Regions Science and Technology 22 (1994) 221-233 Delichatsios, M.A. and Probstein, R.F., 1975. Coagulation in turbulent flow: Theory and experiment. J. Colloid Interface Sci., 51 (3): 163-174. Gosink, J.P. and Osterkamp, T.E., 1983. Measurements and analyses of velocity profiles and frazil crystal rise velocities duhng periods of frazil ice formation in rivers. Ann. Glaciol., 4: 79-84. Hammar, L. and Shen, H.T., 1991. A mathematical model for frazil ice evolution and transport in channels. In: Proceedings from Workshop on Ice, Ottawa. Hanley, T.O'D. and Michel, B., 1977. Laboratory formation of border ice and frazil slush. Can. J. Civ. Eng., 4 (2): 153160. Mercier, R., 1984. The reactive transport of suspended particles: Mechanics and modeling. Ph.D. Dissertation, Joint Committee on Oceanographic Engineering, Massachusetts Institute of Technology, Cambridge, MA. Michel, B., 1963. Theory of formation and deposit of frazil

233

ice. In: Eastern Snow Conference, Proc. Annual Meeting, Quebec. Nyberg, L., 1986. Ice formation in rivers. In: C.A. Brebbia and S.A. Orszag (Editors), Numerical Simulation of Fluid Flow and Heat/Mass Transfer Processes. Lecture Notes in Engineering, No. 18. Springer-Verlag, Berlin, pp. 108121. Nybrant, G., 1945. Bidrag till teorin f'6r isbildning i/ilvar. (A contribution to the theory of ice formation in rivers. In Swedish. ) Tekniska Skrifter nr 120, Tekniska Tidskriftsf'6rlaget, Stockholm. Omstedt, A., 1985a. On supercooling and ice formation in turbulent sea water. J. Glaciol., 31 ( 109 ): 263-271. Omstedt, A., 1985b. Modelling frazil ice and grease ice in the upper layers of the ocean. Cold Reg. Sci. Technol., 11 : 8798. Omstedt, A. and Svensson, U., 1984. Modeling supercooling and ice formation in a turbulent Ekman layer. J. Geophys. Res., 89(C1 ): 735-744.