Polynyas in a dynamic-thermodynamic sea-ice model

Polynyas in a dynamic-thermodynamic sea-ice model E. Ö. Ólason and I. Harms Centre for Marine and Atmospheric Research, Institute of Oceanography, University of Hamburg, Hamburg, Germany Received: 28 October 2009 – Accepted: 9 November 2009 – Published: 24 November 2009 Correspondence to: E. Ö. Ólason (einar.olason@zmaw.de) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
One of the most interesting features of sea ice, in a climate context, is the fact that it acts as an insulating layer between the atmosphere and the ocean.As the ice grows, the transfer of heat from ocean to atmosphere is greatly reduced.This slows further ice formation down and allows the atmosphere to cool much more than it otherwise would.
This blanket of ice is subject to atmospheric and oceanic momentum forcing, which can create openings in the cover, allowing direct heat exchange between the ocean and the atmosphere.Where this happens, the heat transfer may increase up to hundredfold, causing vigorous ice formation, brine release and heating of the atmosphere near the opening.This paper focuses on large openings, known as polynyas, but smaller openings are referred to as leads.
Correspondence to: E. Ö. Ólason (einar.olason@zmaw.de)Polynyas and leads are an important part of the climate system at high latitudes.Maykut (1982), for instance, estimates that about 50% of the total atmosphere-ocean heat exchange over the Arctic Ocean in winter occurs through polynyas and leads.During summer, these openings admit shortwave radiation into the ocean, warming it up and thus impacting the heat and mass balance of the ice and ocean (Maykut and Perovich, 1987;Maykut and McPhee, 1995).Arctic polynyas also play a large role in halocline and deep water formation and Winsor and Björk (2000) estimate a mean ice production from all Arctic polynyas of 300±30 km 3 yr −1 .The resulting salt flux is about 30% of the estimated flux needed to maintain the halocline.
In terms of general circulation models, polynyas are modelled using dynamic-thermodynamic sea-ice models.This has been done successfully by a number of researchers; e.g.Marsland et al. (2004), Kern et al. (2005) and Smedsrud et al. (2006).Not all researchers use the same criterion to define a polynya in dynamic-thermodynamic models.The most straightforward approach would seem to be to use ice concentration alone, like Kern et al. (2005).However, Smedsrud et al. (2006)  This appears to be due to a fundamental difference in the model results of Kern et al. (2005) and Marsland et al. (2004) on one hand and that of Smedsrud et al. (2006) on the other.In the former studies a polynya can be characterised as an area of low ice concentration surrounded by ice of high concentration.In the latter case the polynya is an area of thin ice (but high concentration) surrounded by thick ice.
These differences serve as an incentive to take a closer look at how polynyas form in dynamic-thermodynamic seaice models.Highly idealised setups are useful when trying to understand the basic processes involved and so we choose to revisit a study originally made by Bjornsson et al. (2001).In their study, Bjornsson et al. (2001) compared the granular Published by Copernicus Publications on behalf of the European Geosciences Union.model of Tremblay and Mysak (1997) to the polynya flux model of Willmott et al. (1997) in an idealised basin.In an idealised setup, comparison with measurements is not possible and so the polynya flux model was used to validate the dynamic-thermodynamic model results.
Here we expand on the work done by Bjornsson et al. (2001) and compare the granular model to the more common viscous-plastic model of Hibler (1979) and the lesser known modified Coulombic yield curve by Hibler and Schulson (2000) in a setting similar to that used by Bjornsson et al. (2001).The granular model results are used to assess the outcome from the other two yield curves.Secondly, we consider formulations by Hibler (1979) and Mellor and Kantha (1989) for the thickness of newly formed ice.The former formulation was used by Kern et al. (2005) and Marsland et al. (2004) and the latter by Smedsrud et al. (2006).Finally, we use the collection depth parametrisation of Winsor and Björk (2000) to parametrise the new-ice thickness.Thus we address the important points of a polynya simulation; first the behaviour of the consolidated ice, which is determined by the rheology, and secondly ice formation inside the polynya, determined by the new-ice thickness parametrisation.
The layout of this paper is as follows: in Sect. 2 we discuss polynya formation and how to interpret the results of dynamic-thermodynamic models in light of what we know about polynyas.That section also includes a short description of the dynamic-thermodynamic model and a presentation of the control run by which the following experiments are assessed.In Sect. 3 the response of the model using different yield curves is presented.Section 4 presents the effects different formulations for the thickness of newly formed ice have on the model results.Section 5 contains a discussion of the model results followed by the conclusions of this study.

Polynyas in a dynamic-thermodynamic model
We will now discuss wind-driven polynyas and how they are modelled using dynamic-thermodynamic sea-ice models.The discussion focuses on understanding the processes involved in polynya formation and how to relate those to the results of the dynamic-thermodynamic model.We find that when it comes to understanding the model behaviour inside the polynya it helps to keep some of the assumptions of the flux polynya models in mind.Polynya flux models are simplified physical models which underline the important processes in polynya formation.They have been proven to be useful in simulating a variety of situations (see e.g.Morales Maqueda et al., 2004).In Subsects.2.1 and 2.2 we briefly describe the dynamic-thermodynamic model and the control run used to assess other model results.
Wind-driven coastal polynyas form where the ocean is initially covered by ice and a wind starts blowing off the coast.This causes the ice to move off-shore, opening a polynya at the coast (or fast ice edge).Inside the polynya the ocean is at the freezing point causing frazil ice to form and be herded downstream by the wind and waves.The frazil ice then consolidates at the edge of the initial ice.The polynya remains open as long as the off-shore wind component remains strong enough to maintain it.
The ice in and near a polynya can be divided into three distinct regimes: The thick initial ice, the consolidated ice and the frazil ice inside the polynya.The polynya edge is the interface between the polynya and the consolidated ice.This threefold separation is the basis of flux polynya models.They calculate the location of the polynya edge based on the drift velocity of the consolidated (and thick) ice, the ice formation rate inside the polynya and the thickness of the consolidated ice (H , also referred to as collection depth).
In the first flux polynya model, proposed by Pease (1987), the frazil ice inside the polynya is immediately transferred to the edge of the consolidated ice, where it piles up.The model by Ou (1988), and later models, assume a constant (and finite) velocity for the frazil ice, but this must always be greater than the velocity of the consolidated ice.In reality the frazil ice drifts faster than the consolidated ice because frazil ice, near or at the surface, experiences less water shear stress than the consolidated ice.The water velocity inside the polynya is also different from that under the consolidated ice, but this can be difficult to account for in a simplified setup.Finally the initial ice pack may not drift at the (local) free drift speed, as the wind that creates the polynya is nonuniform and may be weaker further off shore.Islands and other coast lines may also slow down the drift of the initial ice pack.
Polynya flux models focus on the frazil ice representation and the parametrisation of the collection depth.The velocity of the thick initial ice must, however, always be prescribed.Dynamic-thermodynamic models are, on the other hand, designed to model the movement of this thick ice, but they generally do not include any frazil ice parametrisations.When ice forms over open water in dynamic-thermodynamic models solid ice of a predetermined thickness, h 0 (see Sect. 2.1), is immediately formed, more akin to pancake ice than frazil ice.
In a polynya modelled by a dynamic-thermodynamic model the ice in the polynya interior drifts towards the initial ice pack, forming the consolidated ice.The consolidated ice consequently has thickness close to h 0 , which, in this particular setup, is then effectively the model's "collection depth".Assuming that the ice in the polynya interior drifts faster than the consolidated ice, like the Ou (1988) model demands, the dynamic-thermodynamic model should also show the three ice regimes the polynya flux models do (see Fig. 1).In addition Bjornsson et al. (2001) showed that the transition from freely drifting ice in the polynya interior to consolidated ice can occur over a few grid cells.This transition region is then analogous to the polynya edge predicted by flux polynya models.
We also note that since the consolidated ice and the ice in the polynya interior are modelled using the same drag coefficients their free drift speed will be the same.Absent any other forcing, or a divergence in the wind forcing, no sharp polynya edge will form.This is because the polynya that opens up fills with ice drifting at the same speed as the consolidated ice, resulting in linearly increasing ice concentration inside the polynya.Bjornsson et al. (2001) noted this behaviour (see their Fig. 4).In their study a polynya edge forms because the drift of the consolidated ice is slowed down by one of the side walls of their ideal basin.This approach is also used here.
A final point here is that we assume the polynya being modelled to be large enough to cover a substantial number of model grid points.In particular we demand that the model resolve the three ice regimes and the polynya edge.A single grid cell not completely covered with ice may often be interpreted to contain a polynya, especially if the grid resolution is low.In this study, however, we only consider polynyas that are properly resolved by the model grid.

The dynamic-thermodynamic model
The ice model is a two class (ice and open water) dynamicthermodynamic sea-ice model which was written to be coupled with the VOM ocean model (Backhaus, 2008).In this paper the ice model is coupled to a stationary slab ocean.In the present context the most important points to discuss regarding the model are the facts that one can choose between three different viscous-plastic rheologies and two different ways in representing the newly formed ice.In this section we briefly describe the model focusing on the rheology and new-ice thickness formulations.
The ice is modelled as a continuum using an Eulerian perspective.It moves in a horizontal plane, subject to both external and internal forces.Temporal evolution of the seaice cover is described using two continuity equations and the momentum equation.The continuity equation for mass is where m is the sea-ice mass per unit area, S m is a thermodynamic source/sink term and v is velocity.The source/sink term is a simple mass conservation formula: where ρ i is the ice density, A the fractional ice concentration, h represents the thickness changes in the ice already present in the grid cell and ρ i (1−A) h ow represents mass addition through ice formation over open water.An equation for the evolution of the ice thickness distribution within each cell is also needed.equation of conservation of sea-ice concentration.That takes the same basic form; with S A being a source/sink term.The average ice thickness over ice-covered area, h can be derived using m=Ahρ i .The source/sink term S A will be discussed below.In addition the condition A≤1 is imposed.This can be interpreted as a ridging condition since m (and thus h) can increase even if A does not (Hibler, 1979).Together these equations describe the advection of the ice in a given velocity field.
The momentum equation used is (Coon et al., 1974) Here k is a unit vector normal to the surface, τ a and τ w are air and water stresses, f is the Coriolis factor, g is the gravitational acceleration, H is the sea surface height and σ is the sea-ice stress tensor.Of the forcing terms on the right hand side ∇•σ describes forces due to internal stress while the other terms are all external factors.The term on the left hand side is the material derivative; D/Dt=∂/∂t+v•∇.It can in most cases be safely ignored (Rothrock, 1975) where C dw and C da are drag coefficients, ρ w and ρ a are air and water densities and v w and v a are the near surface water and wind velocities.This assumes that the wind velocity is much larger than the ice velocity; i.e. that |v−v a |≈|v a |.
The last term in Eq. ( 4) is the force due to the divergence of the internal stress tensor σ .Stress and strain rate (and thus ice velocity) are related through the sea-ice rheology.The three rheologies tested here are those suggested by Hibler (1979), Tremblay and Mysak (1997) and Hibler and Schulson (2000), which are all isotropic viscous-plastic rheologies.
In isotropic viscous-plastic models the relation between stress and strain can be represented using the stress invariants where ζ and η are the non-linear bulk and shear viscosities, εI and εII are the strain rate invariants and P is a pressure term.The viscosities, ζ and η, depend on the strain rate and P .We briefly describe each rheology here, but refer the reader to the original publications for more detail.
The rheology suggested by Hibler (1979) is by far the most commonly used viscous-plastic rheology.It uses an elliptical yield curve (see Fig. 2) and a normal flow rule.For typical strain rates, plastic flow occurs, while for very small strain rates the ice becomes linear viscous and enters so called "creep flow".This creep flow occurs for all three yield curves at small strain rates.The pressure term is calculated via where P * and C are constants, which must be chosen empirically.The elliptical yield curve reproduces basic sea-ice characteristics, i.e. the ice is weak in tension, strong in shear and strongest in compression.
A different approach to the elliptic yield curve was suggested by Tremblay and Mysak (1997), who proposed a model based on a granular material rheology.This is equivalent to dynamic friction between two dry surfaces where the frictional shear force is proportional to the normal force.For stress ratios of shear force vs. normal force smaller than a given coefficient of friction the ice behaves like an elastic solid, while for larger stress ratios it flows like a fluid.The pressure term is found by perturbing the last known solution to the momentum equation so that the resulting velocity field has a divergence determined by a given angle of dilation.This is done via an additional iterative solver and the numerical performance of the granular model is therefore worse than that of the other two formulations.For pressure larger than the maximum, determined by Eq. ( 8), the ice compresses but otherwise behaves like an incompressible Newtonian fluid with non-linear shear viscosity.The bulk viscosity (ζ ) is therefore always zero.For the granular model P * in Eq. ( 8) is replaced with P * /1.5 in accordance with Tremblay and Hakakian (2006).The resulting yield curve is a triangle (see Fig. 2).
Finally we include a yield curve which can be thought of as a combination of the other two.This is the so-called modified Coulombic yield curve, introduced by Hibler and Schulson (2000).Even though the term "modified Coulombic yield curve" (coined by Heil and Hibler, 2002) is used here, this yield curve is really a modification of the elliptic yield curve proposed by Hibler (1979).The curve gives friction-based failure up to a limiting compressive stress while the elliptic formulation holds for higher stresses.This limit is set at pure shear deformation, in accordance with the results of laboratory experiments.The pressure term is calculated using Eq. ( 8).The yield curve also includes a small amount of tensile stress.The resulting shape is elliptical for large stresses and triangular for small stresses (see Fig. 2).
The other main focus of this study is on the formulation of the new-ice thickness, which is expressed through the source/sink term S A in Eq. ( 3).The two methods for calculating S A considered, were formulated by Hibler (1979) and Mellor and Kantha (1989).
Hibler's approach can be written as where h 0 is a demarcation thickness separating thick and thin ice.Hibler (1979) assumed a constant h 0 , but in Sects.4 and 5 we discuss ways to parametrise it.Ice forms over open water (i.e.h ow >0) when the ocean cools below the freezing point, but the ocean temperature is calculated via the thermodynamics routines from VOM (which are the same routines as in Harms et al., 2003).The amount of ice formed over open water is calculated from the energy needed to bring the The Cryosphere, 4, 147-160, 2010 www.the-cryosphere.net/4/147/2010/ocean surface from its super-cooled state to freezing.When freezing h ow >0 and S m >0 so that Eq. ( 9) becomes simply This means that the newly formed ice covers an area S A so that its thickness is greater than or equal to h 0 .The other approach to calculating S A considered here is the one originally proposed by Mellor and Kantha (1989).
Based on an empirical formula by Nikiferov (1957) they formulate S A as where is an empirically determined function.Mellor and Kantha (1989) differentiate between melting and freezing by giving different constant values; We can easily recover Eq. ( 10) from Eq. ( 11) by setting F =h/ h 0 , which gives h 0 =h/ F .Equation ( 11) therefore states that during freezing the newly formed ice will have a thickness equal to h/ F .In particular, this means that when ice forms in a grid cell that previously had no ice (i.e.h=0) this cell will become fully covered with ice with thickness h ow .

Control run
The model domain is a bay, 135 km long and 75 km wide at 2.5 km resolution (see Fig. 3), similar to the setup Bjornsson et al. (2001) used.At such a high resolution we must assume that on average the ice floes being modelled are no larger than 250 m in diameter.This is because a scale of approximately 10 grain widths can generally be modelled without resolving each individual element using a granular model (Savage, 1998).McNutt and Overland (2003) state that at the multifloe scale (approximately 2-10 km) sea ice behaves like a granular material, so the granular model should be ideal for a simulation at that scale.We can assume the ice floes in the polynya interior are no larger than pancake ice, which is not much larger than 3 m in diameter.The ice in the polynya interior, and by extension the consolidated ice, is therefore well within the maximum allowed floe size.However, assuming that individual floes are no larger than 250 m in diameter may not always be valid for the thick initial ice, depending on the geographical location of the polynya as well as the time of year.This is only a minor concern for this study since it focuses on the steady state solution where the thick initial ice does not play a role (as it has drifted out of the model domain).
A polynya is created by having a 15 m/s wind blow uniformly at a 30 • angle to the direction along the bay.The polynya forms at the bottom of the bay and the excess ice flows out the open boundary at the mouth of the bay.The atmospheric temperature is kept constant at T air =−20 • C and the oceanic temperature is kept at the freezing point for a salinity of S=32.The water velocity is always zero.The model is initialised with ice concentration A=0.9 and thickness h=1 m.For the solid boundaries, a no-slip condition is used, while for the open boundary zero gradient Neumann boundary conditions are applied to all variables.The Neumann condition is also used for the ice pressure P , which Bjornsson et al. (2001) set to zero at the open boundary.This is done because using the Neumann condition improves the model behaviour near the open boundary by eliminating the excessive ice speed observed there by Bjornsson et al. (2001).
For the ice pressure constants C and P * in Eq. ( 8) we use the same values Bjornsson et al. (2001) used; i.e.C=30 and P * =30 kN/m 2 .Bjornsson et al. (2001) showed that these parameters have little effect on the model results using the granular model and we have found the same to be true for the other yield curves.A list of the relevant constants is included in Table 1.
To illustrate the temporal evolution of the polynya, Fig. 4 shows a Hovmöller diagram of the ice concentration field taken along a section at y=37.5 km.The response to the applied wind stress is immediate and a discernible polynya edge starts to form during the first day of the model integration.After two days the polynya has a clear structure and can be considered fully formed.A practically steady state has been reached after eight days.When the polynya has fully formed a band of large gradient in the concentration field analogous to the polynya edge always exists.For further reference Fig. 5 shows the ice concentration in the basin after eight days of model integration.As Figs. 4 and 5 show, the edge in this simulation is at a concentration of between A=0.6 and A=0.9.
Ice formation rates in the model are closely linked to the fractional ice concentration.For open water the ice formation E. Ö. Ólason and I. Harms: Polynyas in an ice model rate is F (A=0)=13.83 cm/day.In the polynya interior the ice is between 30 cm and 34 cm thick.If we assume all that ice is 32 cm thick we can calculate the ice formation rate in the polynya interior as a weighted average of F (A=0) and F (A=1, h=32cm)=3.31cm/day; i.e.
This approximation is correct to within 0.05 cm/day for A 0.8, but starts breaking down as the consolidated ice gets thicker.Defining the polynya as all points for which A<0.8 (this choice will be discussed further below), the mean ice formation rate in the polynya is F =11.1 cm/day after two days and F =10.7 cm/day after eight days.
According to Ou (1988), ice velocity in the model should fall into two categories; that of free drift in the polynya itself and that of the consolidated ice.In the dynamicthermodynamic model this velocity change gives the ice drifting in the polynya interior a barrier of slower consolidated ice to pile up against.Figure 6 shows the velocity field and speed in the control experiment after eight days.
The speed does indeed fall into two categories: The free drift speed |v f |=32.6 cm/s and the speed of the consolidated ice |v c | 25 cm/s, depending on the distance away from the y=75 km boundary.More importantly the cross channel velocity, v, changes from v f =15.7 cm/s to v c 3 cm/s in the consolidated ice.Bjornsson et al. (2001) showed that it is the cross channel velocity that is most important for the size and shape of the polynya.In the polynya interior the ice drifts with the wind but the consolidated ice slides along the y=75 km boundary with a small cross channel velocity component due to ridging at the boundary.This setup exhibits the three fold separation of ice mentioned previously; free drift ice in the interior of the polynya, consolidated ice at the polynya edge and thick initial ice beyond that.This can be seen in the ice concentration field since the concentration is low (A 0.6) in the polynya interior and high (A≈1) in the consolidated ice.The transition between the two takes place over approximately 5 grid cells, a region referred to here as the polynya edge.A similar transition is seen in the velocity field; the ice in the interior is in free drift while the consolidated ice drifts slower.The transition between the consolidated ice and the thick initial ice can only be seen in the ice thickness field.The ice in the polynya interior is almost as thick as the consolidated ice so the polynya edge is not apparent in the ice thickness field.In this control run the polynya is most readily defined as an area of low ice concentration and the polynya edge as the area of a large ice concentration gradient.In the rheologies used here Eq. ( 8) dictates the ice strength under compression as a function of ice thickness and concentration.When this term is small the rheology term becomes small and the ice is in free drift.Since the ice in the polynya interior is in free drift we propose that the location of the polynya edge, in this particular setup, can be approximated by A=0.8.This is because P (A=0.8)/P(A=1)≈0.01so the rheology will play a negligible role at lower concentrations.This choice is valid for the control run since the A=0.8 contour is within the high gradient region where the polynya edge is found (see Fig. 5).It is also valid for all the following experiments done here, except when using a minimum on the bulk viscosity (ζ ) and when using the newice thickness parametrisation of Mellor and Kantha (1989).These two cases will be discussed separately.

Different yield curves
An important part of the motivation for this study was to compare the results of Bjornsson et al. (2001), using the granular model, to a similar setup using different rheologies.The elliptic yield curve of Hibler (1979) is the most popular yield curve in use today, but it was designed for much lower resolution.It is therefore important to see how it fares in this high resolution setup.The modified Coulombic yield curve was, on the other hand, designed to model ice at high resolution so it will be instructive to see its performance here as well.In this section the focus is on the model response in the consolidated ice, since the rheology does not play a role inside the polynya itself.
In his model Hibler (1979) used a minimum for the bulk viscosity; ζ min =4×10 8 kg/s "in order to insure against any nonlinear instabilities" noting that that value is "several orders of magnitude below typical strong ice interaction values and effectively yields free drift results" (Hibler, 1979, p. 823).The lower bound should only be necessary when the material derivative is included in the momentum Eq. ( 4) since in that case having no lower bound will typically result in a noisy solution (Griffies and Hallberg, 2000).Most modern ice models ignore the material derivative and do not include a lower bound on ζ .However, as the resolution increases the material derivative becomes more important.The scaling argument made by Rothrock (1975) shows that the material derivative may need to be included as the shortest significant length scale becomes smaller than about 5 km.
In this idealised setup the material derivative has very little effect on the final solution, regardless of which rheology is used.The difference between the results with and without the material derivative is about 7 mm/s in the grid cells at the x=0 boundary and about 2 mm/s at the y=75 km boundary and at the polynya edge.Everywhere else in the domain the difference is less than 0.1 mm/s.Compared to the free drift speed of |v f |=32.6 cm/s this is small.We also observe no noise, even when the material derivative is included and the lower bound on ζ is set to zero.
Without a lower bound on ζ the results using the elliptic yield curve are nearly identical to those using the granular model.There is a sharp transition from free drift ice to consolidated ice, analogous to the polynya edge, just like in the granular model.This polynya edge forms at nearly the same location as it does in the granular model.Using the modified Coulombic yield curve also gives results very similar to the granular model.This is to be expected, since the modified Coulombic yield curve is in a way a combination of the other two yield curves.The difference between the three model formulations is limited to a small variance in polynya size, which as Fig. 7 shows, is due to a shift of the polynya edge by a few grid boxes.In particular, the commonly used elliptic yield curve is sufficient and can be used safely in this context.
Using a minimum on ζ does, however, give considerably different results from the control run.Setting the minimum to ζ min =4×10 8 kg/s, like Hibler (1979) did, results in a polynya with a very diffuse edge, as Fig. 8 shows.Speed and velocity also fail to meet the criteria for forming a polynya edge; i.e. there is no clear separation between the velocity of the ice in the polynya interior and consolidated ice (see Fig. 8).In addition, the ice in the polynya interior does not flow at a constant speed and its speed is sometimes lower than that of the consolidated ice, which is clearly not plausible.The speed of the ice in the polynya interior when using ζ min =4×10 8 kg/s is also always lower than it is in the control run.A decrease in the ice concentration is also seen by the y=75 km boundary, contrary to what can be seen in the control run.These effects were noted by Hunke (2001) in a different setup, but that discussion focused on the effects seen at the solid boundary, which will not be discussed further here.The polynya edge becomes so diffuse because when imposing such a high minimum on ζ , the viscosity is consistently kept at its minimum value for A 0.9 throughout the simulation.This is because the viscosity is related to the ice pressure via and the pressure to ice concentration via Eq.( 8).The flow where A is sufficiently small is therefore linear viscous, but choosing as a minimum ζ min =4×10 8 kg/s does not result in effectively free drift.Lu et al. (1989) also found that this limit was too high compared to measurements.As previously stated, we observed no non-linear instabilities using ζ min =0 kg/s.Some noise is, however, to be expected in a realistic simulation and so a non-zero ζ min may be required if one wants to include the material derivative.In that case one would need to choose a low, but non-zero value for ζ min .The resolution of Hibler's model was 125 km and since viscosity scales with the distance squared, a choice of ζ min =4×10 4 kg/s seems in order.This yields nearly the same results as with ζ min =0 kg/s; the largest difference in concentration between the two model runs being A=0.006.The maximum concentration gradient when using ζ min =0 is max(|∇A|)=0.230km −1 and it is max(|∇A|)=0.228km −1 when using ζ min =4×10 4 kg/s.When choosing larger values for ζ min , the effects of the capping start to show.For ζ min =4×10 5 kg/s the maximum gradient is max(|∇A|)=0.212km −1 and the difference in concentration between that run and the one with zero ζ min is A=0.08.For ζ min =4×10 6 kg/s the maximum gradient is max(|∇A|)=0.164km −1 and the concentration difference is A=0.3.

New-ice thickness
As we have already seen, the ice rheology affects the initial ice pack and the consolidated ice.The interior of the polynya, on the other hand, is primarily affected by the newice thickness parametrisation.This determines the thickness, and thus the concentration of the ice formed inside the polynya.
The most popular method for parametrising the new-ice thickness is probably the one suggested by Hibler (1979) (see Sect. 2.1).Put simply the new-ice thickness is not allowed to drop below a certain minimum, h 0 .If the total mass of newly formed ice is not enough to cover the open water fraction of the grid cell at that thickness then the concentration of newly formed ice is adjusted accordingly.If more ice is formed then the new ice is simply thicker than h 0 .
The choice of h 0 is not obvious and appears to range from 10 to 50 cm or even more in some cases.Bjornsson et al. (2001) argued for using h 0 =30 cm and that is the value used here so far.Their argument is based on the assumption that the ice that forms in the polynya immediately forms pancake ice.However, Bjornsson et al. (2001) state that pancake ice thickness is closer to 10 cm than 30 cm.We therefore include a model run with h 0 =10 cm.
The main result of using h 0 =10 cm is that with a lower h 0 the polynya fills up much faster.The newly formed ice is thinner, therefore has a larger surface area, which results in faster ice concentration growth in the polynya itself and causes the polynya edge to form closer to the bottom of the bay than before.The polynya edge is also not as sharp in the |v|, bottom) using h 0 =10 cm after eight model days.The resulting polynya is smaller and has a higher ice concentration than the control run.The dash-dotted line shows the isoline for A=0.8 from the control run.concentration field as when using h 0 =30 cm.It is, however, still sharp in the velocity field, as can be seen in Fig. 9.
The mean ice formation rate in the polynya is about 1% lower here than in the control run.The total ice formation is therefore reduced almost only because the polynya is smaller.Using h 0 =30 cm gives polynya area of A p =4.5×10 3 km 2 after eight days, but using h 0 =10 cm the polynya area is A p =2.4×10 3 km 2 after eight days; an approximately 50% reduction in size.Finally, the consolidated ice is thinner since its thickness equals h 0 .
The other approach to determining the thickness of newly formed ice described in Sect.2.1 is the one proposed by Mellor and Kantha (1989).There the thickness of newly formed ice is based on the thickness of the ice already present in the grid cell.Mellor and Kantha (1989) argued that the thickness of newly formed ice should be a quarter of the old ice thickness.This also means that when there is no ice in the grid cell when new ice forms, the ice spreads uniformly over the entire cell, potentially very thinly.A polynya in such a model may therefore be hard to recognise by the change in concentration and researchers using this approach often consider ice below a certain cut-off thickness to represent the polynya.Smedsrud et al. (2006), for instance, use Ah=30 cm for this cut-off thickness.Fig. 10.The ice concentration (A, top) and speed and velocity (v and |v|, bottom) using the new-ice thickness parametrisation by Mellor and Kantha (1989) after eight model days.The resulting polynya is very small with no discernible edge in the velocity field.The dash-dotted line shows the isoline for A=0.8 from the control run.
Using this approach results in a "polynya" that is hardly recognisable in the concentration field, as Fig. 10 shows.Even after eight days there is only a thin sliver of an opening along the x=0 km and y=0 km boundaries and the A=0.8 isoline is a grid box or two away from the shore.More seriously perhaps the velocity field, also shown in Fig. 10, shows no sign of the discontinuity deemed necessary for proper polynya formation.The ice slows down gradually moving away from the x=0 boundary, contrary to our previous assumptions about how a polynya is formed.
Ice thickness near the coast is indeed lower than the initial ice thickness, but as Fig. 11 shows there is no real polynya edge to be found in the ice thickness field.In that respect there is no conceptual difference between using h or Ah.As before the thick initial ice drifts out of the basin, but in this case the ice that replaces it does not have a uniform thickness.It is very thin at the coast with linearly increasing thickness towards the thick initial ice.
We have already mentioned that considerable effort has been put into parametrising the collection depth in polynya flux models.Given the large variation between the results already presented in this section we find it worth considering whether the polynya flux model parametrisations can be applied in the dynamic-thermodynamic model.Mellor and Kantha (1989) taken along a section at y=37.5 km.The vertical axis (x) is the along channel distance and the horizontal axis (t) is the model time.The dash-dotted line shows the isoline for A=0.8 from the control run.The polynya edge after eight days using a constant h 0 and parametrised h 0 according to Eq. (15).A=0.8 is used as a marker for the polynya edge and the edge is plotted for 10, 20 and 30 m/s wind speed.
The parametrisation by Winsor and Björk (2000) lends itself well to immediate inclusion in the dynamicthermodynamic model.It is based only on the wind speed and not the polynya width, frazil ice speed or other quantities not accessible to the dynamic-thermodynamic model.By using this parametrisation we aim at improving the modelled consolidated ice thickness and thus also the size of the polynya.Winsor and Björk (2000) assumed the collection depth to be a function of wind speed as where |v a | is the surface wind velocity, a=1 m, b=0.1 s and c=15.In particular, H ≈7 cm for |v a |=0 and H =30 cm for |v a |=35 m/s so this parametrisation is well within the range of plausible values for h 0 .Equation ( 15) is then used to calculate h 0 =H in each grid point.
Using this parametrisation results in smaller polynyas at low wind speeds, compared to h 0 =30 cm or larger polynyas at high wind speeds, compared to h 0 =10 cm. Figure 12 shows the polynya using h 0 =30 cm and the parametrisation for the wind speeds 10, 20 and 30 m/s.At lower winds the polynya edge starts to become diffuse, which is to be expected.For further reference Fig. 13 shows the size of the resulting polynya as a function of wind speed.The mean ice formation rate in the polynya increases from F =10 cm/day for |v a |=10 m/s to F =11 cm/day for |v a |=35 m/s and choosing different values for h 0 contributes to about 1% change in the ice formation rate.At the same time the polynya size grows approximately four times as the wind strength grows from 10 m/s to 35 m/s.It is therefore clear that variations in polynya size control the variations in total ice formation in the polynya.Bjornsson et al. (2001) have already shown that the granular model can be used to model polynyas in an idealised setting by comparing their results to a polynya flux model.We have shown that this is also the case when using the modified Coulombic yield curve of Hibler and Schulson (2000) and when using the elliptic yield curve of Hibler (1979).This is important since the numerical performance of the granular model is considerably worse than that of the other two models.We find that the granular model requires almost twice the computing time the other two rheologies require.The elliptic yield curve of Hibler (1979) is also already in use in the vast majority of sea-ice models.

Discussion
All three yield curves give nearly identical results (with the exception of using a capped ζ as discussed below).Looking at the stress states (Fig. 14) we see that when using the granular model the σ I values for points in the consolidated ice are all clustered around σ I =−P /1.5.This means that at these points the ice cover yields or is very close to yielding under compression.When this is the case the behaviour of the ice is controlled by Eq. ( 8), which is also one of the main equations governing the behaviour of the other two yield curves.The difference between the granular model and the other two is then almost entirely explained by the different formulation of shear strength between the model formulations.As Fig. 2 shows, the granular model and the ellipse have a lower shear strength than the modified Coulombic yield curve, which is why using the modified Coulombic yield curve gives a smaller polynya.
Using the other two yield curves, the stress states are much more evenly distributed along the σ I axis.For the modified Coulombic yield curve the stress states that lie on the Coulombic slope are all inside the polynya while the stress states in the consolidated ice are all on the elliptic part of the yield curve.This happens because in the consolidated ice the divergence (ε) is always negative, so according to Eq. ( 7) σ I ≤−P /2.Using the elliptic and modified Coulombic yield curves therefore yields similar results for the consolidated ice, where both yield curves have an elliptic shape.In the polynya interior the ice is in free drift so the shape of the yield curve has no effect there.
Using a large ζ min , in particular ζ min =4×10 8 kg/s as Hibler (1979) suggests, does, however, give results that are not plausible.Using this original formulation results in a polynya that is smeared out with no proper edge and a velocity field that has little relation to polynya formation.This happens because the capping of ζ in the model turns the viscous-plastic formulation into a linear viscous model for ice concentration A<0.9.
Hibler (1979) used the lower bound on ζ to dampen grid scale noise which can arise when the grid Reynolds number constraint is not satisfied.Considering the limited case of Burger's equation in one dimension: A must be bounded by A> 1 2 v x (Griffies and Hallberg, 2000).In free drift and ignoring the sea surface tilt term, the momentum equation becomes Burger's equation and for the one dimensional case A = ζ /ρ i .Using the free drift speed of less effects on the simulation than using a non-zero ζ min we conclude that ignoring the material derivative is preferable to including it and a non-zero ζ min .
With regards to the ice behaviour inside the polynya we considered three ways in which to parametrise the thickness of ice forming over open water.These are the methods suggested by Hibler (1979), Mellor and Kantha (1989) and an adaptation of the collection depth parametrisation by Winsor and Björk (2000).Hibler's method was used when investigating the dynamic aspects (Sect.3) with an ice demarcation thickness h 0 =30 cm after Bjornsson et al. (2001).That value may be too high and so we also ran the model using h 0 =10 cm.
Using a lower h 0 resulted in a smaller polynya but little change in mean ice formation rates (about 1%).The ice concentration in the polynya interior was higher and as a result the concentration field did not show a sharp polynya edge.The velocity field, on the other hand, still showed a clear discontinuity at the polynya edge.The width of the polynya did decrease, but that was to be expected and can be understood in relation to the Lebedev-Pease width of a polynya (Pease, 1987) www.the-cryosphere.net/4/147/2010/The Cryosphere, 4, 147-160, 2010 where L is the polynya width and HU is the flux of consolidated ice.The collection depth, H , is analogous to h 0 in the dynamic-thermodynamic model.Lowering h 0 from 30 cm to 10 cm results in approximately 1% reduction in the ice formation rate (F ).The polynya width is therefore bound do decrease.
Results obtained using the formulation of Mellor and Kantha (1989) were, however, very different from those obtained in the control run.Using Mellor and Kantha's formulation there are two ice regimes; thick ice and thin ice, which may be characterised as nilas.This replaces the threefold separation of thick ice, consolidated ice (of uniform thickness) and frazil/free drift ice, seen in polynya flux models and the control run.The polynya edge is considered to be what separates the consolidated ice and frazil/free drift ice, but this distinction is lost when using the Mellor and Kantha (1989) approach.The new-ice thickness formulation by Mellor and Kantha (1989) is therefore not suitable for modelling polynyas.
On the whole, the approach by Hibler (1979) is also more reasonable from a physical standpoint.This is because wind and waves, which cannot be resolved by ocean or atmosphere models, will transform the frazil ice in the polynya into pancake ice, similar to the ice formation in that scheme.What is unrealistic about Hibler's approach is that solid ice forms inside the polynya, even where in reality the ice is mainly frazil ice.The thin ice formed using Mellor and Kantha's approach is more akin to grease ice or nilas which form in calmer conditions.
Wind speed is therefore an important factor in determining the new-ice thickness and it is consequently an important part of collection depth parametrisations for polynya flux models.Winsor and Björk (2000) parametrised the collection depth in the Pease (1987) model based only on wind speed and we found that parametrisation easily adoptable for inclusion in the dynamic-thermodynamic model.
Using the Winsor and Björk (2000) parametrisation gives results in the range between the results when using a constant h 0 =10 cm and h 0 =30 cm.We have already expressed a preference for the Hibler (1979) parametrisation for the new-ice thickness and using the Winsor and Björk (2000) parametrisation enables us to choose a sensible value for h 0 .As a result the polynya size should depend on wind strength in a more realistic manner than when using a constant h 0 .Using the parametrisation should also give more realistic ice thickness for the consolidated ice.
On a more general note, such a small value for h 0 may not be suitable for models describing the central pack ice as well.In such a situation the approach of Mellor and Kantha (1989) may give better results since the thick pack ice appears to require a larger h 0 .
It is trivial to combine all three approaches to new-ice thickness parametrisation discussed here into one: with =4, a=1 m, b=0.1 s and c=15, as before.This approach modifies the previously constant h 0 of Hibler (1979) so that for thick ice the approach of Mellor and Kantha (1989) is used and for thinner ice the parametrisation of Winsor and Björk (2000) is used.Using one value for thick ice and one for thin is appropriate since h 0 is only analogous to the collection depth in a polynya or the marginal ice zone.
Where the ice is thicker the ice behaviour should be similar to that described by Mellor and Kantha (1989), given the empirical origin of their formulation.
In our setup the result of Eq. ( 18) will always be the same as that of Eq. ( 15) since the ice in the polynya is always thin in the sense that h< (a+|v a |b)/c.Further testing of this new parametrisation can therefore not be done here but should be carried out in a realistic simulation.
In conclusion we note that in this idealised setup the polynya is best defined as the area where A<0.8.This is really the concentration where the ice behaviour starts changing from free drift, for A 0.8, to being heavily influenced by internal stresses, for A≈1.This separation is based on Eq. ( 8) and on the choice of C. The A=0.8 isoline is also consistently within the high gradient region for A in all experiments, except when using too high a minimum for the bulk viscosity and when using the new-ice thickness parametrisation from Mellor and Kantha (1989).But both these cases were found to give implausible results.

Conclusions
We have used an idealised setup to test three different sea-ice rheologies and three different formulations for the thickness of newly formed ice during polynya formation.These tests were done using a dynamic-thermodynamic sea-ice model in an idealised channel, similar to what Bjornsson et al. (2001) did.
We were able to reproduce the results of the granular model using both the modified Coulombic yield curve of Hibler and Schulson (2000) and the elliptic yield curve of Hibler (1979).This is important, since the numerical performance of the granular model is substantially worse than that of the other two models and also since the elliptic yield curve is already in popular use.We also found that including the material derivative and setting a minimum on the bulk viscosity is not a viable alternative to ignoring the material derivative.
The formulation of new-ice thickness suggested by Hibler (1979) turned out to give much better results than that of Mellor and Kantha (1989).Using Mellor and Kantha's formulation failed to give a clear polynya edge, both in the concentration and velocity field.We conclude therefore that this approach does not enable us to properly model polynyas.Hibler's approach, on the other hand, gave a clear separation of the consolidated ice and the polynya itself.Using Hibler's new-ice thickness parametrisation and any of the rheologies tested here should give realistic results when modelling polynyas.Polynyas that are fully resolved by the model grid can then be recognised as areas of low concentration enclosed by compact ice and/or land.We also suggest using A<0.8 as a criterion for low concentration.Hibler (1979) assumed a constant demarcation thickness (h 0 ).We suggest, however, using the collection thickness parametrisation of Winsor and Björk (2000) to parametrise h 0 .This results in a value for h 0 which is dependent on wind strength and in the range already deemed acceptable for h 0 .As an aside, a combination of this parametrisation with the approach of Mellor and Kantha (1989) is proposed.This should give a parametrisation for h 0 applicable for both the marginal ice zone and the central ice pack.
use a combination of ice concentration and thickness and Marsland et al. (2004) use a combination of ice concentration and freezing rate.

Fig. 3 .
Fig. 3.The size of the basin (in km) and wind direction during the polynya experiments.The figure also shows the three ice regimes one expects; the polynya interior (where the ice is in free drift), consolidated ice and thick initial ice.

Fig. 4 .Fig. 5 .
Fig. 4. A Hovmöller diagram of the ice concentration field (A) in the control experiment taken along a section at y=37.5 km.The vertical axis (x) is the along-channel distance and the horizontal axis (t) is the model time.The dash-dotted line shows the h=1 m isoline separating the thick initial ice and consolidated ice.

Fig. 6 .
Fig. 6.Ice velocity (v) and speed (|v|) in the control experiment after eight days of model integration.The polynya edge is visible as a sharp decrease in the ice velocity.

Fig. 7 .
Fig. 7.The polynya edge (A=0.8) using the elliptic and modified Coulombic yield curves and the granular model after eight model days.The differences between different model formulations are minor.

Fig. 8 .
Fig. 8. Sea-ice concentration (top) and speed and velocity (bottom) using Hibler's original formulation for the elliptic yield curve, after eight days of model integration.Neither figure shows a discernible polynya edge.The dash-dotted line shows the isoline for A=0.8 from the control run.

Fig. 9 .
Fig. 9.The ice concentration (A, top) and speed and velocity (v and

E
Fig. 11.A Hovmöller diagram of the ice thickness field (h) using the new-ice thickness formulation byMellor and Kantha (1989) taken along a section at y=37.5 km.The vertical axis (x) is the along channel distance and the horizontal axis (t) is the model time.The dash-dotted line shows the isoline for A=0.8 from the control run.
Fig.12.The polynya edge after eight days using a constant h 0 and parametrised h 0 according to Eq. (15).A=0.8 is used as a marker for the polynya edge and the edge is plotted for 10, 20 and 30 m/s wind speed.

Fig. 13 .
Fig. 13.The polynya size (A p ) after eight days as a function of wind speed (|v a |).The size is calculated as the sum of the size of all model points where A<0.8.

Fig. 14 .
Fig. 14.Stress states using the granular model (top), the elliptic yield curve (centre) and modified Coulombic yield curve (bottom) after eight model days plotted in stress invariant space.The colour scale indicates the ice concentration at each point (A) and points with high concentration are drawn on top of those with a lower concentration.Points with high values of A (A 0.8) are all located in the region σ I ≤−P /2.Points near the open boundary are excluded from the figure.
The model uses two ice classes; i.e. ice and open water and so this becomes an The three yield curves discussed in the text plotted in stress invariant space; the elliptic, the modified Coulombic and the granular model's yield curves, depicted as solid, dashed and dotted lines, respectively.

Table 1 .
Main physical parameters and constants used in the simulation.