Abstract

Cracks in the sea ice cover break the barrier between the ocean and atmosphere, exposing the ocean to the cold atmosphere during the winter. These cracks are known as leads within the continuous sea ice pack and polynyas near land or ice shelves. Sea ice formation starts with frazil ice crystals in supercooled waters, which grow and precipitate to the ocean surface to form grease ice, eventually consolidating into a layer of solid sea ice that grows downward. In this study, a numerical model is formulated to simulate the formation of sea ice in a lead or polynya from frazil ice to a layer of new sea ice. The simulations show the refreezing of a lead within 48 h of its opening. Grease ice covers the lead typically within 3–10 h and consolidates into sea ice within 15–30 h. This study uses its model to simulate an observed polynya event in the Laptev Sea showing the vertical distribution of frazil ice and water supercooling. Sensitivity studies are used to investigate the dependence of ice growth on the ambient environment with the surface wind speed shown to be of greatest importance to lead exposure time and total ice growth. The size and distribution of frazil crystals and the time taken for the lead to freeze over is shown to be highly dependent upon the ambient forcing and lead geometry.

1. Background

Leads and polynyas are ice-free areas within the sea ice cover in which the ocean is in contact with the cold atmosphere in winter. They can form due to warm-water upwelling (sensible heat polynya), katabatic winds or ocean currents that drive newly formed ice away (latent heat polynyas), or when the ice breaks due to internal stresses (leads). Leads are typically long thin features 10 m to 1 km wide and up to 100 km long (Wilchinsky et al. 2015), whereas polynyas are defined as “any nonlinear shaped opening enclosed in (sea) ice” (Maqueda et al. 2004, p. 1). Leads and latent heat polynyas are areas of rapid frazil ice production in winter due to the large heat fluxes from the exposed ocean surface, that is, close to its freezing temperature, to the air, which is normally a lot colder.

When pure water is cooled in laboratory conditions, the homogeneous nucleation of frazil ice crystals requires supercooling of up to 40°C (Mossop 1955). For the case of the ocean surface considered in this study, the supercooling will be quenched by the nucleation of frazil crystals onto impurities in the water, organic or inorganic, or ice crystals from water vapor sublimating as it encounters cold air (Osterkamp et al. 1974) or drops of water, which are ejected into the air and freeze before they drop back into the water (Gosink and Osterkamp 1986). Field observations have found that the greatest supercooling observed before the onset of frazil ice is about 0.1°C (Carstens 1966; Osterkamp et al. 1973; Ushio and Wakatsuchi 1993).

The shape of frazil ice crystals has been observed and studied (Kumai and Itagaki 1953; Arakawa 1955; Williams and Chalmers 1967) with the majority of crystals taking the form of a disc. While they can start out in the shape of six-pointed stars, hexagonal plates, spheres, or dendritic ice, they all evolve into disk shapes in turbulent water (Daly 1984). The size of the crystals has been observed to be in the range of 0.05 mm to several millimeters (Morse and Richard 2009). Once there is an initial seed, the number of frazil crystals may grow quickly through secondary nucleation, also known as collision breeding. This process breaks smaller ice crystals off from larger ones during collisions, either between two larger crystals or between crystals and hard surfaces. The new, smaller crystals then act as nuclei for crystal growth. This creates a positive feedback process, since a higher number of crystals will create more collisions (Clark and Doering 2009).

Previous efforts in modeling the formation of frazil or new ice in leads and latent heat polynyas fall into two categories: process models, such as the one presented in this paper, and components of sea ice climate models. The simplest representations of frazil or new ice formation come in sea ice climate models where a thin layer of ice is assumed to form in the open ocean, which then grows (Maykut 1986; Leppäranta 1993). Pease (1987) presents a simplified model of the wind-driven opening of a latent heat polynya in order to investigate the sensitivity of polynya size to environmental conditions. Although the formation and collection of frazil ice crystals in this study are very simply parameterized, it enables the prediction of the time and length scales associated with the formation of a latent heat polynya.

The dynamics of multiple-sized frazil ice crystals at the surface of the ocean was first considered by Svensson and Omstedt (1994) and then Svensson and Omstedt (1998), using a model, which is an ancestor of the model presented in this paper. The second of these studies gives the vertical distribution of frazil crystals in the upper ocean. This multicrystal size class formulation was used by Smedsrud and Jenkins (2004) for the consideration of underice shelf water plumes and then further developed for the same task by Holland and Feltham (2005, hereinafter HF), who added a depth dependence to the consideration of frazil crystals. The HF model includes a prescription of frazil crystals precipitating onto an ice shelf, a method that we show in this study is also valid for frazil crystals precipitating onto the ocean surface. The dynamics of the ocean under a lead during a refreezing event has been considered by Skyllingstad and Denbo (2001) in a two-dimensional, large-eddy simulation, which includes a frazil ice concentration, although this concentration is of a single-sized frazil crystal. Skyllingstad and Denbo focus on changing ocean turbulence rather than the refreezing of the lead.

In this study, we endeavor to expand the knowledge of frazil ice formation in leads and latent heat polynyas. In particular, we use the work of Svensson and Omstedt (1998) and HF to describe the vertical structure of frazil formation and ocean supercooling immediately following the opening of a lead and within a wind-driven polynya and expand it to describe the subsequent collection and consolidation of new ice at the ocean surface and the salt content of both the ocean and new sea ice cover. We investigate the sensitivity of both the frazil ice and resulting sea ice over to environmental conditions and model parameters. We consider leads and wind-driven polynyas: regions of open water in the sea ice cover that are wide enough (>10 m; Smedsrud and Skogseth 2006) to allow the wind-driven collection of frazil ice crystals ice across it to form a grease ice cover. We do not explicitly consider sensible heat polynyas; though the model presented can be applied to them in future study.

This paper documents a model of ice formation in a lead or polynya of greater complexity and scope than previous models. To represent the changing physical processes occurring during the refreezing of a lead, the model is split into three phases:

  • the supercooling of open waters and initial frazil ice formation based on the model of HF made applicable to a lead or polynya through the novel consideration of a varying salinity profile and sea ice–specific boundary conditions (section 2a);

  • the herding of precipitated frazil ice to form a grease ice layer as described by the model of Smedsrud (2011; section 2c); and

  • the thickening of the resulting sea ice cover treated as a mushy layer (Feltham et al. 2006), allowing us to calculate the congelation of new ice and drainage of brine from the sea ice cover (section 2d) using the model introduced by Wells et al. (2011).

Our paper is structured as follows: In section 2 and the  appendix, we present a description of our model. In section 3a, we use our model to simulate the first 48 h of the refreezing of a lead. In section 3b, we use our model to simulate an observed polynya event in the Laptev Sea, and in section 3c, we present a sensitivity study of the refreezing process to model parameters, ambient conditions, and geometry. Our concluding remarks are presented in section 4.

2. Model description

a. Frazil ice formation

The frazil ice component of the model is based on HF, whose model equations are presented in the  appendix. The model considers the concentration of frazil crystals of multiple size classes within the water column. Melt, growth, and secondary nucleation of the frazil crystals are represented by the movement of crystals from one size class to another. As discussed in the introduction, the spontaneous formation of crystals with the ocean is unlikely, so all the frazil ice within our model comes from an initial seeding of ice crystals. There is no mechanism within the model for the growth or formation of frazil ice crystals in an area of zero total crystal concentration. All crystal growth must come from smaller ice crystals, either from an initial seeding or secondary nucleation. We present modifications to the model of HF: the inclusion of water salinity as a variable, modified secondary nucleation, the use of an updated parameterization for the crystals rise velocity (Morse and Richard 2009), and adjusted boundary conditions to apply the model to the ocean surface.

We consider a mixture of frazil ice and seawater with density ρm defined by

 
formula

where C is the frazil ice volume concentration of the mixture, ρI is the density of ice, and ρO is the density of seawater. The values used for these constants, along with all other parameters, are given in Table 1. We define Nice different size classes of crystals. By discretizing the crystal sizes, we simplify the problem of growth and melting of crystals to transfer between size classes. The total concentration is

 
formula

where Ci is the volume concentration of frazil ice in size class i. The temperature of the ice crystals is assumed to be at the freezing temperature of the fluid part of the mixture. HF found that the results of simulations were qualitatively similar for Nice = 10 and Nice = 200, and so we continue to use Nice = 10. Following numerical experiments, we set our radii range from 0.01 to 2 mm (0.01, 0.02, 0.05, 0.1, 0.2, 0.3, 0.5, 0.8, 1, and 2 mm). As we are assuming that the crystals have a constant thickness to radius ratio (ar = 0.02), we can define the crystal volume , where ti = arri is the thickness of a crystal with radius ri.

Table 1.

List of model parameters and constants.

List of model parameters and constants.
List of model parameters and constants.

We use the ice concentration [(A1)] and temperature [(A2)] equations of HF along with a variable salinity S that alters the freezing temperature of seawater through the equation of state given by Millero (1978) for salinities between 4 and 40 g kg−1 with

 
formula

where η is the depth below mean sea level, a = −0.0573°C (g kg−1)−1, b = 0.0832°C, and c = −7.61 × 10−4 °C m−1. The salinity equation is derived the same way HF derives their temperature equation [(A2)] with

 
formula

where u is the ocean fluid velocity, νT is the turbulent diffusivity, and −g′ is the net discharge of water released due to crystal melting {equal to minus the rate of frazil formation +g′ in the temperature equation [(A2)]}.

We consider a uniform lead or polynya that is laterally invariant across its width and length. The equations of HF are simplified by assuming that vertical processes dominate with ∂/∂x = ∂/∂y = 0. As with HF, we assume the crystal rising-driven advection and turbulent diffusion of the mixture are much larger than the fluid advection giving rise to wCi/∂z = wT/∂z = wS/∂z = 0. With these assumptions, the balance equations of frazil ice concentration [(A1)], water heat [(A2)], and salt [(3)] become

 
formula
 
formula
 
formula

where is the frazil crystal interaction term of HF described in the  appendix, and wi is the buoyancy-driven rise velocity for a crystal of size class i.

The chosen profile of turbulent diffusivity νT considers waves at the surface, wind-driven ocean currents, and water density–driven convection (Haidvogel and Beckmann 1999). We assume the effects of the first two processes are even throughout the mixed layer, giving a constant turbulent diffusivity for a stable density profile with νT1 = 1 × 10−2 m2 s−1 for z < Hmix and νT1 = 1 × 10−4 m2 s−1 for zHmix, where z is the depth from the surface, and Hmix is the depth of the mixed layer (Haidvogel and Beckmann 1999). The vertical gradient in νT is therefore zero except at z = Hmix. As the term is bounded above by in our mixed layer structure profile and it acts at a point only, its effect is negligible to the frazil heat and salt balances. To keep the overall density profile stable, we use the parameterization of turbulent diffusivity given by Haidvogel and Beckmann (1999) with a localized high mixing rate of νT = 1 m2 s−1 wherever the local density profile is unstable with

 
formula

where erf( ) is the error function. The approximation for the mixing and turbulent diffusivity is crude to allow us to focus this study on the complex nonlinear frazil ice dynamics. Improvements to the mixing scheme, including gradients in νT and mixing induced by the movement of frazil crystals, can be included in future study.

The water density profile is derived using the vertical derivative of the equation of state given by HF, calculating the nondimensional density profile from the temperature and salinity profiles with

 
formula

where βT = 3.87 × 10−5 °C−1 and βS = 7.86 × 10−4 (g kg−1)−1, as given by HF. The localized high mixing rates introduce gradients into the νT profile that will be small compared to the local high mixing rate, and we assume will have negligible impact.

1) Secondary nucleation

Secondary nucleation occurs when crystals collide and an ice fragment breaks off, becoming the nucleus of a new crystal. Here, we present formulation for this process derived through the consideration of a moving volume of frazil crystals of a particular size colliding with all other frazil crystals present in the fluid. An assumption made in this formulation, as made by HF is that any new crystals formed through collisions are inserted into the smallest size class. HF gives the secondary nucleation rate from each size class being proportional to its concentration Ci. We propose that the rate of secondary nucleation from a crystal class i is proportional to both the concentration Ci. and the total concentration C [as in (1)] to represent collisions between crystals of different size classes. We define the secondary nucleation rate similar to Svensson and Omstedt (1994) as

 
formula

where Wi represents the velocity of crystal motion within the fluid with , where ϵ is the turbulent dissipation rate, ν0 is the molecular viscosity of seawater, wi is the rising velocity (as described in the following section), anuc is a calibration parameter, and is the effective crystal radius—the radius of a sphere with the same volume of the crystal. Note that Wi represents all turbulent motion of crystals within the fluid, whereas Wi discussed below represents the net vertical buoyancy-driven velocity. The formulation of is for 2 ≤ iNice, with . The value of anuc has been studied by Radia (2014) in a simplified well-mixed model. Following this study, we use a base value of anuc of order 10−6 with variations to the parameter discussed later.

2) Frazil rise velocity

We use the formulation of a frazil crystal’s rise velocity wi (mm s−1) in relationship to its diameter di = 2ri in millimeters, as given by Morse and Richard (2009) with

 
formula

There is a corrigendum to this paper, which results in a formulation giving rise velocities approximately 70% of those given in (6). Considering the corrigendum, we still use the formulation as given by (6) as it has better fit with alternative formulations and observations of rise velocities, as discussed by Morse and Richard (2009), and is within the range of uncertainty for the crystal radii considered in this study (B. Morse and M. Richard 2015, personal communication).

3) Frazil precipitation

For the consideration of frazil crystals precipitating onto the underside of an ice shelf, HF describe a viscous sublayer where viscous stresses dominate and a buffer layer below where both viscous and turbulent stresses are important. The sizes of these layers in dimensionless distance is defined by , where is the friction velocity, z is the coordinate normal to the surface, and ν is the molecular viscosity. The friction velocity is defined by , where Dc = 1.5 × 10−3 is the quadratic drag coefficient and Ua is the wind speed. A viscous boundary layer at the air–ocean interface is present but can be narrower than next to a solid surface, with Wu (1984) finding the layer thickness of z+ = 4 to 8, depending on wind velocity compared to z+ = 7 as used by HF. We use the same methodology as HF to model the frazil precipitation at the ocean surface. We note, however, that the transfer of wind stress to the ocean is not properly accounted for in our model, which uses a wind speed–independent turbulent diffusivity in the mixed layer.

The crystal volume–mass balance in the buffer layer is applied as a surface boundary condition at z = 0 to the frazil equations in the mixed layer (see Fig. 1). As with HF, we integrate vertically across the nondimensional buffer layer with the simplifications

 
formula

as the concentration at the top of the buffer layer is zero and the rate of change in crystal concentration and crystal interaction terms on the rhs now correspond with the depth-integrated values of these quantities, with units of meters per second. The boundary condition for the ith crystal size class at the surface is thus

 
formula

where is the precipitation term per unit length, as given by HF and defined in the  appendix. Note that there is no inconsistency in applying (7) derived from the depth integral of the mass balance as a boundary condition to the frazil equation [(4a)] since the buffer and viscous layers are collapsed to zero in these equations.

Fig. 1.

Schematic of the surface of a lead of width Llead. The buffer layer is for ice formation in open water with all precipitated frazil ice p′ blown to the edge of the lead by a wind of speed Ua to collect as a grease cover with a shape defined by hg(xg) and L, where xg is the distance from the edge of the lead.

Fig. 1.

Schematic of the surface of a lead of width Llead. The buffer layer is for ice formation in open water with all precipitated frazil ice p′ blown to the edge of the lead by a wind of speed Ua to collect as a grease cover with a shape defined by hg(xg) and L, where xg is the distance from the edge of the lead.

b. Surface heat balance

For the surface of the lead, both the open ocean and grease ice, we use the surface energy balance

 
formula

(e.g., Ebert and Curry 1993) where FLW is the incoming and is the outgoing longwave radiation, T0 = T|z=0 is the ocean and grease ice surface temperature (assumed to be equal) and ϵ0 is the emissivity, is the net shortwave radiation that is set to zero for this study (where i0 is the fraction of radiation that is not absorbed near the surface), and α is the albedo; Fsens is the sensible heat emitted with

 
formula

where ρa and ca are the density and specific heat capacity of dry air, CT is an atmospheric turbulence transfer coefficient for both heat and moisture further described by Ebert and Curry (1993), and υa and Ta are the 5-m wind speed and air temperature. The term Flat is the latent heat with

 
formula

where Hυ is the latent heat of vaporization, qa is the 5-m humidity, and q0 is the surface humidity with

 
formula

where pυ is the surface atmospheric vapor pressure for surface temperature T0 at the lead or new ice cover surface, and patm is the atmospheric surface pressure, as given by Ebert and Curry (1993). A sensitivity study into the role of atmospheric pressure was performed by Radia (2014), with no significant changes to the model result when varying from the values given by Ebert and Curry. All parameters used in the surface flux equations are from Ebert and Curry (1993) and listed in Table 1. As we consider the water and ice part of the ocean surface separately and assume that the grease ice layer is at the ocean freezing temperature (as described in the following section 2c), the ocean surface temperature boundary condition is

 
formula

where acover is the fraction of the lead covered by grease ice, and c0 is the heat capacity of pure water.

c. Grease ice

The frazil crystals, which are precipitated into the viscous sublayer, are swept laterally by the wind and collect at one side of the lead or polynya as grease ice of volume Vg. We calculate the volume of precipitated frazil ice in all categories , which entered the viscous sublayer, where we assume it stays and cannot reenter the water column. We use the parameterization of Smedsrud (2011), which has been validated by the experiment of Naumann et al. (2012) to calculate the vertical cross-sectional shape of grease ice in a lead and how it affects the ice cover and hence the heat loss at the surface.

As with Smedsrud (2011), we assume that grease ice is 25% solid ice and 75% water. Taking the precipitated ice from the viscous sublayer , we can calculate the total volume of precipitated ice Vp and the total volume of grease ice Vg per unit width of the lead by

 
formula

The grease ice is assumed to be pushed against the pack (floe) ice edge by the wind, and by using the grease ice layer collection parameterization of Smedsrud (2011), we calculate the width of the grease ice cover L from the volume of grease ice and wind speed Ua with

 
formula

where ρa is air density, Ca = Dc is the open-ocean drag coefficient, and Kr ≈ 100 N m−3 is the grease ice resistance force constant. This formulation assumes that the grease ice is in hydrostatic equilibrium and has a uniform cross-sectional shape along the length of the lead. We define acover = L/Llead as the proportion of the lead covered in grease ice. For numerical stability the point at which the lead closes in the model is set as acover = 0.95.

The formed grease ice consolidates, that is, increases its solid fraction, through the heat loss to the atmosphere, while leaving its total volume unchanged with a solid fraction ϕg = ϕg0Vg/Vp, with ϕg0 = 0.25 as the initial grease ice solid fraction.

d. Congelation ice growth and brine drainage

Once the grease ice cover reaches acover = 0.95, the model enters a new phase with all the remaining frazil crystals precipitated onto the grease ice. The amount of ice left in the water was found to be small compared to the amount of ice that had precipitated up until this point. In the absence of observational evidence or a physically based theory, we decided to keep the bulk salinity of the grease ice layer constant during the consolidation phase for model simplicity. The grease ice solid fraction increases as latent heat is lost to the atmosphere until ϕg = 0.7, at which point it is treated as a developing sea ice cover as discussed below.

The developing sea ice cover is modeled using the mushy layer equations of Feltham et al. (2006) for the vertical temperature profile, and the bottom ablation model of Notz et al. (2003) is used for the ice–ocean interface. The initial ice state is the grease ice distributed evenly over the lead or polynya with thickness Vp/0.25, as the total precipitated volume was created only by precipitating crystals forming grease ice with a solid fraction of ϕg = ϕg0 = 0.25. The consolidating grease ice described by (10) increases the values of ϕg and Vg with Vp unchanged. The bulk salinity is Sbulk = (1 − ϕg)SO, where SO is the ocean surface salinity. An initial constant vertical temperature profile of TI = Tf [(2)] is prescribed here, and the temperature profile of the ice domain then develops according to

 
formula

where ceff and keff are the effective volumetric-specific heat capacity and thermal conductivity of sea ice, and AR represents the absorption of shortwave radiation that is zero for winter simulations. As with Feltham et al. (2006), the values ceff and keff are calculated from

 
formula

where Tf(Sbulk) is the freezing temperature for the bulk salinity Sbulk, Tf(0) is the freezing temperature of pure water, θ = TITf(0), ci is the volumetric heat capacity of pure ice, and kbi, kb, and ki are the conductivities of bubbly ice, brine, and pure ice given by

 
formula

where ka = 0.03 W m C−1 is the conductivity of air, and Va = 0.025 is the fractional volume of air in sea ice (Feltham et al. 2006). It is possible that the parameter values given by Feltham et al. (2006) may be improved upon to better describe the recently consolidated grease ice layer with 30% liquid content, but the investigation and tuning of these parameters is beyond the scope of this study.

At the surface of the ice the boundary condition is a modified version of (9) with . The boundary condition at the bottom of the ice is linked to the ocean temperature through the ice–ocean interface temperature Tint and salinity Sint with TI = Tint = T, where Tint = Tf(Sint), as in (2). The interface conditions are solved along with the rate of growth of the sea ice dh/dt through a modified version of the Stefan condition (extending Notz et al. 2003) with

 
formula

where ϕ = 1 − Sbulk/Sint; αs and αh are turbulent exchange coefficients for salt and heat, as given by Notz et al. (2003) and listed with other parameter values in Table 1; is the friction velocity; Focean is the heat flux from the ocean; and T and S are for the ocean far below the sea ice, which we take at 60-m depth, where we assume the temperature and salinity conditions will remain constant at the values given by the frazil ice model. Equations (13) with Tint = Tf(Sint) comprise a quadratic equation in Sint, Tint, or dh/dt, which is solved exactly. The ice growth salt release is calculated from this quadratic system and gives the salinity boundary conditions for (4c).

We allow brine to drain out of the ice cover, lowering the bulk salinity with dSbulk/dt = −Fbrine, using the formulation of Wells et al. (2011) that describes brine convection through brine channels with

 
formula

where Rm is the porous medium Rayleigh number, and its critical limit is Rc = 10. The porous medium Rayleigh number is

 
formula

where γ = 0.023C0/(CEC0); β is the haline expansion coefficient; ΔC = CES, where CE = T/a is a salinity calculated from the temperature at the vertical midpoint of the sea ice; μ is the dynamic viscosity and ka is the thermal diffusivity of the liquid fraction; and Π0 is the permeability of sea ice as given by Golden et al. (2007), as Π0 = 3 × 10−8(1 − ϕ)3.

3. Model simulations

The models presented in the preceding sections were coded into a FORTRAN computer program [full model code available in supplemental material, with an earlier version in Radia (2014)], using Numerical Algorithms Group (NAG) routines D03 PCA and D03 PZF to solve the system of water equations [(4)] with boundary conditions of (7), (9), and ∂S/∂z = 0 for the ocean surface and ∂/∂z = 0 for the deep ocean at 100-m depth. Once the lead has frozen over, (12) is solved using the NAG routines to give the temperature profile within the thickening ice cover with (13) giving the thickness of the ice cover and the interface temperature. The spatial discretization was 1000 points at 0.1-m spacing for the ocean and 200 points at variable uniform spacing for the congealing sea ice. The model time step is 1 s over 2 days of model time.

a. Reference run

Here, we consider the formation of ice and evolution of the ocean mixed layer in response to the sudden formation of a lead in the ice cover. We have chosen a mixed layer depth of Hmix = 30 m and salinity of 33.0 g kg−1 to represent average ocean conditions for the winter Eurasian and Makarov Basins (Peralta-Ferriz and Woodgate 2015). We simulate the adjustment from initial conditions characteristic of the ocean beneath an intact ice cover. To initiate the model, we set a small amount of frazil ice, C0 = 4 × 10−6, to be present at t = 0, evenly distributed over all crystal size classes in the top 30 m of the model domain (see Fig. 2). This is in accordance with the well-mixed box investigated by Tsang and Hanley (1985) and modeled by Radia (2014). The initial profiles of temperature and salinity were chosen to best replicate a realistic upper layer of ocean with a 30-m mixed layer (Fig. 2). We initialize the model with a water temperature close to the freezing temperature (TTf) in the mixed layer and higher in the ambient water below. We set the initial salinity to be lower in the mixed layer (S ≈ 33.0 g kg−1) and higher in the ambient water (S ≈ 34.0 g kg−1).

Fig. 2.

Initial profiles of (left) frazil ice concentration, (center) water ambient and freezing temperatures, and (right) salinity. The plotted frazil ice concentration represents a single size class, with all size classes having equal concentration and a total frazil ice concentration of 4 × 10−6.

Fig. 2.

Initial profiles of (left) frazil ice concentration, (center) water ambient and freezing temperatures, and (right) salinity. The plotted frazil ice concentration represents a single size class, with all size classes having equal concentration and a total frazil ice concentration of 4 × 10−6.

Here, we present the results of a reference run with the following parameters: FLW = 168 W m−2, FSW = 0 W m−2, Ua = 5 m s−1, and Tair = 249°K or −25°C, as given by Ebert and Curry (1993) for the Arctic winter. We consider ice growth within a lead of width Llead = 150 m. The initial conditions are shown in Fig. 2, and the results are shown in Figs. 3 to 8. The first stage of the model has the water column adjusting from the initial conditions to the imposed forcing. The seeding of frazil crystals partially melts and rises to the near surface where the ocean begins to supercool. There is a delay before a “blooming” of frazil crystals. During this “bloom” the initial seeding of crystals is now contained within the top 2 m of the mixed layer, and their concentration increases from ≈10−5 to ≈10−3 over less than 0.9 h and start to precipitate to the ocean surface. The model is now in a quasi-steady state, and we define the beginning of this as the time where the total precipitation rate first exceeds 10−6 m2 s−1. The crystals continue to precipitate to the surface until time , where the lead is covered by grease ice. The grease ice then consolidates until time , where its solid fraction reaches 70%. From this point to the end of the run at 48 h, the ice cover continues to consolidate and thicken.

Fig. 3.

Vertical profiles of total frazil ice crystal concentration, water temperature, and salinity during the first = 5.5 h of the reference run. The left-hand plot is for the total frazil ice concentration with the lines at selected intervals chosen to show the changes to the mixed layer during this period. The center temperature plot shows the ocean temperature at selected intervals (solid lines) and the freezing temperature at the dashed line.

Fig. 3.

Vertical profiles of total frazil ice crystal concentration, water temperature, and salinity during the first = 5.5 h of the reference run. The left-hand plot is for the total frazil ice concentration with the lines at selected intervals chosen to show the changes to the mixed layer during this period. The center temperature plot shows the ocean temperature at selected intervals (solid lines) and the freezing temperature at the dashed line.

Fig. 4.

The distribution of frazil ice crystals in the open water in a newly formed lead. (a),(c) The mean concentration over the mixed layer of 10-m depth of each individual size class with colors shown at the bottom of the figure. (b) The concentration weighted average crystal radius of the crystals in the mixed layer (solid line) and precipitated to the ocean surface (dashed line). (d) The precipitation rates of each individual size class. The vertical dashed line in each plot represents the time as discussed in section 3a.

Fig. 4.

The distribution of frazil ice crystals in the open water in a newly formed lead. (a),(c) The mean concentration over the mixed layer of 10-m depth of each individual size class with colors shown at the bottom of the figure. (b) The concentration weighted average crystal radius of the crystals in the mixed layer (solid line) and precipitated to the ocean surface (dashed line). (d) The precipitation rates of each individual size class. The vertical dashed line in each plot represents the time as discussed in section 3a.

Fig. 5.

Solid fraction and thickness of the new ice cover. Before time the solid fraction (solid line) is as calculated from (10) with the early high fraction due to the total volume of precipitated grease ice being very low. From to the end of the run, the solid fraction increases due to brine drainage. The mean thickness (dashed line) before is given as Vp/0.25, from onwards it is calculated from the consolidating sea ice cover. The maximum grease ice thickness (dotted line) is calculated from the equations of Smedsrud (2011), similarly to (11).

Fig. 5.

Solid fraction and thickness of the new ice cover. Before time the solid fraction (solid line) is as calculated from (10) with the early high fraction due to the total volume of precipitated grease ice being very low. From to the end of the run, the solid fraction increases due to brine drainage. The mean thickness (dashed line) before is given as Vp/0.25, from onwards it is calculated from the consolidating sea ice cover. The maximum grease ice thickness (dotted line) is calculated from the equations of Smedsrud (2011), similarly to (11).

Fig. 6.

Total volume of precipitated grease ice (solid) and grease ice cover (dashed) for the reference run.

Fig. 6.

Total volume of precipitated grease ice (solid) and grease ice cover (dashed) for the reference run.

Fig. 7.

Internal temperature of the mushy layer sea ice from to the end of the run. The vertical axis shows the negative distance from the sea ice surface. Note that the later contours extend further as the ice has thickened by this point. The dashed line indicates the initial temperature profile set at t = .

Fig. 7.

Internal temperature of the mushy layer sea ice from to the end of the run. The vertical axis shows the negative distance from the sea ice surface. Note that the later contours extend further as the ice has thickened by this point. The dashed line indicates the initial temperature profile set at t = .

Fig. 8.

Salt release from the thickening congelation ice from to the end of the run. The top figure is the salt release rate and the bottom figure is total salt released from brine pockets (solid) and due to ice growth (dashed).

Fig. 8.

Salt release from the thickening congelation ice from to the end of the run. The top figure is the salt release rate and the bottom figure is total salt released from brine pockets (solid) and due to ice growth (dashed).

For the first phase of the model with t < , the open water in the lead is rapidly cooled from above and becomes supercooled, as seen in the profiles for 0.3 h in Fig. 3. During this period, the crystals in the initial ice state rise to the surface due to their increased buoyancy, reducing the average crystal size in the mixed layer (Fig. 4b). There is a low (≈10−8 m s−1) rate of surface precipitation that temporarily consolidates, as seen in Fig. 5. From 0.5 h the concentration of 0.15–0.4-mm size crystals starts to increase, initially in the top 1 m of the water column and then deeper as the run continues. The ice growth quenches some of the surface supercooling, although the surface remains supercooled by 0.001°C for the rest of the run.

By the model time of = 0.9 h, the surface ice concentration has increased to allow the model to enter a quasi-steady state. In this state the average crystal size of the frazil ice precipitated on the surface and within the mixed layer remains constant (Fig. 4b). From 1–2 h of model run, the total concentration of frazil ice within the mixed layer continues to increase, both in magnitude and vertical extent, with crystals slowly mixing downward with significant concentrations reaching 15-m depth by = 5.5 h in Fig. 3. The mixed layer is at its freezing temperature for the same vertical extent as the significant concentration of frazil crystals. At = 5.5 h, the top 4 m is supercooled by 0.001°C, and the salinity profile is near vertically uniform. The concentration of each individual size class, although being different in overall magnitude, has the same vertical profile as that shown in Fig. 3. The overall rate of ice precipitation is well correlated to the concentration of frazil ice in the mixed layer, although the representation of the different crystal size classes changes. The 0.3–0.4-mm-sized crystals are the most abundant in both the mixed and precipitated layers. The fraction of the largest crystal sizes (1–2 mm) is low within the mixed layer but makes a significant contribution to the precipitated grease ice (Fig. 4d). The precipitated frazil crystals collect to form the grease ice layer that initially consolidates rapidly due to its low thickness (Fig. 5), although the solid fraction returns to near the 25% initial value after the bloom of frazil crystals. Both the extent of the grease ice cover and the total volume of precipitated crystals increase near linearly (Fig. 6). By = 5.5 h, the ice lead is covered by grease ice with an average thickness of 0.19 m with a solid fraction of 30%.

In the second model phase from to , the grease ice cover consolidates into solid ice. The grease ice fraction increases, removing freshwater from the water fraction that thus becomes more saline. No salt is released to the ocean during this phase, which is a model simplification that could be revisited subsequent to a detailed investigation of the processes involved in grease ice consolidation. During the final model phase of congelation ice growth, the state of the ice cover is determined from the surface boundary conditions and the interface conditions given by equations (13). The temperature profile of the ice given by (12) becomes linear within 2 h of , as shown in Fig. 7 with a surface temperature that approaches −19°C and interface temperature at the ocean freezing temperature of −1.87°C. The heat flux given by (8) for this surface temperature has component values of approximately −70 W for the net longwave radiation, −70 W for the sensible heat flux, and −10 W for the latent heat flux, with a total heat flux of −150 W. This is in comparison to the heat flux of −500 W for the open lead at the beginning of the run. The ice cover thickens at a rate of ≈5.6 × 10−7 m s−1 to 0.227 m by the end of the run, which results in salt release into the ocean from the advancing ice–ocean interface [as given by equations (13)] that is near constant at 1.2 × 10−5 g (kg s−1)−1 (see top of Fig. 8). The bulk salinity of the sea ice decreases from 22.8 to 5.7 g kg−1 due to the brine drainage rate [as given by (14)] that peaks at 20.4 × 10−5 g (kg s−1)−1 and decays to 0.2 × 10−5 g (kg s−1)−1 by the end of the model run, by which time the ice growth salt release has released 0.75 kg and the brine drainage has released 2.8 kg of salt to the ocean per unit area of new sea ice (see bottom of Fig. 8). The salt release from brine pockets and ice growth increases the mixed layer salinity by, on average, 0.5 g kg−1 by the end of the model run, which is an order greater than the 0.04 g kg−1 increase during the initial formation of frazil ice.

b. Recreating the polynya event observed by Dmitrenko et al. (2010)

Dmitrenko et al. (2010) present a study of the formation of the Laptev Sea polynya in the Arctic. Using a range of measurements including radar from Envisat and satellite imagery from TerraSAR-X, they captured a polynya event from start to finish and made measurements of the water salinity and temperature profiles from a mooring attached to the fast ice edge. Here, we recreate the polynya event on 28–30 April 2008 discussed in Dmitrenko et al. using an adapted version of our model. This is achieved by assuming that all the frazil ice precipitated to the ocean surface is blown away by the wind with acover = 0 for the whole run. We set the air temperature at −12°C and wind speed at 8 m s−1, giving a heat flux of 550 W m−2 as observed by Dmitrenko et al. and the air humidity is taken from ERA reanalysis data. The mixed layer depth is set to 6 m, with high turbulent diffusivity of νT = 1 × 10−2 m2 s−1 there to account for high mixing. Below the mixed layer the turbulent diffusivity is lower (νT = 1 × 10−5 m2 s−1), as used in the reference run in section 3a. The initial profiles of temperature and salinity are taken from observations on 24 April (see Fig. 9b), which we assume are similar to the conditions on 26 April when the polynya event started.

Fig. 9.

Temperature and salinity conditions associated with the polynya event from 28 to 30 Apr 2008 documented by Dmitrenko et al. (2010). (a) The observed time series of temperature (red) and salinity (green) in the mixed layer with the polynya event highlighted in blue and (b) the observed profile of temperature (dashed) and salinity (solid) on 24 Apr both from Dmitrenko et al. (2010). (c) Modeled mixed layer temperature (red) and salinity (green) at a depth of 3 m representing the polynya event from 28 to 30 Apr.

Fig. 9.

Temperature and salinity conditions associated with the polynya event from 28 to 30 Apr 2008 documented by Dmitrenko et al. (2010). (a) The observed time series of temperature (red) and salinity (green) in the mixed layer with the polynya event highlighted in blue and (b) the observed profile of temperature (dashed) and salinity (solid) on 24 Apr both from Dmitrenko et al. (2010). (c) Modeled mixed layer temperature (red) and salinity (green) at a depth of 3 m representing the polynya event from 28 to 30 Apr.

Figure 9a shows the time series of temperature and salinity measured at 4.5-m depth during 23 days in April and May, with the 3 days corresponding to the polynya event highlighted. We assume that the mixed layer is well mixed and these values are true for the entire mixed layer, as with the profile in Fig. 9b from 2 days before the event. The observed change in water temperature and salinity from −1.42° to −1.52°C and from 26 to 28 g kg−1 is closely matched by the model results in Fig. 9c, although our model was run with constant forcing causing the changes in temperature and salinity to be linear. The constant forcing results in the model achieving a quasi-steady state after an hour, with a profile of ice concentration shown in Fig. 10, which remains unchanged during the remainder of the model run. The profile of the temperature and freezing temperature also has a similar shape and amount of supercooling (of 0.007°C) over the model run, with average mixed layer temperature decreasing linearly.

Fig. 10.

Frazil ice concentration and temperature profiles within the mixed layer 1 day into the Dmitrenko et al. comparison model run. (left) The crystal classes C1–4 with the remaining classes C5–10 of negligible amount. (right) The water temperature (solid) and freezing temperature (dashed) profiles with a supercooling of 0.007°C at the ocean surface.

Fig. 10.

Frazil ice concentration and temperature profiles within the mixed layer 1 day into the Dmitrenko et al. comparison model run. (left) The crystal classes C1–4 with the remaining classes C5–10 of negligible amount. (right) The water temperature (solid) and freezing temperature (dashed) profiles with a supercooling of 0.007°C at the ocean surface.

The quasi-steady state has a two-layer form. The lower layer from −2- to −4-m depth is well mixed with near-constant ice concentration and with the ocean at its freezing temperature. The ice concentration is split over the three smallest crystal classes with a small contribution from C4. Each size has a peak concentration at different depths showing crystals growing in size as they rise to the surface. This distribution within the crystal classes presents a steady state that does not develop in the reference run in section 3a. The upper layer has a lower concentration of crystals because they are being precipitated to the ocean surface (out of the domain) and then swept away by the wind. The removal of crystals from the upper layer reduces further ice growth, which can only occur on existing crystals, and reduces the quenching of supercooling. In addition, the upper layer has large crystals (average radius of 0.12 mm compared with 0.06 mm in the lower layer) due to their greater buoyancy and thus higher rise speed, and these larger crystals grow more slowly, further reducing their ability to quench the supercooling.

The amount of supercooling (0.007°C) is within the range Dmitrenko et al. derive from their observations during the polynya event, which are variable and peak at 0.02°C. Our results suggest that during the polynya event discussed supercooling is confined to the top 2 m of the ocean with the majority of the frazil crystals found below.

c. Sensitivity studies

We perform sensitivity studies to investigate the role of forcing and model geometry in ice formation and to show the numerical behavior of the model. The forcing is investigated by varying the wind speed Ua and air temperature Ta; the model geometry is investigated by varying lead width Llead, mixed layer depth Hmix, and turbulent diffusivity νT1. The sensitivities of the frazil crystal formation are further investigated by changing the initial concentration of crystals C0 and the secondary nucleation parameter anuc. Important results from the sensitivity runs are displayed in Table 2. For all the sensitivity runs, the model forcing and parameter values are the same as in the reference run unless stated.

Table 2.

Sensitivity studies for varying model parameters: is the time taken for the precipitation rate to exceed 10−6 m2 s−1, is the time taken for the lead to be covered with grease ice, and is the time taken for the grease ice cover to consolidate with a solid fraction of 70%. Max Sc is the maximum ocean surface supercooling TfT for the whole run. The term is the average total frazil ice concentration over the mixed layer, is the average suspended frazil crystal radius, depth Ci is the vertical distance below the ocean surface where the total frazil ice concentration exceeds 10−5, Pr. is the average total surface precipitation rate, and Pr. is the average radius of a precipitated ice crystal all from the period to . The term h is the thickness of the sea ice covering the lead at 48 h of model run time.

Sensitivity studies for varying model parameters:  is the time taken for the precipitation rate to exceed 10−6 m2 s−1,  is the time taken for the lead to be covered with grease ice, and  is the time taken for the grease ice cover to consolidate with a solid fraction of 70%. Max Sc is the maximum ocean surface supercooling Tf − T for the whole run. The term  is the average total frazil ice concentration over the mixed layer,  is the average suspended frazil crystal radius, depth Ci is the vertical distance below the ocean surface where the total frazil ice concentration exceeds 10−5, Pr.  is the average total surface precipitation rate, and Pr.  is the average radius of a precipitated ice crystal all from the period  to . The term h is the thickness of the sea ice covering the lead at 48 h of model run time.
Sensitivity studies for varying model parameters:  is the time taken for the precipitation rate to exceed 10−6 m2 s−1,  is the time taken for the lead to be covered with grease ice, and  is the time taken for the grease ice cover to consolidate with a solid fraction of 70%. Max Sc is the maximum ocean surface supercooling Tf − T for the whole run. The term  is the average total frazil ice concentration over the mixed layer,  is the average suspended frazil crystal radius, depth Ci is the vertical distance below the ocean surface where the total frazil ice concentration exceeds 10−5, Pr.  is the average total surface precipitation rate, and Pr.  is the average radius of a precipitated ice crystal all from the period  to . The term h is the thickness of the sea ice covering the lead at 48 h of model run time.

The model’s greatest sensitivity is to wind speed Ua. A high wind speed causes the grease ice to bunch up thicker, and a higher volume is needed to cover the lead. This causes the runs with increasing wind speed to take longer to form a grease ice cover (increased ), with thicker ice at 48 h compared to the reference run. Increasing the width of the lead Llead has the same result, with a greater average thickness of grease ice required to cover a wider lead. While a thinner grease ice cover consolidates faster with an earlier , subsequent growth of congelation ice is slow compared to the rate of grease ice accumulation. This results in the higher wind speeds and wider lead having later and and a thicker sea ice cover at the end of the simulation. Slower wind speed (Ua = 2.5 m s−1) or a thinner lead (Llead = 50 m) has the opposite effect. Our models’ greatest sensitivity to wind speed is in contrast to that of Pease (1987), who found that their model was more sensitive to the surface heat flux. This is almost certainly because the Pease (1987) model used a constant grease ice collection depth, unaffected by the wind speed and in contradiction to observations (Naumann et al. 2012). It is difficult to directly compare our results to those of Pease (1987), as we have a prescribed lead or polynya width, though this can be considered in a future study.

Increasing wind speed also causes a greater surface heat loss, as does a lower air temperature Ta. Both these cases (Ua = 10 m s−1 and Ta = −35°C) cause an increase in maximum surface supercooling and average frazil ice concentration. The increased average frazil concentration is due to both an increased crystal growth rate and the frazil crystals extending deeper into the mixed layer. For all the model runs, the precipitated crystals are on average larger than the suspended crystals. The development of the mixed layer with a bloom of frazil crystals happens earlier for increased surface heat flux. Once the grease ice cover forms, an increased surface heat flux causes it to consolidate faster (earlier ), and a thicker sea ice layer is formed due to having more time for the congelation ice to grow and due to the ice growth rate being higher. The ice growth salt release increases, and the brine pockets salt release is unchanged. For the runs with a thicker grease ice layer (Ua = 10 m s−1 and Llead = 250 m), the brine pocket salt release increases due to greater collected grease ice volume. The cases of lower surface heat flux (e.g., Ta = −15°C) have the opposite results: decreased supercooling, frazil crystal concentration, and depth covered; later times , , and ; and a thinner ice cover at the end of the run.

For certain model setups, the quasi-steady state described in the reference run does not occur: for example, a thin mixed layer of Hmix = 10 m and a low surface heat flux with Ta = −15°C, shown in Table 2 as combined (1). In this particular model setup, 90% of the initial seeding of frazil crystals precipitate to the surface before the bloom of frazil crystals can occur. There remains a 1-m layer of low concentration (≈10−6) of ice crystals in the largest size class. The layer slowly increases in concentration, and by 5.3 h the precipitation rate has become significant and continues to rise until the lead freezes over.

For the majority of sensitivity runs with Hmix = 30 m, the lead freezes over before frazil crystals mix down to fill the layer (see frazil depth column of Table 2). One exception to this is for Ua = 10 m s−1, where there is both an increased surface heat flux and greater time for the frazil crystals to mix down through the mixed layer. The other exception is for an increased mixed layer turbulent diffusivity νT1 → 2νT1, where the frazil crystals are mixed throughout the mixed layer from 3 h onward. The increased mixing delays the start of the frazil crystal bloom, as the near-surface crystal concentration takes longer to increase to a level where secondary nucleation can begin. Reducing the turbulent diffusivity νT1 → 0.7νT1 also delays the frazil bloom as more crystals precipitate out of the ocean in the first hour of the model run. After the bloom the crystals cover a shallower portion of the mixed layer. Reducing νT1 by more than 0.7 causes the model to become unstable. For the case of a 20-m-deep mixed layer, despite the distance being greater than the depth of crystals in the reference run, the crystals eventually cover the entire mixed layer that is at its freezing temperature. There is not enough model time or great enough surface heat flux for this to happen in the reference run. For all the sensitivity studies with equal surface heat flux and grease ice collection depth, those with a deeper coverage of frazil crystals (increased ν, Hmix = 20 m) have a slight increase in ice thickness at the end of the model run due to the addition of more suspended crystals to the grease ice cover at . For the given reference run, increasing the mixed layer depth to Hmix = 50 m makes no difference to the model results, as the frazil ice crystals similarly do not fill the mixed layer. Also decreasing the mixed layer salinity to 28 g kg−1, representing a fresh Beaufort Gyre, for example, does not change the model results apart from a small decrease in the brine channel drainage rate in the consolidated new ice cover. These runs are not included in Table 2.

Increasing the secondary nucleation parameter anuc changes the distribution of frazil crystals in the mixed layer with both suspended and precipitated crystals, having on average a smaller diameter. For the case of anuc = 10−3, the smaller crystals have a lower buoyancy and so mix down through the mixed layer quicker and initially precipitate more slowly to surface with a later . However, the lower buoyancy causes an increase in average crystal concentration in the mixed layer that results in an increased overall precipitation and earlier . For the case of anuc = 10−1, the decrease in crystal size has reached the limit of the smallest crystal sizes with 90% of the suspended and precipitated crystal sizes being of the first two size classes. The crystals that grow from the first and smallest size class to the second size class quickly cause secondary nucleation and return to the smallest size class. Despite the large changes in the distribution of crystals, the resulting ice cover at the end of the runs with altered anuc differs little from the reference run, which is a consequence of the same atmospheric forcing.

The model also has a small sensitivity to the choice of limits for acover and ϕg that control the timing of model phases. Taking a smaller value of acover = 0.8, for example, causes an earlier and thinner initial covering of grease ice. The grease ice takes less time to consolidate (earlier ), but the final thickness of ice cover is thinner. As the salt content of the consolidating ice cover is conserved, changing model phase at a different value of ϕg only changes the timing of when the sea ice starts consolidating. Changing model phase at ϕg = 0.5, for example, causes only a slight increase in thickness (0.01 m) and total salt release (0.05 kg) by the end of the model run due to the increased time spent consolidating compared to the reference run. The majority of salt release comes within the first 10 h of brine pocket release, unaffected by the choice in limit for ϕg.

4. Concluding remarks

Leads and polynyas are areas of rapid ice formation and negative buoyancy production (through salt release) to the ocean in winter. Accurately calculating the time taken to freeze over a lead has been shown to have significant impact on the mass balance of the entire sea ice pack (Wilchinsky et al. 2015). A one-dimensional model of ice formation in a lead or polynya has been created that accounts for supercooling of the mixed layer, frazil ice production, precipitation of frazil ice to grease ice, the herding of grease ice, and the consolidation of grease ice followed by congelation ice growth. Our model is the only model to address all these processes, and it combines several submodels of different growth stages and preserves the heat, salt, and mass budgets.

Our reference simulation gives insight into the complex behavior of growth, precipitation, and nucleation of frazil crystals within the mixed layer under a newly formed lead in unprecedented detail. Our model successfully recreates the temperature and salinity conditions in the Laptev Sea polynya as observed by Dmitrenko et al. (2010). We show the vertical structure of supercooling and frazil crystal formation; the top 2 m of the mixed layer are supercooled by the atmosphere with frazil crystals growing and quickly precipitating to the surface.

We have used our model to show how the size distribution of frazil crystals is highly sensitive to external conditions, the size of the lead, and atmospheric conditions above it. Our sensitivity studies show that the greatest controlling factor for the refreezing of a lead is the surface wind speed Ua. High wind speeds result in greater heat loss to the atmosphere and cause the precipitated frazil crystals to collect into a thicker grease ice layer so that the lead takes longer to freeze over and ultimately results in a thicker layer of new sea ice. It is important to understand whether the phenomena of all frazil crystals being precipitated out of the mixed layer beneath a lead is a physical possibility or an artifact of our model setup.

Extension of our model to a two-dimensional cross section of a lead, or introduction of a new parameterization, would allow the influx of seeding crystals or particulates from the sides. The inclusion of a horizontal dimension would also allow a greater sophistication in the way frazil crystals precipitate and collect into a grease ice layer. Such a model would also be able to consider the collection and consolidation of the grease ice layer continuously without the need of distinct model stages. While our model has generated quantitative insight into the process of lead refreezing, it is important to validate the results of this study through observations of the mixed layer and ice cover in a newly formed lead.

Acknowledgments

We would like to thank other CPOM members and the anonymous reviewers for helping to improve this study. This work was funded under NERC GRANT NE/K011510/1.

APPENDIX

Model Equations of HF 

Here, we present the model equations used in this study that are described in HF) for the water column beneath an ice shelf that is represented as a mixture of water and ice crystals. For the full mathematical expansion and derivation of these equations, we refer the reader to HF and Radia (2014). Applying the Boussinesq approximation to the conservation of mass of the mixture of ice and water, HF derive volume balances for the frazil ice and water fractions of the mixture. These two equations are combined to give the volume balance equation for the ith class of frazil crystals with

 
formula

where u is the velocity of the seawater and ice mixture, νT = μt/ρI is the turbulent diffusivity with μt the turbulent eddy viscosity, and is the interaction term between the different crystal size classes. The temperature balance for the fluid part of the mixture is derived by balancing the rate of change in temperature, advection of heat by both the fluid part and crystal rising, turbulent mixing, and a heat flux QT dependent on crystal freezing and melting. The term QT is expanded to represent both latent and sensible heat through the net rate of change in crystal volume g′, and the equation

 
formula

is derived with Hf the latent heat of ice fusion, and c0 is the specific heat capacity of pure water.

a. Frazil interaction term

We use the frazil interaction term of HF to describe in (A1) and (4). The crystals are modeled as being circular disks with a constant ratio of radius to thickness, so we can define their size with one parameter only, radius r. As the concentration is split into Nice size classes, we can model crystal growth , melting , and secondary nucleation as transfer between different size classes.

The interaction term of size class i is defined by the difference between the sources, melting from size class i + 1 and growth from size class i − 1, and the sinks, melting and growth from size class i. Since the different size class crystals have different volumes, we must account for this when transferring crystals, and so if a crystal in size class (i) of volume υi grows and enters class (i + 1), then crystal growth in term Ci must be at least to transfer one crystal to size class (i + 1). The interaction term is defined as

 
formula

where H = He(TfT), where He is the Heaviside step function to account for the cases of T < Tf (supercooled) and T > Tf (not supercooled), and and are as defined by HF with boundary conditions of . The net discharge of water g′ is calculated from summing through all the interaction terms with .

b. Frazil ice precipitation

The rate of frazil precipitated onto the surface of the ocean as grease ice is parameterized as with HF using the relation

 
formula

where and are the precipitation from turbulent and laminar flows, Ri is the Richardson number to distinguish between the flows with the critical limit RiC = 0.25, and d = 8 is a transition smoothing parameter. Using the parameterization of McCave and Swift (1976), the turbulent precipitation is given as

 
formula

where He is the Heaviside step function, U is the depth-mean velocity in the mixed layer, and UT is the root-mean-square tidal velocity. The term UCi is the critical deposition velocity for each crystal size class, which can be calculated from the Shields criterion and the quadratic law:

 
formula

The laminar precipitation is given as

 
formula

where the step function is used to ensure that .

To determine whether the flow is turbulent or laminar, HF use a modified form of the Richardson number to provide a single dimensionless quantity that represents the effects of frazil–seawater mixture viscosity, shear, and stability on turbulence with

 
formula

where

 
formula

is calculated using the predicted velocity profile in the vicinity of z+ = 35 by the law of the wall with depth-mean velocity U and Von Karman’s constant κ = 0.41. The molecular viscosity is given as a function of the ice concentration with

 
formula

where ν0 is the pure water viscosity when C = 0.

REFERENCES

REFERENCES
Arakawa
,
K.
,
1955
:
Studies of the freezing of water III
.
J. Fac. Sci., Hokkaido Univ., Ser. 2
,
4
,
355
358
. [Available online at http://hdl.handle.net/2115/34220.]
Carstens
,
T.
,
1966
:
Experiments with supercooling and ice formation in flowing water
.
Geophys. Norv.
,
26
,
1
18
.
Clark
,
S. P.
, and
J. C.
Doering
,
2009
:
Frazil flocculation and secondary nucleation in a counter-rotating flume
.
Cold Reg. Sci. Technol.
,
55
,
221
229
, doi:.
Daly
,
S. F.
,
1984
: Frazil Ice Dynamics. CRREL Monogr., No. 84-1, U.S. Army Corps of Engineers Cold Regions Research and Engineering Laboratory, 56 pp.
Dmitrenko
,
I. A.
, and Coauthors
,
2010
:
Observations of supercooling and frazil ice formation in the Laptev Sea coastal polynya
.
J. Geophys. Res.
,
115
,
C05015
, doi:.
Ebert
,
E. E.
, and
J. A.
Curry
,
1993
:
An intermediate one-dimensional thermodynamic sea ice model for investigating ice-atmosphere interactions
.
J. Geophys. Res.
,
98
,
10 085
10 109
, doi:.
Feltham
,
D. L.
,
N.
Untersteiner
,
J. S.
Wettlaufer
, and
M. G.
Worster
,
2006
:
Sea ice is a mushy layer
.
Geophys. Res. Lett.
,
33
,
L14501
, doi:.
Golden
,
K. M.
,
H.
Eicken
,
A. L.
Heaton
,
J.
Miner
,
D. J.
Pringle
, and
J.
Zhu
,
2007
:
Thermal evolution of permeability and microstructure in sea ice
.
Geophys. Res. Lett.
,
34
,
L16501
, doi:.
Gosink
,
J. P.
, and
T. E.
Osterkamp
,
1986
: Frazil ice nucleation by ejecta from supercooled water. Proc. Int. Association of Hydraulic Research, Symp. on Ice, Iowa City, Iowa, Institute of Hydraulic Research, The University of Iowa, 249–264.
Haidvogel
,
D. B.
, and
A.
Beckmann
,
1999
: Numerical Ocean Circulation Modelling. Series on Environmental Science and Management, Vol. 2, World Scientific, 344 pp.
Holland
,
P. R.
, and
D. L.
Feltham
,
2005
:
Frazil dynamics and precipitation in a water column with depth-dependent supercooling
.
J. Fluid Mech.
,
530
,
101
124
, doi:.
Kumai
,
M.
, and
K.
Itagaki
,
1953
:
Cinematographic study of ice crystal formation in water
.
J. Fac. Sci., Hokkaido Univ., Ser. 2
,
4
,
235
246
.
Leppäranta
,
M.
,
1993
:
A review of analytical models of sea-ice growth
.
Atmos.–Ocean
,
31
,
123
138
, doi:.
Maqueda
,
M.
,
A. J.
Willmott
, and
N. R. T.
Biggs
,
2004
:
Polynya dynamics: A review of observations and modeling
.
Rev. Geophys.
,
42
,
RG1004
, doi:.
Maykut
,
G. A.
,
1986
: The surface heat and mass balance. The Geophysics of Sea Ice, N. Untersteiner, Ed., Springer, 395–463.
McCave
,
I. N.
, and
S. A.
Swift
,
1976
:
A physical model for the rate of deposition of fine-grained sediments in the deep sea
.
Geol. Soc. Amer. Bull.
,
87
,
541
546
, doi:.
Millero
,
F. J.
,
1978
: Freezing point of seawater. Eighth Rep. of the Joint Panel on Oceanographic Tables and Standards, UNESCO Tech. Papers in Marine Science 28, Annex 6, 29–31.
Morse
,
B.
, and
M.
Richard
,
2009
:
A field study of suspended frazil ice particles
.
Cold Reg. Sci. Technol.
,
55
,
86
102
, doi:; Corrigendum,
61
,
143
145
, doi:.
Mossop
,
S. C.
,
1955
:
The freezing of supercooled water
.
Proc. Phys. Soc.
,
68B
,
193
, doi:.
Naumann
,
A. K.
,
D.
Notz
,
L.
Håvik
, and
A.
Sirevaag
,
2012
:
Laboratory study of initial sea-ice growth: Properties of grease ice and nilas
.
Cryosphere
,
6
,
729
741
, doi:.
Notz
,
D.
,
M. G.
McPhee
,
M. G.
Worster
,
G. A.
Maykut
,
K. H.
Schlünzen
, and
H.
Eicken
,
2003
:
Impact of underwater-ice evolution on Arctic summer sea ice
.
J. Geophys. Res.
,
108
,
3223
, doi:.
Osterkamp
,
T. E.
,
C. S.
Benson
,
R. E.
Gilfilian
, and
T.
Ohtake
,
1973
: Winter history of an Alaska stream. Annual Rep. of the Geophysical Institute, University of Alaska, Fairbanks.
Osterkamp
,
T. E.
,
T.
Ohtake
, and
D. C.
Warniment
,
1974
:
Detection of airborne ice crystals near a supercooled stream
.
J. Atmos. Sci.
,
31
,
1464
1465
, doi:.
Pease
,
C. H.
,
1987
:
The size of wind-driven coastal polynyas
.
J. Geophys. Res.
,
92
,
7049
7059
, doi:.
Peralta-Ferriz
,
C.
, and
R. A.
Woodgate
,
2015
:
Seasonal and interannual variability of pan-Arctic surface mixed layer properties from 1979 to 2012 from hydrographic data, and the dominance of stratification for multiyear mixed layer depth shoaling
.
Prog. Oceanogr.
,
134
,
19
53
, doi:.
Radia
,
N. V.
,
2014
: Frazil ice formation in the polar oceans. Ph.D. dissertation, University College London, 167 pp.
Skyllingstad
,
E. D.
, and
D. W.
Denbo
,
2001
:
Turbulence beneath sea ice and leads: A coupled sea ice/large-eddy simulation study
.
J. Geophys. Res.
,
106
,
2477
2497
, doi:.
Smedsrud
,
L. H.
,
2011
:
Grease-ice thickness parameterization
.
Ann. Glaciol.
,
52
,
77
82
, doi:.
Smedsrud
,
L. H.
, and
A.
Jenkins
,
2004
:
Frazil ice formation in an ice shelf water plume
.
J. Geophys. Res.
,
109
,
C03025
, doi:.
Smedsrud
,
L. H.
, and
R.
Skogseth
,
2006
:
Field measurements of Arctic grease ice properties and processes
.
Cold Reg. Sci. Technol.
,
44
,
171
183
, doi:.
Svensson
,
U.
, and
A.
Omstedt
,
1994
:
Simulation of supercooling and size distribution in frazil ice dynamics
.
Cold Reg. Sci. Technol.
,
22
,
221
233
, doi:.
Svensson
,
U.
, and
A.
Omstedt
,
1998
:
Numerical simulations of frazil ice dynamics in the upper layers of the ocean
.
Cold Reg. Sci. Technol.
,
28
,
29
44
, doi:.
Tsang
,
G.
, and
T. O.
Hanley
,
1985
:
Frazil formation in water of different salinities and supercoolings
.
J. Glaciol.
,
31
,
74
85
, doi:.
Ushio
,
S.
, and
M.
Wakatsuchi
,
1993
:
A laboratory study on supercooling and frazil ice production processes in winter coastal polynyas
.
J. Geophys. Res.
,
98
,
20 321
20 328
, doi:.
Wells
,
A. J.
,
J. S.
Wettlaufer
, and
S. A.
Orszag
,
2011
:
Brine fluxes from growing sea ice
.
Geophys. Res. Lett.
,
38
,
L04501
, doi:.
Wilchinsky
,
A. V.
,
H. D. B. S.
Heorton
,
D. L.
Feltham
, and
P. R.
Holland
,
2015
:
Study of the impact of ice formation in leads upon the sea ice pack mass balance using a new frazil and grease ice parameterization
.
J. Phys. Oceanogr.
,
45
,
2025
2047
, doi:.
Williams
,
R. B.
, and
B.
Chalmers
,
1967
:
Morphology of ice solidified in undercooled water
.
J. Phys. Chem. Solids
, Crystal Growth (Suppl.),
739
743
.
Wu
,
J.
,
1984
:
Viscous sublayer below a wind-disturbed water surface
.
J. Phys. Oceanogr.
,
14
,
138
144
, doi:.

Footnotes

Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JPO-D-16-0224.s1.

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

Supplemental Material