2.1 Stefan's ice growth formulation

Stefan's (Reference Stefan1889) approach to thermal ice growth considers a single ice layer bounded by the ice/air and the ice/water interfaces. The air/ice interface temperature is assumed to equal the air temperature T _A and the ice/water interface temperature is assumed to equal the freezing temperature of water T _F. Therefore, the temperature difference across the ice thickness creates a temperature gradient, which induces heat to flow from the higher temperature T _F towards the lower temperature T _A. The heat is transferred through the ice and snow medium with conduction and in air and water mostly with convection. Convection is the heat flux transferred by the molecular or bulk motion of a fluid. In snow, the high porosity and air content can cause, in addition to conduction, heat transfer by convection. The air temperature T _A is in practice the temperature outside the thermal boundary layer above the ice or snow surface. Thus, T _A is the screen level air temperature measured near the surface from meteorological stations. For example, the height from the ground is taken as 1.25 to 2 m with a tolerance of 0.05 m in ISO, 17714:2007-07.

Stefan's model neglects both convective fluxes by equalising the air temperature with the upper interface temperature and omitting the oceanic heat flux from the model. The ice formation is assumed to occur only vertically, therefore the horizontal dimension is omitted resulting in a 1-D growth model. The thermal growth problem is further simplified by assuming that the equilibrium temperature distribution is reached instantly, i.e. ∂T(z)/∂t = 0 where T(z) is the temperature distribution in ice, z is the vertical coordinate with an origin in the ice/atmosphere interface and t is the time. This means that the heat released during the ice temperature decrease due to thermal inertia is disregarded. Therefore, the assumption of a linear temperature gradient leads to Fourier's law of steady-state conduction of the heat (Incropera and others, Reference Incropera, Dewitt, Bergman and Lavine2007). The heat conducted through the ice layer to the atmosphere φ _I is expressed as (in units of W m⁻²):

(1)$$\varphi _{\rm I} = {-}k_{\rm i}\displaystyle{{{\rm d}T} \over {{\rm d}H}} = {-}k_{\rm i}( {T_{\rm A}-T_{\rm F}} ) /H_{\rm I}, \;$$

where k _i is the thermal conductivity of the solid ice (unit: W m⁻¹ °C⁻¹) and ΔT = (T _A − T _F) is the temperature gradient between the upper and lower interfaces of ice, the thickness of which is equal to H _I. The thermal conductivity of sea ice for temperatures above −8.2°C is mainly affected by the salinity and for lower temperatures by the sea-ice density, which depends on the air content (Schwerdtfecer, Reference Schwerdtfecer1963). Pure ice conductivity k _pi in W m⁻¹ K⁻¹ as a function of temperature in K is given by Yen (Reference Yen1981) as:

(2)$$k_{{\rm pi}} = 9.828{\rm e}^{( {-0.0057T} ) }.$$

From Eqn (2), the conductivity coefficient for pure ice at the freezing temperature of 273.15 K is equal to 2.07 W m⁻¹ K⁻¹. Additionally, the proposed thermal conductivity equations for the so-called ‘bubbly’ ice k _bi, which is the conductivity of ice that contains air inclusions, and sea ice k _i, which is the conductivity of ice that contains air and brine inclusions are given by Schwerdtfecer (Reference Schwerdtfecer1963) as:

(3)$$k_{{\rm bi}} = \displaystyle{{2k_{{\rm pi}} + k_{\rm a}-2v( {k_{{\rm pi}}-k_{\rm a}} ) } \over {2k_{{\rm pi}} + k_{\rm a} + v( {k_{{\rm pi}}-k_{\rm a}} ) }}k_{{\rm pi}}, \;$$

(4)$$k_{\rm i} = k_{{\rm bi}}-( {k_{{\rm bi}}-k_{\rm b}} ) \displaystyle{{S\rho _{\rm s}} \over {xT_{\rm i}}}.$$

Sea-ice thermal conductivity depends on the bulk salinity S and the ice temperature T _i. The coefficient of proportionality x between the fractional salt content in the ice mass and temperature is equal to −0.0182 °C⁻¹ and is extracted from the sea-ice phase diagram (Assur, Reference Assur1960). k _a and k _b are air and brine conductivity coefficients, respectively, while ν is the ratio between air and total volume of the air–ice matrix. Brine conductivity k _b (kcal m⁻¹ °C⁻¹) as a function of temperature (°C) is given as:

(5)$$k_{\rm b} = ( {1.25 + 0.03T_{\rm i} + 0.00014T_{\rm i}^2 } ) \times 10^{{-}2}.$$

Alternative sea-ice conductivity formulations are given by Untersteiner (Reference Untersteiner1964), and also by Pringle and others (Reference Pringle, Eicken, Trodahl and Backstrom2007).

The ice growth ΔH _I in the time interval Δt is estimated from the continuity of the heat fluxes in the ice/water interface. Therefore, the heat released from the ice formation φ _F in the ice/water interface (in W m⁻² or usually expressed as kJ m⁻² d⁻¹), is assumed to equal the heat flux conducted through the ice layer φ _I into the atmosphere and which starts from the lower surface of the ice. The heat flux continuity requires:

(6)$$\varphi _{\rm I} = k_{\rm i}\displaystyle{{( {T_{\rm F}-T_{\rm A}} ) } \over {H_{\rm I}}} = \rho _{\rm i}L_{\rm i}\displaystyle{{{\rm d}H_{\rm I}} \over {{\rm d}t}}, \;$$

where ρ _i and L _i are the ice bulk density in kg m⁻³ and latent heat of formation in kJ kg⁻¹ (or W d kg⁻¹) respectively. From the initial condition of the ice thickness H ₀ = 0 at t = 0, the integration of continuity Eqn (6) leads to the ice thickness formula:

(7)$$H_{\rm I} = \sqrt {\displaystyle{{2k_{\rm i}} \over {\rho _{\rm i}L_{\rm i}}}\theta } , \;$$

where cumulative freezing degree-days (variable θ) is calculated as:

(7a)$$\theta = \int_{t_0}^t {( T_{\rm F}-T_{\rm A}) } .$$

Stefan's constant, or the ice growth rate coefficient (e.g. in units of mm °C^−0.5 d^−0.5) is:

(7b)$$\alpha = \sqrt {\displaystyle{{2k_{\rm i}} \over {\rho _{\rm i}L_{\rm i}}}} .$$

In Eqn (7), the growth rate coefficient omits all the uncertainties of the model and neglected parameters such as the snow layer, atmospheric coupling, oceanic heat flux, net radiation heat flux, wind action, etc. Therefore, the growth rate varies significantly, as shown earlier. Comfort and Abdelnour (Reference Comfort and Abdelnour2013) estimated the growth rate using 436 ice thickness measurements. The growth rate coefficient was determined to be in a range from 1.1 up to 31.4 mm °C^−0.5 d^−0.5 with an average value of the growth rate equal to 17 mm °C^−0.5 d^−0.5. The theoretical value of the growth rate coefficient for pure ice is 34.12 mm °C^−0.5 d^−0.5 for k _pi determined from Eqn (2), with ice density equal to 917 kg m⁻³ and latent heat equal to 3.88 W d⁻¹ kg⁻¹.

2.2 Atmospheric ice/air coupling

Later studies introduced the atmospheric air/ice coupling into Stefan's approach (see e.g. Ashton, Reference Ashton1986, Reference Ashton1989; Ettema and Huang, Reference Ettema and Huang1990; Andres and van der Vinne, Reference Andres and van der Vinne2001). The difference between the air/ice interface temperature T _IA and the air temperature T _A depends on the convective heat flux coefficient h _ia. The air/ice interface temperature is estimated from the heat flux continuity. Therefore, the convective heat flux that is flowing outwards from the ice surface towards the air is equal to the ice conductive heat flux. This is generated from the temperature gradient between ice/air and ice/water interfaces. The air/ice interface temperature T _IA is obtained from the heat flux continuity at the upper boundary as follows:

(8)$$h_{{\rm ia}}( {T_{{\rm IA}}-T_{\rm A}} ) = k_{\rm i}\displaystyle{{( {T_{\rm F}-T_{{\rm IA}}} ) } \over {H_{\rm I}}}$$

where the ice/air interface temperature is obtained as:

(9)$$T_{{\rm IA}} = \displaystyle{{k_{\rm i}T_{\rm F} + h_{{\rm ia}}H_{\rm I}T_{\rm A}} \over {H_{\rm I}h_{{\rm ia}} + k_{\rm i}}}$$

The value of the air convective heat transfer coefficient h _ia varies between 10 and 25 W m⁻² °C⁻¹ and is usually assumed to depend on the wind speed (Ashton, Reference Ashton1986, Reference Ashton1989). The equation for heat flux continuity at the lower boundary is:

(10)$$\rho _{\rm i}L_{\rm i}\displaystyle{{{\rm d}H_{\rm I}} \over {{\rm d}t}} = {\rm \;}\;k_{\rm i}\displaystyle{{( {T_{\rm F}-T_{{\rm IA}}} ) } \over {H_{\rm I}}}$$

Substituting T _IA in Eqn (10) gives the ice growth rate expressed as:

(11)$$\displaystyle{{T_{\rm F}-T_{\rm A}} \over {H_{\rm I}/k_{\rm i} + 1/h_{{\rm ia}}}} = \rho _{\rm i}L_{\rm i}\displaystyle{{{\rm d}H_{\rm I}} \over {{\rm d}t}}.$$

Assuming that the initial ice thickness at time t ₀ is zero and at time t the ice thickness is equal to H _I, then the integration of the above expression yields the ice thickness equation at time t as:

(12)$$H_{\rm I} = \sqrt {{\left({\displaystyle{{k_{\rm i}} \over {h_{{\rm ia}}}}} \right)}^2 + \displaystyle{{2k_{\rm i}} \over {\rho _{\rm i}L_{\rm i}}}\theta } -\displaystyle{{k_{\rm i}} \over {h_{{\rm ia}}}}.$$

This ice thickness equation yields more accurate results for the initial ice growth compare to Stefan's (Reference Stefan1989) equation (Ashton, Reference Ashton1986, Reference Ashton1989). However, it neglects several components, such as the snow layer, radiation and thermal inertia.

2.3 Snow/air coupling

As long as the snow mass is smaller than the mass needed to depress the ice below the water level (Archimedes principle), snow remains dry and acts as insulation, decreasing the ice growth rate. The insulative effect on the ice growth and the snow/air coupling were included in Stefan's model by, for example, Nakawo and Sinha (Reference Nakawo and Sinha1981), Leppäranta (Reference Leppäranta1983), Ashton (Reference Ashton1986), Leppäranta (Reference Leppäranta1993), Andres and van der Vinne (Reference Andres and van der Vinne2001), Comfort and others (Reference Comfort, Gong, Singh and Abdelnour2003), Comfort and Abdelnour (Reference Comfort and Abdelnour2013). The two-layer model consists of three interface temperatures: T _F, snow–ice interface temperature T _IS and snow-air upper boundary temperature T _SA, the latter two being variables. The assumption of a linear temperature gradient through snow and ice permits the determination of the interface temperatures from the heat flux continuity at the interfaces. Thus, the snow–atmosphere interface temperature T _SA from the heat flux continuity at snow/air interface is expressed as:

(13)$$h_{{\rm sa}}( {T_{{\rm SA}}-T_{\rm A}} ) = k_{\rm s}\displaystyle{{( {T_{{\rm IS}}-T_{{\rm SA}}} ) } \over {H_{\rm S}}}.$$

This gives the snow–atmosphere interface temperature T_SA with the following expression:

(14)$$T_{{\rm SA}} = {\rm \;}\displaystyle{{k_{\rm s}T_{{\rm IS}} + h_{{\rm sa}}H_{\rm S}T_{\rm A}} \over {H_{\rm S}h_{{\rm sa}} + k_{\rm s}}}, \;$$

where h _sa is the air convective heat transfer coefficient in W m⁻² °C⁻¹. As mentioned in the previous section, the air convective transfer coefficient, h _ia depends mainly on wind conditions. The wind has a major effect on snow layer spatial thickness distribution, density and porosity. Considering that the effects of wind acting on the ice and snow layer are different, consequently, h _sa may be different in value from h _ia. The thermal conductivity of the snow layer k _s varies with the bulk density of snow. Empirical formulations, which may be applied in the growth models are given, for example, by Mellor (Reference Mellor1977), Sturm and others (Reference Sturm, Holmgren, König and Morris1997), Calonne and others (Reference Calonne2011) and Riche and Schneebeli (Reference Riche and Schneebeli2013). The most commonly used snow conductivity k _s (in units of W m⁻¹ °C⁻¹) in ice growth models is given as (Yen, Reference Yen1981):

(15)$$k_{\rm s} = 2.22362( \rho _{\rm s}) ^{1.885}, \;$$

where ρ _s is in units of Mg m⁻³. For a snow density ρ _s equal to 300 kg m⁻³, the snow thermal conductivity k _s is equal to 0.23 W m⁻¹ °C⁻¹.

The snow–ice interface temperature T _IS is determined by substituting the snow/air interface temperature T _SA (Eqn 14) into the heat flow continuity equation at the ice/snow interface:

(16)$$k_{\rm s}\displaystyle{{( {T_{{\rm IS}}-T_{{\rm SA}}} ) } \over {H_{\rm S}}} = {\rm \;}k_{\rm i}\displaystyle{{( {T_{\rm F}-T_{{\rm IS}}} ) } \over {H_{\rm I}}}.$$

This results in the following expression of the ice/snow interface temperature:

(17)$$T_{{\rm IS}} = {\rm \;}\displaystyle{{k_{\rm i}T_{\rm F}h_{\rm e} + T_{\rm A}H_{\rm I}} \over {H_{\rm I} + k_{\rm i}h_{\rm e}}}.$$

The parameter h _e was introduced by Riska and others (Reference Riska2019) and represents the equivalent heat transfer coefficient from the snow layer to the atmosphere and is given as:

(18)$$h_{\rm e} = {\rm \;}\displaystyle{{H_{\rm S}} \over {k_{\rm s}}} + \displaystyle{1 \over {h_{{\rm sa}}}}.$$

The latent heat of freezing is equal to the heat flux conducted through the ice and snow layer to the atmosphere:

(19)$$k_{\rm i}\displaystyle{{( {T_{\rm F}-T_{{\rm IS}}} ) } \over {H_{\rm I}}} = {\rm \;}\rho _{\rm i}L_{\rm i}\displaystyle{{{\rm d}H_{\rm I}} \over {{\rm d}t}}.$$

Substitution of the snow/ice interface temperature T _IS, calculated by Eqn (17), into Eqn (19) gives the growth rate expressed as:

(20)$$\displaystyle{{T_{\rm F}-T_{\rm A}} \over {H_{\rm I}/k_{\rm i} + h_{\rm e}}} = {\rm \;}\rho _{\rm i}L_{\rm i}\displaystyle{{{\rm d}H_{\rm I}} \over {{\rm d}t}}$$

Integration of Eqn (20) yields the ice thickness equation:

(21)$$H_{\rm I} = \sqrt {k_{\rm i}^2 h_{\rm e}^2 + \displaystyle{{2k_{\rm i}} \over {\rho _{\rm i}L_{\rm i}}}\theta } -k_{\rm i}h_{\rm e}.$$

In this formulation, a layer of dry snow decreases the rate of heat that is transferred from the lower surface of the ice to the atmosphere, causing the ice temperature to increase, even at very low air temperatures (Comfort and others, Reference Comfort, Gong, Singh and Abdelnour2003). For example, the maximum snow thickness of 140 mm observed in a prairie pond by Andres and van der Vinne (Reference Andres and van der Vinne2001) in the winter of 1992–93, reduced the maximum ice thickness by 37 mm compared to the bare ice growth in the same pond, which was maintained snow-free.

2.4 Snow ice

Apart from the insulation effect of the snow, the snow layer can also form a snow-ice layer. Snow-ice occurrence is mainly observed in the early winter, or at the time when the ice layer is relatively thin and heavy snowfall occurs, as well as at the end of winter when water is available from the melting and refreezing cycles of snow and precipitation.

The snow-ice formation in the early winter or from heavy snowfall occurs when the snow mass overcomes the buoyancy of ice and submerges the upper ice surface. In the presence of thermal or mechanical cracks, brine or air inclusions, the water travels to the ice/snow interface and forms a slush layer. The ice freeboard submerges, and the slush thickness increases until the hydrostatic equilibrium is restored. The slush layer thickness may also be affected by the water rising upwards in the snow caused by capillarity. The excessive amount of water driven by capillarity in the snow layer was found to depend on snow density and crystal size (Ager, Reference Ager1962).

After flooding and slush formation, surface ice growth begins. The surface ice growth is defined as slush transformation to snow ice. Snow-ice density varies usually between 850 and 900 kg m⁻³ (Ager, Reference Ager1962). The surface ice growth caused by flooding was predicted by Leppäranta (Reference Leppäranta1983). The snow thickness that depresses the ice layer below the water level is determined by comparing the weight of snow with the buoyancy of ice. The weight of snow is not enough to push ice under the water if the following is valid:

(22)$$\rho _{\rm s}H_{\rm s}\;\le {\rm \;}( {\rho_{\rm w}-\rho_{\rm i}} ) H_{\rm I}.$$

In this case, slush formation does not occur, the dry snow layer contributes only as an insulating layer, and the ice growth can be predicted using Eqn (21).

If the snow weight is larger than the buoyancy of the ice layer, then the ice surface can be flooded, and slush formation can occur. The initial slush thickness H _SL can be determined by a simple mass balance if two main assumptions are made: (1) routes (e.g. cracks) for the water are available for the immediate flooding of the snow/ice interface and (2) the effect of water capillarity on the snow is negligible (Leppäranta, Reference Leppäranta1983). The mass conservation states that the snow weight reduction is equal to ice buoyant weight increase, which for the initial slush thickness is given as:

(23)$$\rho _{\rm s}H_{\rm S}-\rho _{\rm s}H_{{\rm SL}} = {\rm \;}( {\rho_{\rm w}-\rho_{\rm i}} ) H_{\rm I} + ( {\rho_{\rm w}-\rho_{{\rm sl}}} ) H_{{\rm SL}}, \;$$

where ρ and H represent the snow, ice, slush and water, densities and thicknesses, respectively, and are indexed as s, i, sl and w.

If the buoyancy condition (Eqn 22) is not fulfilled and the slush layer is formed, its thickness can be determined from the mass conservation (Eqn 23). This is valid for the initial slush formation when the model consists of water, ice and dry snow which initiates the flooding. The initial slush layer thickness, which is derived from Eqn (23) is given as follows:

(24)$$H_{{\rm SL}} = {\rm \;}\displaystyle{{H_{\rm S}\rho _{\rm s}-H_{\rm I}( {\rho_{\rm w}-\rho_{\rm i}} ) } \over {\rho _{\rm s} + \rho _{\rm w}-\;\rho _{{\rm sl}}}}, \;$$

where H _S ρ _s is the weight of the snow layer and H _I(ρ _w − ρ _i) is the buoyancy B per unit area of the ice layer under the waterline.

If the slush thickness increases continuously due to snowfall before the slush-snow-ice transformation occurs the subsequent thickness increase of the slush layer is further estimated by adding in Eqn (23) the new buoyancy of layers under the waterline, which considers the buoyancy of ice and the initial slush thicknesses. Therefore, the following slush thickness increase is derived from the new mass conservation model that consists of water, ice, old slush layer and dry snow that continues to submerge. The balance between snow weight reduction and buoyancy increase is calculated as:

(25)$$\rho _{\rm s}H_{\rm S}-\rho _{\rm s}H_{{\rm SL\;}} = ( {\rho_{\rm w}-\rho_{\rm i}} ) H_{\rm I} + ( {\rho_{\rm w}-\rho_{{\rm sl}}} ) H_{{\rm SL\;old}} + ( {\rho_{\rm w}-\rho_{{\rm sl}}} ) H_{{\rm SL\;}}.$$

The slush thickness increase is as follows:

(26)$$H_{{\rm SL\;}} = \displaystyle{{H_{\rm S}\rho _{\rm s}-H_{\rm I}( {\rho_{\rm w}-\rho_{\rm i}} ) + H_{{\rm SL\;old}}( {\rho_{\rm w}-\rho_{{\rm sl}}} ) } \over {\rho _{\rm s} + \rho _{\rm w}-\;\rho _{{\rm sl}}}}, \;$$

where the term [H _I(ρ _w–ρ _i) + H _{SL old}(ρ _w–ρ _sl)] is the buoyancy B before the slush thickness increase occurs.

When the slush layer freezes completely, and the model takes into account the water, ice, snow ice and dry snow that initiates the slush formation then Eqn (25) is transformed to the following expression:

(27)$$\rho _{\rm s}H_{\rm S}-\rho _{\rm s}H_{{\rm SL}} = {\rm \;}( {\rho_{\rm w}-\rho_{\rm i}} ) H_{\rm I} + ( {\rho_{\rm w}-\rho_{{\rm si}}} ) H_{{\rm SI}} + ( {\rho_{\rm w}-\rho_{{\rm sl}}} ) H_{{\rm SL}\;}.$$

The slush thickness initially formed is then calculated as:

(28)$$H_{{\rm SL}} = {\rm \;}\displaystyle{{H_{\rm S}\rho _{\rm s}-H_{\rm I}( {\rho_{\rm w}-\rho_{\rm i}} ) + H_{{\rm SI}}( {\rho_{\rm w}-\rho_{{\rm si}}} ) } \over {\rho _{\rm s} + \rho _{\rm w}-\;\rho _{{\rm sl}}}}, \;$$

where ρ _si and H _SI are the density and thickness of the snow ice. In addition, the term H _I(ρ _w − ρ _i) + H _SI(ρ _w − ρ _si) is the buoyancy before the slush thickness continues to increase due to more snow is depressed from the incoming snowfall. The subsequent slush increase in a model that consists of water, ice, sea-ice, old slush and dry snow can be determined from the below buoyance balance:

(29)$$\rho _{\rm s}H_{\rm S}-\rho _{\rm s}H_{{\rm SL}} = {\rm \;}( {\rho_{\rm w}-\rho_{\rm i}} ) H_{\rm I} + ( {\rho_{\rm w}-\rho_{{\rm si}}} ) H_{{\rm SI}} + {\rm \;}( {\rho_{\rm w}-\rho_{{\rm sl}}} ) H_{{\rm SL\;old}} + ( {\rho_{\rm w}-\rho_{{\rm sl}}} ) H_{{\rm SL\;}}.$$

The slush thickness increase is as follows:

(30)$$H_{{\rm SL}} = \displaystyle{{H_{\rm S}\rho _{\rm s}-H_{\rm I}( {\rho_{\rm w}-\rho_{\rm i}} ) + H_{{\rm SI}}( {\rho_{\rm w}-\rho_{{\rm si}}} ) + H_{{\rm SL\;old}}( {\rho_{\rm w}-\rho_{{\rm sl}}} ) } \over {\rho _{\rm s} + \rho _{\rm w}-\;\rho _{{\rm sl}}}}, \;$$

where the term [H _I(ρ _w − ρ _i) + H _SI(ρ _w − ρ _si) + H _{SL old}] is B before the slush thickness addition. Based on Eqns (24, 26, 28, and 30), and assuming the slush thickness increase varies with the buoyancy before the slush addition, then a general formulation of the slush thickness is expressed as:

(31)$$H_{{\rm SL}} = \displaystyle{{W_{\rm s}-B} \over {\rho _{\rm s} + \rho _{\rm w}-\;\rho _{{\rm sl}}}}, \;$$

where W _s is the snow weight. Equation (31) has been previously used to estimate the slush thickness due to flooding by Saloranta (Reference Saloranta2000), Shirasawa and others (Reference Shirasawa2005), Cheng and others (Reference Cheng2014) and Wang and others (Reference Wang, Cheng, Wang, Gerland and Pavlova2015).

The ice, snow and water densities are usually assumed constants when estimating the ice thickness development. In addition, the difficulty in accurately measuring the slush density and the lack of data from field measurements can make the slush thickness estimation prone to error. Theoretical estimations of the maximum slush density were previously proposed by Adolphs (Reference Adolphs1998). The density was estimated based on the following assumptions: (1) the snow and the slush consist of two phases – ice/air and ice/water, (2) the snow air voids are replaced by water when the snow transforms to slush and (3) the air density is negligible in comparison with the water and ice densities. The snow consists of an ice fraction equal to V _i = ρ _s/ρ _i and an air fraction equal to V _a = (1 − ρ _s/ρ _i), where V _a is assumed to be equal to the volume fraction of water V _w. The slush mass is expressed as a sum of ice and water fraction masses and is taken as:

(32)$$\rho _{{\rm sl\;max}}V_{{\rm sl}} = {\rm \;}\rho _{\rm i}V_{\rm i} + \rho _{\rm w}V_{\rm w}.$$

Substituting the respective volume fractions yields the maximum slush density:

(33)$$\rho _{{\rm sl\;max}}{\rm \;} = \rho _{\rm w}-\rho _{\rm s}\left({\displaystyle{{\rho_{\rm w}} \over {\rho_{\rm i}}}-1} \right)\;.$$

As illustrated, if we consider ρ _s = 300 kg m⁻³, ρ _i = 917 kg m⁻³ and ρ _w = 1000 kg m⁻³, then the resulting ρ _{sl max} is 972.9 kg m⁻³.

The temperature is considered in the slush layer to equal the freezing point of the slush–water mixture, which varies with salinity. For air temperatures below freezing, the slush layer begins to freeze from above, at the snow/slush interface. The snow ice can also form in the slush/ice interface due to the stored heat capacity of ice. However, this phenomenon is assumed to have a low contribution to the snow-ice growth, and consequently, the ice layer is assumed to reach the freezing point instantly. During the snow-ice formation, the growth at the lower ice/water interface is assumed to be zero (Comfort and others, Reference Comfort, Gong, Singh and Abdelnour2003, Comfort and Abdelnour, Reference Comfort and Abdelnour2013). This implies that the temperature gradient between the upper and lower boundaries of ice is zero for the period that the ice has a temperature equal to T _F. The snow-ice formation from the freezing of slush can be determined from the continuity of the heat fluxes in the snow/slush interface, as illustrated in Table 1(e). The heat flux released from slush freezing is equal to the heat flux through the snow to the atmosphere, as follows:

(34)$$k_{\rm s}\displaystyle{{( {T_{\rm F}-T_{{\rm SA}}} ) } \over {H_{\rm S}}} = \left({1-\displaystyle{{\rho_{\rm s}} \over {\rho_{\rm i}}}} \right)\rho _{{\rm si}}L_{\rm i}\displaystyle{{{\rm d}H_{{\rm SI}}} \over {{\rm d}t}}, \;$$

where the term (1 − ρ _s/ρ _i) is the volume fraction of water V _w in the slush layer. Substitution of the surface temperature (Eqn 14), and integrating the above equation yields the snow-ice thickness, which is assumed to stop when all slush has refrozen:

(35)$$H_{{\rm SI}} = \displaystyle{{\theta \;} \over {( H_{\rm s}/k_{\rm s} + 1/h_{{\rm sa}}) V_{\rm w}\rho _{{\rm si}}L_{\rm i}}}.$$

Alternative estimates of the snow-ice formation include the so-called ‘snow to slush to snow-ice’ compression coefficient β (Leppäranta and Kosloff, Reference Leppäranta and Kosloff2000; Shirasawa and others, Reference Shirasawa2005; Ohata and others, Reference Ohata, Toyota and Shiraiwa2016). β is an empirical coefficient defined as the sum of the thickness change of snow when transformed to slush and when subsequently transformed to snow ice. The compression coefficient was estimated to be 1.5 from measurements conducted in Lake Pääjärvi where snow ice constituted 10–43% of the total thickness (Leppäranta and Kosloff, Reference Leppäranta and Kosloff2000). The value was estimated to be 2.0 for measurements conducted in Lake Abashiri, where the snow ice constituted 29–73% of the total thickness (Ohata and others, Reference Ohata, Toyota and Shiraiwa2016).

After the slush layer freezes, the subsequent growth of the ice is modelled considering the snow ice as a layer overlying the sea/ice, as illustrated in Table 1(e). From the continuity of heat fluxes at the interface of each layer: air/snow; snow/snow ice; snow ice/ice, and the integration of heat flux balance between the heat conducted through ice and the latent heat of freezing yields the following ice growth equations:

(36)$$\displaystyle{{{\rm d}H_{\rm I}} \over {{\rm d}t}} = {\rm \;}\displaystyle{1 \over {\rho _{\rm i}L_{\rm i}}}\left({\displaystyle{{T_{\rm F}-T_{\rm A}} \over {H_{\rm I}/k_{\rm i} + H_{{\rm SI}}/k_{{\rm si}} + H_{\rm s}/k_{\rm s} + 1/h_{{\rm sa}}}}} \right)\;, \;$$

(37)$$h_{{\rm ek}} = \displaystyle{{H_{\rm I}} \over {k_{\rm i}}} + \displaystyle{{H_{{\rm SI}}} \over {k_{{\rm si}}}} + \displaystyle{{H_{\rm s}} \over {k_{\rm s}}} + \displaystyle{1 \over {h_{{\rm sa}}}}, \;$$

(38)$$H_{\rm I} = \sqrt {k_{\rm i}^2 h_{{\rm ek}}^2 + \displaystyle{{2k_{\rm i}} \over {\rho _{\rm i}L_{\rm i}}}\theta } -k_{\rm i}h_{{\rm ek}}.$$

As an example of this method, the omission of the slush occurrence and snow-ice formation from the sea-ice thickness simulations for the monthly mean values of the winters of 1979–90, led to a total ice thickness difference of 25% (Saloranta, Reference Saloranta2000). There are, however, several sources of errors affecting the accuracy of the growth model. These include the assumption of the slush layer with an instant formation together with the slush to snow-ice transition, the thermal and physical properties of slush (Leppäranta, Reference Leppäranta1983; Saloranta, Reference Saloranta2000), as well as the availability of the cracks or other water routes through the ice layer (Ashton, Reference Ashton2011). These errors were noticed when thickness measurements were compared with the predicted values. The main model discrepancies were related to the time of the snow-ice formation in the individual annual predictions (Saloranta, Reference Saloranta2000), as illustrated in Figure 1. However, it is important to mention that the ice thickness presented by Saloranta (Reference Saloranta2000) was estimated based on a numerical approach, which computes the surface temperature from the surface heat budget (including the net radiation, sensible and latent heat fluxes).

Fig. 1. Observed (o) and modelled (m) snow (HS), snow ice (HSI) and total sea-ice thickness (HI) are illustrated with blue, red, black, dashed and solid lines, respectively. The measurements were conducted on fast sea ice in the Baltic Sea during the winter (1988 and 1989). Modified from Saloranta (Reference Saloranta2000).

Other approaches have focused on the effect of snow on lake ice development. A thermodynamic model that considered the thermal inertia through the ice was used to simulate the ice and snow thickness for the three consecutive winters (2009–12) (Cheng and others, Reference Cheng2014). The model solved the partial differential equations of conduction through the snow and ice layers and the surface heat fluxes were also estimated numerically. Compared to previous models i.e. models developed by Leppäranta (Reference Leppäranta1983) and Saloranta (Reference Saloranta2000), the snow and bottom ice thickness simulations showed higher accuracy. The freeboard simulations, however, had a higher deviation from the measurements, as shown in Figure 2. The freeboard error was judged to be a consequence of omitting the effect of water capillary on the snow layer.

Fig. 2. Modelled thickness of lake ice (HI (m)), snow (HS (m)) and freeboard (F (m)), for winter (2011–12) are illustrated with black, blue and red solid lines. The meteorological station data were used as input values for the thickness estimations. The measured thicknesses are shown with three different patterns. Modified from Cheng and others (Reference Cheng and Liang2014).

A partly numerical approach has been introduced recently aiming to study the effect of snowfall on the initial ice growth stage (Toyota and others, Reference Toyota, Ono, Tanikawa, Wongpan and Nomura2020). The initial ice growth was estimated by including the surface temperature in Stefan's analytical ice growth approach. The surface temperature was estimated numerically by considering the net radiation, sensible and latent heat fluxes. Observations showed that a 10 mm snowfall formed a slush layer that froze rapidly on a snow-ice layer above the initial 2 mm columnar ice layer. The study suggested that the snowfall contributed to the total ice thickness and affected the ice microstructure. However, freezing of a thin snow-slush layer at the initial ice formation did not have a significant effect on the subsequent ice growth rate.

2.5 Surface heat budget

The heat flux exchanged between the ice or snow surface with the atmosphere, denoted earlier as φ _IA or φ _SA, combines the effects of all heat fluxes that affect the surface temperature. An accurate forecast of the ice or snow surface temperature is essential for the ice growth at the lower surface of the ice, which is driven by the temperature gradient between boundaries. The surface temperature is a function of all the heat fluxes that flow in and out from the surface. The net radiation flux has a strong effect on this heat balance. Therefore, a numerical approach was developed to estimate the heat budget on the upper boundary of the system (Untersteiner, Reference Untersteiner1961, Reference Untersteiner1964), as in Table 1(f).

All the heat fluxes that flow towards the surface (down from the atmosphere) are considered positive and the heat fluxes flowing away from the surface (towards the atmosphere) are considered negative. The heat budget equation commonly used to estimate the surface temperature is:

(39)$$ \eqalign{& \varphi _{{\rm SE}\updownarrow } + \varphi _{{\rm L}\updownarrow } + \varphi _{\rm I} + ( {1-a} ) \varphi _{{\rm SW}\downarrow } + \varphi _{{\rm LW}\downarrow }-I_{0_z} \cr& \quad -\varphi _{{\rm LW}\uparrow } = 0, \;} $$

where φ _I, φ _SE, φ _L are the conducted, sensible and latent heat fluxes per unit area, respectively. In addition, φ _SW↓, φ _LW↓ and φ _LW↑ are the incoming shortwave radiation, incoming longwave radiation, and the outgoing longwave radiation heat fluxes, again per unit area.

The net longwave radiation influences the snow and ice-covered surface temperatures, particularly during the winter season when the incoming solar radiation is low. Stefan–Boltzmann's equation of the black-body radiation is widely used to parametrise the incoming and outgoing longwave radiation (φ _LW↓ and φ _LW↑), which is emitted from the atmosphere and snow or ice surfaces respectively. The outgoing longwave radiation is given by:

(40)$$\varphi _{{\rm LW}\uparrow } = \varepsilon _{\rm s}\sigma T_{\rm s}^4 , \;$$

where ɛ _s is the surface emissivity fraction, σ is the Stefan–Boltzmann constant (5.67 × 10⁻⁸ W m⁻² K⁻⁴) and T _s is the ice or snow surface temperature given in K. Earlier, this temperature was denoted for ice/atmosphere and snow/atmosphere interface as T _IA or T _SA, respectively. The incoming longwave radiation φ _LW↓ is linearly correlated to snow or ice surface emissivity that varies depending on the physical and optical properties of the surface (Rees and James, Reference Rees and James1992; Hori and others, Reference Hori2006; Cheng and Liang, Reference Cheng2014). The surface and sky emissivity can be simply estimated from the ratio of the measured longwave radiation with black-body radiation of the same temperature (Konzelmann and others, Reference Konzelmann1994).

The sky emissivity depends on vertical temperature profile, humidity, vapour pressure, precipitated water and cloud occurrence (Niemelä and others, Reference Niemelä, Räisänen and Savijärvi2001). Equations of the downwelling radiation based on different sky conditions are summarised in Table 2.

Table 2. Parametrisation of downwelling longwave radiation for clear and overcast skies

T _A is the air temperature (K); e ₀ is the vapour pressure (hPa); σ is the Stefan–Boltzmann constant equal to 5.67 × 10⁻⁸ W m⁻² K⁻⁴; c is the cloud fraction in tenths; w _p is the precipitated water (cm) and the notation cl stands for a clear sky.

The incoming shortwave radiation in spring and summer causes ablation and can retard the freezing process in the autumn and early winter (Perovich and Grenfell, Reference Perovich and Grenfell1981). The shortwave radiation that is absorbed on the surface I _s and influences the surface temperature (Eqn 39) is expressed as:

(41)$$I_{\rm s} = {\rm \;}( {1-a} ) \phi _{{\rm SW}\downarrow }-I_{0_z}, \;$$

where φ _SW↓ is the incoming shortwave radiation and a is the albedo which is defined as the part of the incident irradiance that is reflected in the atmosphere. It is considered as the ratio of the upwelling irradiance to the down-welling irradiance (Perovich, Reference Perovich1990). Albedo values change from high values with low-spatial variation during winter due to snow-covered surfaces to lower values with higher-spatial variation during spring and summer due to possible melting ponds, bare ice and open water (Perovich and others, Reference Perovich, Roesler and Pegau1998, Reference Perovich, Grenfell, Light and Hobbs2002a, Reference Perovich, Tucker and Ligett2002b). Additionally, a vertical variation of the total albedo is possible due to changes in sea-ice thickness, microstructure and brine content. Then, for an ice cover thicker than 0.8 m, the albedo approaches a constant value (Perovich, Reference Perovich1996). The first term in Eqn (41) represents the amount of incoming shortwave radiation that is not reflected. The radiation component $I_{0_z}$ is the fraction of shortwave radiation that penetrates deeper into the snow or ice layers, which can also be defined as transmitted radiation (Grenfell and others, Reference Grenfell, Light and Perovich2006). This can be thought of as the ratio between the measured light intensity in certain depths and the measured incoming radiation. In addition, the transmitted radiation $I_{0_z}$ (intensity measured) at each depth and the absorbed shortwave radiation I_z are related by the Bouguer–Lambert equation. The solar radiation absorption is assumed to decrease exponentially with depth (from the surface coordinate z) as follows:

(42)$$I_{0_{z}} = I_z \times {\rm e}^{{-}k_{{\rm abs}}z}, \;$$

where k _abs is the radiation absorption coefficient. The absorbed radiation is sensitive to ice microstructure (Perovich, Reference Perovich1996) and the bulk ice density (Jin and others, Reference Jin, Stamnes, Weeks and Tsay1994). The air content scatters the light, and the brine inclusions both scatter and absorb the light. Therefore, higher ice bulk densities absorb more radiation.

In addition to the effect of net radiation on the surface heat budget (Eqn 39), the latent heat φ _L that is the energy gained or released during surface evaporation and sublimation, and the sensible heat flux that is the heat exchanged between water or ice surface with the atmosphere φ _SE, were found to be an important heat source, especially in the initial phase of the ice growth (Toyota and Wakatsuchi, Reference Toyota and Wakatsuchi2001; Toyota and others, Reference Toyota, Ono, Tanikawa, Wongpan and Nomura2020). The sensible and latent heat fluxes can be determined from equations proposed by Maykut (Reference Maykut1978):

(43)$$\varphi _{{\rm SE}} = \rho _{\rm a}c_{\rm a}k_{{\rm se}}w( {T_{\rm A}-T_{\rm S}} ) , \;$$

(44)$$\varphi _{\rm L} = L_{{\rm sb}}k_{\rm e}w( {q_{\rm a}-q_{\rm s}} ) .$$

The sensible heat flux φ_SE represents the convective heat transfer between i.e. air and snow. Its value depends on the temperature differences between air and snow or ice surface, the air density ρ _a, the air specific heat capacity c _a, the sensible heat transfer coefficient k _se and the wind speed w. Based on Eqn (43), a wind speed equal to zero yields a sensible heat flux equal to zero, which implies that there is no heat exchange between the ice surface and atmosphere. The sensible heat flux was estimated, for the initial ice growth from 0.05 to 0.2 m, to vary from 29 to 362 W m⁻² (Maykut, Reference Maykut1978).

Latent heat flux represents the heat released due to the ice formation on the surface by sublimation. L _sb is the latent heat of sublimation; k _e is the evaporation coefficient; q _a and q _s are the specific humidity in the air and the surface; the sensible and latent heat fluxes are substituted in the surface heat budget (Eqn 39). The latent heat flux was estimated to vary from 12 to 36 W m⁻² for the initial ice growth from 0.05 to 0.2 m (Maykut, Reference Maykut1978). Consequently, the heat fluxes suggested by Maykut (Reference Maykut1978) for an ice growth from 0 to 0.2 m are considerable. In addition, Toyota and others (Reference Toyota, Ono, Tanikawa, Wongpan and Nomura2020) showed that for the first 10 mm of the initial ice growth from open water and the initial consolidation of 10 mm slush layer, the sensible heat fluxes were ~27 and 26 W m⁻², respectively. For the same conditions, the latent heat fluxes were 17 and 16 W m⁻². Also, the change in ice thickness increase from 10 to 30 mm decreased the sum of sensible and latent heat fluxes by 3 and 7 W m⁻², respectively. The contribution of the sensible and latent heat fluxes to the initial phase of ice growth may depend on the initial energy state of the open water, i.e. the oceanic flux (which is discussed later), and the energy stored in the water during summer radiation. For example, during brash ice growth, these heat fluxes contribute to the growth of ice in the open water fraction. The influence of sensible and latent heat fluxes on ice thickness is not well defined and should be further investigated.

2.6 Oceanic heat flux

The oceanic heat flux is the heat flux transferred from water to ice at the ice/water interface. It is generated from the temperature gradient developed in the oceanic vertical water column. Three different processes can be attributed to affecting the temperature gradient in the water column: turbulent currents induced mainly by the wind stress; convection energy exchange in the interface of the oceanic water layers which are stratified based on density difference, and the tidal currents (Gabison, Reference Gabison1987).

The oceanic heat flux is generally included in the seasonal or annual ice thickness prediction models as an external (constant) heat source. It is estimated from the observed heat flux balance at the ice–water interface and calculated from the difference between the total energy conducted through the ice to the atmosphere and the ice freezing latent heat flux (Shirasawa and others, Reference Shirasawa, Ingram and Hudier1997; Ohata and others, Reference Ohata, Toyota and Shiraiwa2016):

(45)$$\varphi _{\rm W} = k_{\rm i}\displaystyle{{\Delta T} \over {\Delta z}}-\rho _{\rm i}L_{\rm i}\displaystyle{{{\rm d}h} \over {{\rm d}t}}.$$

Another formulation of the ice/water interface energy budget equation commonly used to estimate the oceanic heat flux is as follows (McPhee and Untersteiner, Reference McPhee and Untersteiner1982):

(46)$$\varphi _{\rm i} + \varphi _{\rm F} + \varphi _{\rm S}-\varphi _{\rm W} = 0, \;$$

where φ _i, φ _F, φ _w and φ _s are the conductive, freezing, oceanic and thermal inertia heat fluxes.

The oceanic heat flux numerical value that is generated from the temperature gradient in the water column can be estimated after an equation proposed by Ohata and others (Reference Ohata, Toyota and Shiraiwa2016):

(47)$$\varphi _{\rm W} = \rho _{\rm w}c_{\rm w}k_{\rm w}W( {T_{\rm w}-T_{\rm F}} ) .$$

where ρ, c and k are the density, specific heat capacity and conductivity of water. W is the water current velocity and (T _w − T _F) is the temperature gradient between the ice/water interface and water under the ice. For the lake ice growth (Ohata and others, Reference Ohata, Toyota and Shiraiwa2016), the oceanic heat flux varied between 1 and 11 W m⁻². Equation (47) shows that if there are no water currents, the oceanic heat flux is zero even when there is a temperature gradient present. In addition, Eqn (47) neglects the depth of the water column by considering the water temperature under the ice, which can reach the freezing point quickly, and consequently assumes that the oceanic heat flux is an important heat source only in the initial phase of ice growth. This can be considered a reasonable approximation when modelling the landfast ice and lake ice growth.

Solar radiation was found to have a significant contribution to the magnitude of the oceanic heat flux measured under a drifting ice pack in the central Arctic (Maykut and McPhee, Reference Maykut and McPhee1995). The seasonal variation of the oceanic heat flux was suggested to be important when estimating the annual thermodynamic ice growth (Wang and others, Reference Wang, Cheng, Wang, Gerland and Pavlova2015). Thickness estimations of Arctic fast ice for three different oceanic heat flux inputs (φ _w = 3, 5 and 10 W m⁻²) resulted in high thickness variations. While recent observations of the turbulent currents under drifting ice located in the Arctic showed that storms could double the oceanic heat flux in the ice/water interface (Peterson and others, Reference Peterson, Fer, McPhee and Randelhoff2017), as shown in Figure 3. In addition, oceanic heat fluxes were found to be as high as 10–15 W m⁻² for the initial ice growth period in the central Arctic (Lei and others, Reference Lei2014). Despite the seasonal changes in the oceanic heat flux, the warm Atlantic streams and the effect of wind, the oceanic heat flux during the winter growing season in the central Arctic for three consecutive winters (2010–14) varied from 2 to 4 W m⁻² (Lei and others, Reference Lei2018). Therefore, the importance of the oceanic flux in the ice growth estimations can vary on the local atmospheric conditions including winds, water currents as well as the river's turbulent current (as commonly observed in the Baltic sea), and the solar radiation that over summer contributes to the water's temperature with interannual variations.

Fig. 3. Relative frequency of the oceanic heat flux φ _w (W m⁻²), during calm and windy conditions respectively the shade bars and the solid lines for two different cases, the Arctic winter (for water depth above 3750 m) and over the Atlantic water (for water depth below 2000 m). Reprinted from Peterson and others (Reference Peterson, Fer, McPhee and Randelhoff2017).

3.1 Brash ice accumulation

Brash ice growth models draw many elements from the level ice models, such as the description of the consolidation process. However, the brash ice models are more complex due to mechanical and hydrodynamical effects. Brash ice models have been developed based on Stefan's growth approach to predict the ice volume accumulated in the ship tracks after each breaking event. The first brash ice model was reported by Ashton (Reference Ashton1974). After each passage, part of the channel's surface is assumed to consist of open water. The channel's surface fraction occupied by open water is denoted as γ. The surface fraction of the channel surface occupied by the brash ice pieces is 1 − γ. The open water fraction was assumed constant for the performed model validation (20–30%). However, in reality, the surface open water fraction in the initial phase of brash ice formation differs from the second accumulation phase where γ reaches a constant value (Ashton, Reference Ashton1974). In the initial phase of brash ice formation, the open water fraction reduces until the brash ice channel surface is covered with ice and the second brash ice development phase begins.

The consolidation between two vessel passages is assumed to occur in the open water fraction in the same manner as in the open sea. The brash ice pieces continue to consolidate from the fraction of the ice pieces that are in the channel after the breaking event. It is assumed that the amount of ice that remains in the channel is equal to the amount of ice before the breaking event. In Ashton's model, the brash ice macroporosity is neglected, and as such, the brash ice average thickness (ice fraction + water fraction) formed after the ship passage is neglected. Therefore, the brash ice fraction in the channel after the ship passage is considered to be the former brash ice consolidated thickness H _C,(i−1), where the index i refers to the number of ship passages. Under these assumptions and neglecting ice volume loss caused by the sideways movement, the consolidation of the brash ice between two vessel passages H _C,i was formulated as follows:

(48)$$H_{{\rm C, \;}i} = [ {H_{{\rm C}, ( {i-1} ) }^2 + {( {1-\gamma } ) }^2\alpha^2\Delta \theta_i} ] ^{0.5} + \gamma \alpha \Delta \theta _i^{0.5} .$$

The first term in Eqn (48) represents the growth of the brash ice fraction and the second term represents the ice growth in the open water surface. Δθ_i = θ_i − θ_{i −1} is the cumulative degree-days between two vessel passages i − 1 and i. The ice growth rate coefficient used for the brash ice estimation was 18 mm °F^−0.5 d^−0.5 in the example case for the Mississippi River.

Another equation was proposed to predict the brash ice growth between two consecutive ship passages by Sandkvist (Reference Sandkvist1980, Reference Sandkvist1981). The initial ship transit was assumed to occur through an intact ice layer with an initial thickness of H ₀. Ship passages were assumed to generate an equivalent brash ice thickness in the channel. The equivalent brash ice layer is the uniform brash ice thickness where all the ice pieces generated after each breaking event are considered to be in a channel that has a width equal to the vessel's beam. In addition to the previous assumptions, the open water fraction is neglected (this means, in effect, that the initial formation phase is neglected) and the brash ice consolidation between two consecutive passages H _C,i is formulated as the following:

(49)$$H_{{\rm C, \;}i} = H_0 + \mathop \sum \limits_{i = 1}^i \alpha ( {\Delta \theta_i} ) ^{0.5}.$$

The ice consolidation rate α proposed for estimating the brash ice accumulation in the fast sea-ice channels in the Bay of Bothnia in the Baltic was based on empirical data fitting and resulted in a value equal to 12 mm °C^−0.5 d^−0.5. The Sankvist formulation was used to estimate the equivalent brash ice thickness generated in the Saimaa Canal, in the Gulf of Finland, during the winters of 1980–82 (Eranti and others, Reference Eranti, Penttinen and Rekonen1983). The consolidation rate coefficient was found to have a good fit with the measured thicknesses when α = 6.5 mm °C^−0.5 d^−0.5.

In addition to the assumptions made for Stefan's ice growth, the above-mentioned models neglect the ice growth caused by the thermal inertia of ice pieces submerged after each breaking event and also any explicit inclusion of porosity. Furthermore, the assumption that all the ice pieces remain in the channel after every ship transit results in an overestimation of the brash ice thickness (Ettema and Huang, Reference Ettema and Huang1990; Riska and others, Reference Riska2019).

Another model was developed to estimate the brash ice consolidation between two vessel passages H _C,i, and the brash ice thickness H _Bi that remains in the channel after each breaking event (Ettema and Huang, Reference Ettema and Huang1990). The model assumed that the brash ice layer after each vessel passage consists of a uniform layer of the brash ice and water mixture, which has a constant porosity p before and after each passage. Also, a constant open water fraction was assumed. In addition to previous models, the ice expulsion sideways was considered. The thermal growth equation was developed based on Stefan's theoretical assumptions discussed previously. The ice consolidation between two vessel passages was calculated as a sum of the open water consolidation component and the thermal growth of an existing brash ice thickness. This thickness was equal to the total thickness of brash ice multiplied by its ice fraction content (1 − p)H _Bi,(i−1). The brash ice consolidation between two vessel passages was formulated as:

(50)$$\eqalign{ & H_{{\rm C}, \;i} = \;H_{{\rm C}, ( {i-1} ) } + \gamma \alpha \left({\mathop \sum \limits_{i = 1}^i \Delta T( P_i-t_{\rm p}) } \right)^{0.5} \cr & \!\! + \!( {1\! - \!\gamma } ) \left[{{\left({( {1 \!- \!p} ) H_{{\rm Bi}, ( {i-1} ) }a^2\mathop \sum \limits_{i = 1}^i \Delta T( P_i-t_{\rm p}) } \right)}^{0.5} \!\!\! -( {1\!-\!p} ) H_{{\rm Bi}, ( {i-1} ) }} \!\right], \;} $$

where H _C,(i−1) is the former thickness of brash ice in the channel. The second term estimates the ice growth in the open water fraction and is modified to estimate the cumulative freezing degree-days between two vessel passages. This is the difference between the two vessel transit times P_i and the time needed for the ice to reach the thermal equilibria t _p before the growth begins. However, the time needed to reach the thermal equilibrium varies with the former brash ice transit period and the cumulative freezing degree-days within that period. Therefore, its determination is not straightforward. The third term is the ice growth of a brash ice volume which has an initial thickness of (1 − p)H _Bi,(i−1).

The brash ice consolidation can, however, be simplified and corrected by assuming that the volume increase is a sum of the growth in the open water fraction γ and the consolidation in the macropores of brash ice, which is described by the porosity p. Based on Stefan's thermal ice growth in open water and the macropores, then the ice consolidation between two vessel passages can be estimated as:

(51)$$ \eqalign{ H_{{\rm C}, \;i} = H_{{\rm C}, ( {i-1} ) } + \alpha \left({\displaystyle{1 \over \gamma }\mathop \sum \limits_{i = 1}^i \Delta T( P_i-t_{\rm p}) } \right)^{0.5} \cr + \alpha \left({\displaystyle{1 \over {\,pH_{{\rm Bi}, ( {i-1} ) }}}\mathop \sum \limits_{i = 1}^i \Delta T( P_i-t_{\rm p}) } \right)^{0.5}, \; } $$

where the first term is the ice volume before the breaking event, the second term is the ice growth in the open water and the third term is the growth of ice within pores of the new brash ice having a uniform thickness.

The brash ice uniform thickness H _Bi after each ship passage is derived from the ice mass conservation equations before and after the shipping transit (Ettema and Huang, Reference Ettema and Huang1990). Thus, it is assumed that the brash ice fraction after each breaking event (1 − p)H _Bi is equal to the sum of the former brash ice fraction before consolidation (1 − γ)(1 − p)H _Bi,(i−1), with the ice volume increase between two passages (H _C,i − H _C,(i−1)). Thus, the brash ice thickness is given as:

(52)$$H_{{\rm Bi}} = {\rm \;}\displaystyle{{( {1-\varepsilon } ) [ {H_{{\rm Bi}, ( i-1}) ( {1-p} ) ( {1-\gamma } ) + ( {H_{{\rm C}, i}-H_{{\rm C}, ( {i-1} ) }} ) } ] } \over {( {1-p} ) }}\;, \;$$

where ɛ represents the fraction of the brash ice expelled sideways and is assumed to be constant.

Another brash ice accumulation model was developed by Riska and others (Reference Riska2019). The brash ice was assumed to consist of three layers; the dry brash ice above the water level H _d, the consolidated brash ice thickness H _C,i and the brash ice thickness submerged in the water below the consolidated layer called wet brash ice H _b. Before a breaking event, the fraction of dry brash ice above water level and the fraction of the ice below were estimated based on the buoyancy principle as follows:

(53)$$H_{{\rm d}, ( {i-1} ) \;} = H_{{\rm B}, ( {i-1} ) }\left({1-\displaystyle{{\rho_{\rm i}} \over {\rho_{\rm w}}}} \right), \;$$

(54)$$H_{{\rm b}, ( {i-1} ) } = H_{{\rm B}, ( {i-1} ) }-\left({H_{{\rm B}, ( {i-1} ) }\left({1-\displaystyle{{\rho_{\rm i}} \over {\rho_{\rm w}}}} \right)-H_{{\rm C}, ( {i-1} ) }} \right), \;$$

where H _B,(i−1) is the total brash ice equivalent thickness before the breaking event, and H _C,(i−1) is the consolidated layer of brash ice formed between consecutive breaking events and is determined from the last step of the breaking cycle. For simplicity, the model was divided into phases as shown in Figure 5. As the first step, the model considered the breaking event, which forms a new brash ice layer with a thickness H _B and porosity p ₀, and a uniform temperature T _eq (equilibrium temperature). In this step, the model estimates the equilibrium temperature and the new equivalent brash ice thickness H _B. T _eq is obtained from the mass and heat conservation before and after the breaking event given as:

(55)$$ \eqalign{T_{{\rm eq}} = & \,\,{{( {( T_{{\rm DI}} + T_{{\rm DA}}) /2} ) ( {1-p_0} ) H_{{\rm d}, ( {i-1} ) } }}\cr& {{ + ( {( T_{\rm F} +1T_{{\rm DI}} ) /2} ) H_{{\rm I}, ( {i-1} ) } + T_{{\rm F\;}} ( {1-p_{i-1}} ) H_{{\rm b}-1}} \over {H_{\rm B} ( {1-p_0} ) }}, \; }$$

where the interface temperatures T _DA and T _DI are determined from the continuity of the heat flux balance at the interfaces as discussed earlier. p_{i −1} and p ₀ are the wet brash ice porosity before the breaking event and the dry brash ice porosity before the breaking event. In addition, the equivalent brash ice thickness after the breaking event is obtained from the mass conservation of the ice before and after the ship passage:

(56)$$H_{\rm B} = H_{{\rm d}, ( {i-1} ) } + \displaystyle{{H_{{\rm I}, ( {i-1} ) }} \over {( {1-p_0} ) }} + H_{{\rm b}, ( {i-1} ) }\displaystyle{{( {1-p_{i-1}} ) } \over {( {1-p_0} ) }}.$$

Fig. 5. Scheme of brash ice models where the brash ice thickness before breaking event is divided into three layers: the dry layer above the water level with air-filled voids, the consolidated ice layer and the water-filled pores brash below solid ice. After the ice-breaking incident, the brash ice profile was assumed to have a vertically uniform temperature and porosity. Ice temperature and porosity are redistributed, and refreezing occurred before another breaking event. Adapted from Riska and others (Reference Riska2019).

As the second phase of the brash ice accumulation, the model considers the wet brash ice porosity change from p ₀ to p _i right after each breaking event. The porosity change is caused by the initial thermal inertia brash ice growth and is assumed to occur instantly after the breaking event. However, the porosity of the dry brash ice after the breaking event remains the same (p ₀). The new brash ice porosity is calculated from the heat balance between ice pieces and the surrounding water. The heat released due to the freezing of water is equal to the heat absorbed by the ice pieces that manifest in a temperature increase. Then the new brash ice porosity p _i is calculated as:

(57)$$p_{\rm i} = p_{0\;}\left[{\displaystyle{{c_{\rm i}\Delta T} \over {L_{\rm i}}} + 1} \right]-\displaystyle{{c_{\rm i}\Delta T} \over L}, \;$$

where ΔT = (T _F − T _eq) is the average temperature difference between the ice pieces and surrounding water. In the last step, the model estimates the thermal consolidation of the brash ice layer based on Stefan's ice growth model. The schematic illustration of this step is given in Figure 6. The thickness equation is determined under the same approach and derivations as the atmosphere/dry snow/ice growth model discussed earlier in Section 2.3. The two main changes compared to Ashton (Reference Ashton1974) and Sandkvist (Reference Sandkvist1980, Reference Sandkvist1981) incorporated in the model are the dry brash layer with its insulating effect, and the porosity p, which is added to the freezing heat flux equation. It is assumed that freezing occurs in macropores filled with water present in the brash ice layer. The model scheme and heat flux continuity equations at each interface are illustrated in Figure 6. The substitution of the interface temperatures at heat flux continuity equations and the integration of the heat flux balance at the lower surface of the consolidated layer (Eqn 3 in Fig. 6) yields the brash ice consolidation thickness equation:

(58)$$\displaystyle{{T_{\rm F}-T_{\rm A}} \over {H_{\rm I}/k_{\rm i} + h_{\rm e}}} = p_{\rm i}\rho _{\rm i}L_{\rm i}\displaystyle{{{\rm d}H_{\rm I}} \over {{\rm d}t}}, \;$$

which can be integrated to give the consolidated layer thickness as:

(59)$$H_{{\rm IC}} = \sqrt {k_{\rm i}^2 h_{\rm e}^2 + \displaystyle{{2k_{\rm i}} \over {\,p_{\rm i}\rho _{\rm i}L_{\rm i}}}\theta } -k_{\rm i}h_{\rm e}, \;$$

(60)$$h_{\rm e} = \displaystyle{{H_{\rm d}} \over {k_{\rm d}}} + \displaystyle{1 \over {h_{{\rm da}}}}, \;$$

where h _e is the equivalent heat transfer coefficient and k _d is the conductivity coefficient of the dry brash ice layer, the value of which depends on the porosity and was assumed equal to 1.31 W m⁻¹ K⁻¹ for p ₀ equal to 0.2.

Fig. 6. Atmosphere–surface coupling growth model of the consolidated brash ice between two consecutive ship passages. H _d, H _C,i and H _b represent, respectively, the thicknesses of the dry brash ice, consolidated brash ice and wet brash layer. Adapted from Riska and others (Reference Riska2019).

3.2 Porosity and piece size distribution

In partially consolidated brash ice or ice ridges, the macroporosity is defined as the ratio between the volume of water and air with the total volume bordered by the outer surfaces of brash ice. Porosity is an important parameter that influences the consolidation of the brash ice and side ridges. The porosity is mainly influenced by the piece size distribution and shape of the brash ice pieces. These are affected by particulars of the ship, breaking frequency, former brash ice consolidation thickness, the time that allows the ice piece to settle t _p (initial consolidation) and the confinement by side ridges. In addition, the packing of the brash ice pieces is also affected by the side ridges' thickness, which after several passages forces the ice to remain in the channel. The porosity also reduces consolidation during t _p. The period between breaks is important when the brash ice channels are navigated very frequently, such that it reduces or prevents the consolidation between passages.

An early study of the full-scale brash ice channel in the Bay of Bothnia showed that the brash ice water content increases with depth from 20% of water content in the surface to 100% water below the brash ice layer (Sandkvist, Reference Sandkvist1980). The water content was determined across the channel's cross section at 0.05 m depth intervals. The water content at each depth interval (i.e. from 0 to 0.05 m; 0.05 to 0.1 m; 0.1 to 0.15 m, etc.) was calculated as a ratio between the number of drilled holes that contained water with the total amount of holes observed across the channels' cross section, i.e. assuming a uniform brash ice layer. The number of measurements across the cross section in each 50 mm depth varied from 30 to 50 observations. In other words, in each drilled hole, the ice or water occurrence is recorded for the 0.05 m depth intervals. For example, in the depth interval from 0.2 to 0.25 m, ten holes (from 30 holes in total) across the channel cross section contain water, and 20 holes contain ice. Subsequently, the water content in the depth interval 0.2–0.25 m is 10/30 which is equal to 33% of water content. The estimated change of water content with depth is shown in Figure 7. Results showed an almost linear increase of water content with time and breaking events. It should, however, be recognised that this is not a true indication of the channel porosity (and thus referred to as water content) since it assumes a uniform brash layer when in reality it is varied across the channel.

Fig. 7. Development of water content, assuming a uniform brash ice layer across the channel, with time and cumulative freezing degree-days. Where 13/12/78-14.5 shows the time of measurement and the cumulative temperature from the last passage. The linear correlation coefficient of porosity as a linear function of depth is given. Observations were conducted in a full-scale brash ice channel in Luleå harbour, Sweden during winter 1978–79 by Sandkvist (Reference Sandkvist1980). Adapted from Sandkvist (Reference Sandkvist1980, Reference Sandkvist1986).

Another study that focused on the ship resistance in brash ice reported a model scale brash ice channel porosity equal to 46% (Kitazawa and Ettema, Reference Kitazawa and Ettema1985). A more recent study of a model scale brash ice channel observed porosity values between 40 and 60% (Bridges and others, Reference Bridges, Riska, Suominen and Haase2020). In addition, the pressure ridge porosity can be a useful reference for brash ice studies. A study that focused on the first-year pressure ridges has reported rubble ice macroporosity that varied between 10 and 27% (Ervik and others, Reference Ervik, Høyland, Shestov and Nord2018). Higher values reported by Høyland (Reference Høyland2007) reached 45%. However, the rubble ice forms naturally from compressing or shearing of level ice (Høyland, Reference Høyland2007), while brash ice continuously changes in shape, size and layer packing due to the freeze-breaking cycle. Brash ice pieces are usually spherical ice blocks and the rubble ice are naturally sharped edge blocks that change in size mainly from consolidation and melting processes. Although the rubble ice macroporosity in ice ridges may differ from the brash ice macroporosity, it can be a useful reference for the brash ice layer and brash ice side ridge macroporosity, considering that limited brash ice channel measurement data are published on this aspect. As such, it is difficult to compare the full-scale brash ice porosity with model-scale brash ice porosity, especially when considering all the factors that affect this parameter. For example, the piece sizes and shapes of the model brash ice are different from the full-scale brash ice, consequently, the porosity may be different, see Matala (Reference Matala2021). This comparison requires porosity measurements in full-scale brash ice channels and brash ice accumulated in port areas.

The physical properties of three different brash ice models were investigated and reported by Matala and Skogström (Reference Matala and Skogström2019) and Matala (Reference Matala2021). The measured brash ice porosity varied between 30 and 55.1%. In addition, the research reported that the piece size distribution of the model brash ice channel was limited by the thickness of the original ice layer. Measurements on a full-scale brash ice channel in the Barents Sea reported a lognormal piece size distribution with sizes varying between 0.06 and 0.7 m (Bonath and others, Reference Bonath, Zhaka and Sand2019), this distribution trend is also supported by Tuovinen (Reference Tuovinen1979), although the piece size measurements were limited to the equipment size that was used to collect the brash ice. The vertical piece size distribution measured in a full-scale testing channel in the Bay of Bothnia reported a mean piece size of 0.7 m (Sandkvist, Reference Sandkvist1982).

Further investigation and research of porosity as a function of piece size, former thickness and ship passages could give a better understanding of the brash ice development and increase the brash ice forecast accuracy.

3.3 Consolidation

Between two consecutive vessel passages, the brash ice consolidates and its bulk temperature cools according to the cumulative freezing temperature. The consolidation process is defined as the downward freezing of ice initially caused by the stored specific heat of ice (initial phase of consolidation) and subsequently by the increase of the cumulative freezing degree-days (Høyland, Reference Høyland2002). Heat exchange between water and cold submerged ice blocks reach a temperature equilibrium, which for water and ice mixture is the freezing temperature that results in ice formation and is defined as the initial phase of consolidation. Consequently, after each vessel passage, the average macroporosity in the brash ice layer continuously decreases with time until the next breaking event.

The consolidation process of side ridges piled on both sides of a ship channel may be similar to the first-year pressure ridge consolidation. The theories and findings of ridge consolidation can be applied in the brash ice side ridge consolidation estimations. However, the formation phase and its duration make an important difference between these two features. The main challenge is estimating the brash ice side ridge accumulation, which is a gradual process, as the brash ice side structures grow in thickness following each ship transit. In addition to ice side way movement, the turbulent currents caused by ship passages could be a heat source that can affect the side ridge consolidation, however these aspects have not been studied.

Similar to ice ridges, consolidation of rubble ice provides a case that can be considered similar to that of brash ice. Studies on the initial phase of rubble ice consolidation have shown that this phase begins with the submergence of ice and ends when a thermodynamic temperature equilibrium between water and ice is reached (Høyland and Liferov, Reference Høyland and Liferov2005). The ice growth and heat exchange between ice and water during the initial consolidation phase were found to depend on the initial rubble temperature, size of the rubble, as well as on the salinity (Høyland and others, Reference Høyland2001; Høyland and Liferov, Reference Høyland and Liferov2005; Chen and Høyland, Reference Chen and Høyland2016; Chen and others, Reference Chen, Høyland and Ji2020).

The influence of the rubble thickness in the initial ice consolidation was investigated by Høyland and others (Reference Høyland2001). Three ice layers were prepared from brackish water having a salinity of 8.5 ppt. The ice layers had an initial temperature of −32°C and thicknesses of 24.5, 48.5 and 96 mm. The ice layers were insulated from the top and the sides thus hindering the effect of air temperatures on the ice growth. The growth in the ice/water interface caused by the thermal gradient between the initial ice temperature and the water temperature was measured. The ice layers grew 61, 41 and 26% for the three different thicknesses, during the first 250 min. This study highlights the influence that the rubble thickness can have on the initial phase of consolidation. Similar to level ice growth, during the initial consolidation the heat conducts faster through the thinner rubble ice and thus resulting in quicker consolidation.

Another laboratory experiment was conducted to compare the effect of salinity on the initial stage of consolidation by Chen and Høyland (Reference Chen and Høyland2016). Two ice samples were prepared with the same initial temperature of −32°C and thickness of 100 mm. The first sample was prepared from fresh water and the second sample had a salinity of 2.65 ppt. The ice blocks were submerged in fresh water with an initial temperature of 0.2°C. The top and side surfaces of the ice rubble were insulated and the growth in the ice/water interface was measured. The freshwater ice grew 15.9% and the saline ice grew 26%. Even though the growth of the saline ice sample was higher, the thermal equilibrium for the freshwater ice reached faster, 500 and 2500 min, respectively. This study highlighted that the saline ice has higher specific heat stored in ice and therefore the initial consolidation phase for saline ice lasts longer and the ice growth is higher.

The same laboratory tests were conducted for five freshwater ice samples and the results have been reported recently by Chen and others (Reference Chen, Høyland and Ji2020). Four samples had an initial temperature of −20°C and an initial thickness of 49, 64, 102 and 205 mm. A fifth sample had an initial ice temperature and thickness of −32°C and 101 mm. The experiments showed that the ice growth was greater for the coldest sample. On the contrary, the thicker ice samples required a longer time to reach the temperature equilibrium between the ice rubble and underlying water. This study highlighted that the initial ice temperature and thickness influence the initial phase of consolidation.

In addition to the effects of ice thickness, temperature and salinity on the initial rubble consolidation, the heat flux exchange between water and ice in the ice/water interface for the insulated systems discussed above was estimated by Chen and Høyland (Reference Chen and Høyland2016) and Chen and others (Reference Chen, Høyland and Ji2020). In this closed system, the temperature gradient between submerged ice rubble and freezing water drives the ice formation in the ice/water surface. The amount of ice formed or the amount of the latent heat released during the ice formation is equal to the difference between the energy stored in the ice rubble (the specific heat capacity of ice which for freshwater ice depends only on the initial temperature) and the heat convected by the water. The water heat flux was determined to be the difference between the specific heat capacity flux of the rubble ice with the latent heat flux. This stored energy (specific heat), in each time interval of ice formation defines the potential for further ice growth. When the thermal equilibrium (between ice and water) is reached the heat capacity of ice theoretically reaches zero. At the beginning of the experiment, the maximum heat fluxes were 280 and 380 W m⁻² for the saline and freshwater ice samples, respectively (Chen and Høyland, Reference Chen and Høyland2016). The value of heat flux for the freshwater ice approached a constant value with time. However, the saline water reached a minimum after 2 h of exposure and increased again, postulated to be due to brine expulsion, as shown in Figure 8.

Fig. 8. Water heat transferred and ice growth evolution with time. The ice samples had different salinities but the same initial temperature (−35°C) and were submerged in fresh water with an initial temperature equal to 0.2°C. SI and FWI are, respectively, saline and freshwater ice. Adapted from Chen and Høyland (Reference Chen and Høyland2016).

Another study estimated the convective heat flux in the initial formation phase of two pressure ridges artificially formed in Svalbard, Norway (Høyland and Liferov, Reference Høyland and Liferov2005). A simple numerical heat flux conservation equation was used:

(61)$$\varphi _{{\rm C}, {\rm r}}-\varphi \;_{\rm W} = \varphi _{\rm F}-\varphi _{{\rm sp}}, \;$$

from where the heat flux was estimated as:

(62)$$\varphi _{\rm W} = \varphi _{{\rm C}, {\rm r}}-\varphi _{\rm F} + \varphi _{{\rm sp}}, \;$$

where φ _C,r is the heat flux conducted in the atmosphere through the ridge, φ _W is the heat flux transferred from the water into the ice, φ _F is the heat flux generated from ice freezing and φ _sp is the initial specific heat capacity stored in the ice given as:

(63)$$\varphi _{{\rm sp}} = \rho _{\rm i}V_{\rm i}c_{\rm i}( {T_0-T_{\rm F}} ) , \;$$

where the term V _i is the rubble ice volume, T ₀ is the initial temperature of the ice rubble and c _i is the specific heat capacity of the ice. The change in heat flux per unit time φ _sp/Δt represents the heat flux caused by the stored energy in ice. It was observed that the salinity content above 4 ppt influences the heat capacity of sea ice. The specific heat capacity estimation of the sea ice was presented by Schwerdtfecer (Reference Schwerdtfecer1963) based on Assur's (Reference Assur1960) sea-ice phase diagram. The equation for heat capacity of the sea water (c _sw) is given as:

(64)$$c_{{\rm sw}} = -\displaystyle{S \over {xT^2}}L_{\rm i} + \displaystyle{S \over {xT}}( {c_{\rm w}-c_{\rm i}} ) + c_{\rm i}, \;$$

where c _w and c _i are the specific heat capacities of freshwater and pure ice, respectively. The coefficient of proportionality between bulk salinity S and temperature in °C is x = −0.0182 °C⁻¹. The heat flux transfer from water to ice was estimated to equal 63 and 235 W m⁻² for initial consolidation thicknesses of 0.1 and 0.32 m, respectively (Høyland and Liferov, Reference Høyland and Liferov2005).

After the initial phase of the pressure ridge consolidation, the rubble macroporosity reduces while consolidation progresses downwards. This consolidation depends on the water permeability from the ocean into the macropores between the ice rubble (Shestov and Marchenko, Reference Shestov and Marchenko2016). The influence of initial macroporosity, water temperature and sea-ice salinity in the porosity reduction was estimated based on measurements conducted in the Barents sea (Shestov and Marchenko, Reference Shestov and Marchenko2016) (see Fig. 9). Results showed that the macroporosity changes were smaller for lower water temperatures, which could be a result of a lower temperature gradient between ice and water. In addition, the estimations showed that an increase in the water salinity from 3 to 6 ppt for the same initial macroporosity and water temperature resulted in a greater reduction in porosity. This implies that an increase in salinity will also increase the heat capacity of ice, therefore the consolidation and thus the porosity reduction is higher in the initial phase of consolidation. This result is also supported by Chen and Høyland (Reference Chen and Høyland2016). Finally, a study on the consolidated layer along the cross section of the channel was obtained by Sandkvist (Reference Sandkvist1980). The thickness of the consolidated layer measured in the so-called ‘North Track’ is illustrated in Figure 10. The scatter of the consolidated layer with the cumulative freezing temperature can be considered to be due to the initial porosity in the location of measurements and the frequency of breaking.

Fig. 9. Macro-porosity reduction Δp versus initial porosity p ₀ for water salinity of 3 and 6 ppt. The water temperature is (1) −1.6°C; (2) −0.8°C and (3) −0.2°C. Adapted from Shestov and Marchenko (Reference Shestov and Marchenko2016).

Fig. 10. Thickness of the consolidated layer was measured in a full-scale brash ice channel versus the cumulative freezing degree-days. Adapted from Sandkvist (Reference Sandkvist1980).

The above discussion supports the importance of parameters that influences the initial ice consolidation after each breaking event due to the temperature gradient between ice and water and its dependency on the initial physical properties such as temperature, size and salinity. In addition, the effect of brash ice macroporosity on consolidation and vice versa should be further investigated and related to physical properties such as piece size distribution and brash ice shape.