A review of level ice and brash ice growth models

Abstract

Brash ice forms in harbours and ship channels from frequent ship passages and the resulting freezing–breaking cycles create a unique ice formation. The brash ice accumulation over the winter season is a result of meteorological, thermodynamical and mechanical processes. A reliable brash ice growth model is an important asset when determining navigation routes through ice conditions and when establishing port ice management solutions. This review aims to describe the brash ice development and its modelling as well as the key parameters that influence the brash ice growth and its estimation. This paper summarises the brash ice growth models and the fundamental theories of level ice growth upon which these models are based, and outlines the main knowledge gaps. The results highlight the importance of porosity and piece size distribution and their effect on the consolidation process. The inclusion of the brash ice lateral movement and the side ridge formation would improve the accuracy of forecast models. Furthermore, the findings of the study identify the effect of omitting meteorological parameters such as snow and radiation, from the brash ice growth models. Their contribution to the level ice thickness suggests a significant influence on the brash ice consolidation process.

1. Introduction

Mid- and high-latitude harbours and fairways are liable to severe incidents due to cold weather and ice cover. Events such as icing of structures, ridge formation, ice drifting, ice collar formation and brash ice accumulation, require proper management and control to ensure safe and reliable operations. A discussion of these topics is given by de Almirόn de Andrés and others (2018). Brash ice accumulation caused by vessel transits in harbours and channels often causes time delays, economical losses, hindering ship movements and in some cases, collar formation in the berthing area (see PIANC, 2004). Therefore, winter navigation relies upon management strategies such as systematic observations, forecasting models and physical ice management with a focus on ice-breaking operations. Frequent navigation in harbours and channels accelerates ice accumulation, and shipping operations must depend on the assistance of icebreakers when brash ice accumulates. The formation of thick brash ice increases ship resistance and does not allow ships to always navigate along the centre of the channel. Consequently, the brash ice is not always symmetrically distributed along the channel (Sandkvist, 1986). Eventually, with frequent breaking, navigation in the brash ice becomes harder than that in level ice (Mellor, 1980; Fransson and Sandkvist, 1985).

Brash ice accumulation depends on the variation of the channel geometry caused by the passage frequency and ship hull geometry and propulsion. The inclined bows move the brash ice mainly under the ship's hull, and thus it goes back into the channel, while if the bow is wedge-shaped ice is pushed mainly sideways (Kitazawa and Ettema, 1985; Ettema and others, 1998). The blocks of ice that are pushed sideways can pile up under the adjacent level ice and form side ridges. The ship channels can be divided into narrow or wide channels. In narrow channels, brash ice experiences higher confinement, which can increase the thickness of brash ice. In narrow channels, the side ridges are also more prominent. The depth of the water is also an important parameter, especially for harbours, close to shore, rivers and lakes where the water is shallower, and the accumulation of brash ice and side ridges can extend vertically down to the seabed. In wider channels such as in harbour basins, the confinement is not that prominent nor the side ridges (Riska and others, 2019; Bridges, 2020).

Thermally grown ice cover contributes to the weather by insulating the ocean from the atmosphere. In particular, the development of fast ice during winter influences the local climate. Vice versa, the weather conditions, i.e. low air temperatures, wind and precipitation, influence the ice growth and brash ice consolidation. Therefore, the ice thickness development and accuracy of forecast models depend on different meteorological conditions in mid- and high-latitudes (Ohata and others, 2016). While thermodynamic models have been developed, differences between predicted and measured values are present, e.g. caused by snow (Cheng and others, 2014). The snow contribution to ice growth has been previously studied for coastal landfast sea ice by Leppäranta (1983), Saloranta (2000), Shirasawa and others (2005) and Wang and others (2015), for lake ice by Leppäranta and Kosloff (2000), Cheng and others (2014) and Ohata and others (2016), lake and river ice by Ashton (2011) and for Antarctica sea ice by Crocker and Wadhams (1989), Adolphs (1998) and Fichefet and Maqueda (1999). From these studies, it can be deduced that one of the main challenges when estimating ice thickness is accounting for the snow-ice formation caused by the interface flooding after heavy snowfall. This implies that the time of slush formation and slush thickness development do not always correspond to the assumption of instant flooding based on Archimedes' principle. The slush formation depends on the existence of fractures that enable water to rise and flood the snow, and also the slush thickness which depends on the snow's ability to absorb water (capillarity). These phenomena are not usually quantified and thus neglected from the ice growth models.

Thermal and physical properties of ice, i.e. thermal conductivity, heat capacity, density and microstructure are affected by the entrapped brine inclusions. The brine content in ice was found to be a linear function of the ice growth rate (Weeks and Lofgren, 1967; Nakawo and Sinha, 1981). However, observations indicate that this does not apply to brackish water (Granskog and others, 2006). In addition, radiation, temperature gradients and salinity are affected by the ice thickness (Maykut, 1978). Furthermore, modelling the ocean–atmospheric interaction is essential to the initial ice growth in the early winter (Gabison, 1987).

Brash ice consolidation between two consecutive vessel passages occurs analogously to the level ice growth, with the addition of the macroporosity and the initial heat capacity of ice blocks that contribute to the consolidation. Stefan's (1889) approach was used earlier to forecast the brash ice volume increase between two consecutive ship passages, which is a consequence of brash ice consolidation and breaking (Ashton, 1974; Berenger and Michel, 1975; Sandkvist, 1980, 1981, 1982, 1986; Vance, 1980). These models relate the volume increase of brash ice only to the cumulative freezing degree-days and the breaking frequency. Effects of important parameters such as the lateral movement of ice, macroporosity and snowfall are neglected, which can lead to brash ice thickness overestimation. An attempt was made to include the lateral movement of the ice pieces by introducing coefficients for expelled ice blocks (Ettema and Huang, 1990). However, parametrisation and estimation of such coefficients require more observations and measurements. Improved theoretical models were reported by Chomatas (2015), Riska and others (2014, 2019) and Karulin and others (2019). The side ridge formation, radiation and snow cover influence on the brash ice development are not fully considered in these models. In addition, the initial brash ice formation and growth when the channel is not completely covered by ice have not been modelled. As such, the accuracy and reliability of existing brash ice models can be improved. Systematic observations and measurements in ports and channels are necessary to validate and improve these models.

The current study aims to identify important parameters that influence the brash ice development and effects of the accuracy of the brash ice growth models. This study initially provides an up-to-date review of the analytical level ice growth models (Section 2) which gives the foundation on the physical, thermal and meteorological parameters that are crucial for ice growth predictions. This is followed by a review of the brash ice growth models and physical properties that affect the brash ice consolidation (Section 3). The study concludes with a brief outline of knowledge gaps concerning the brash ice development and brash ice prediction models (Section 4). Finally, the nomenclature for the symbols used in this paper is summarised (Appendix).

2. Level ice formation and growth

Ice formation and thickness changes of freshwater and sea-water ice occur in the ice/water and ice/air or snow/air interfaces, which are usually defined as lower and upper boundaries. During freshwater ice formation and growth, the temperature equilibrium between ice and water is preserved by ice formation or melting (Stefan, 1889). Initially, when the air temperature is below the freezing point, the water gets colder until the top temperature of water reaches the freezing point and subsequently, the water-ice phase change starts to take place. The maximum density of fresh water at 4°C keeps some of the water column above freezing point and gives rise to oceanic heat flux. The temperature gradient between the top and bottom boundaries causes the heat generated by ice formation, and oceanic heat to be conducted upwards into the atmosphere. The temperature gradient between boundaries is mainly affected by the air temperature. A comprehensive summary of ice formation for fresh water is given by Ashton (1986).

During sea-water ice formation and growth at the water/ice interface, the salt in sea water is partially excluded from the ice structure and partially entrapped in the ice (Bennington, 1963). The amount of salt that is trapped in the ice varies with the freezing rate. A temperature decrease causes ice formation and an increase of salt concentration within brine inclusions, while an increase in temperature causes the melting of the surrounding ice and dilution of the brine (Untersteiner, 1961). The sea-ice surface can also undergo accretion due to snowfall and melting upon downwelling radiation and air temperature increase (Maykut and Untersteiner, 1971).

The thermal growth of ice was predicted by Stefan (1889) and this formulation has been used as a basis for deriving equations that include the ice–atmosphere coupling (Ashton, 1989), dry snow layer as an insulator, and snow–atmosphere coupling (Nakawo and Sinha, 1981). These cases are shown in Table 1(a)–(c) respectively. Additionally, the snow-ice formation in the snow/ice interface is induced by the isostatic equilibrium between snow mass and ice buoyant weight pushing the ice top surface under the water level. The term ‘ice’ refers here to congelation ice formed from water freezing. When the snow layer weight exceeds the ice buoyant capacity, the ice/snow interface is submerged, and water is seeping on top of the ice turning some snow into slush. This process is described in detail later. Subsequent slush to snow-ice transformation was found to be important when estimating the snow and ice thickness especially considering areas and winters with high-snow precipitation, i.e. the Baltic Sea (Leppäranta, 1983). These models including the snow ice are illustrated in Table 1(d) and (e), respectively.

Table 1. Schematic illustration of level ice growth model approaches

The analytical level ice growth models include: the heat flux conducted through the ice, snow ice and snow layer denoted by, respectively, φ _I, φ _SI and φ _S, the convective heat flux that flows through the ice or snow to the atmosphere denoted with φ _IA and φ _SA. Heat flux generated due to water and slush freezing are, respectively, φ _F and φ _FSL. In each interface, the incoming and outgoing heat fluxes are equal, assuming continuity in the heat fluxes. The numerical model considers the sensible and latent heat flux φ _SE and φ _L, downwelling and upwelling shortwave and longwave radiation φ _SW↓, φ _SW↑, φ _LW↓, φ _LW↑ and the conductive heat fluxes through ice and snow layer.

Furthermore, the net radiation heat flux was introduced in the growth models by Maykut and Untersteiner (1971) and thermal inertia by solving the conduction partial differential equations with finite-element techniques (Launiainen and Cheng, 1998). A partially numerical and analytical thermodynamic model was initially tested by Toyota and Wakatsuchi (2001), which was focused on the initial phase of ice growth. A schematic illustration of the heat fluxes that are usually included in numerical estimations is given in Table 1(f).

2.1 Stefan's ice growth formulation

Stefan's (1889) 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, 2007). 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, 1963). Pure ice conductivity k _pi in W m⁻¹ K⁻¹ as a function of temperature in K is given by Yen (1981) 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 (1963) 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, 1960). 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 (1964), and also by Pringle and others (2007).

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 (2013) 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, 1986, 1989; Ettema and Huang, 1990; Andres and van der Vinne, 2001). 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, 1986, 1989). 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 (1989) equation (Ashton, 1986, 1989). 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 (1981), Leppäranta (1983), Ashton (1986), Leppäranta (1993), Andres and van der Vinne (2001), Comfort and others (2003), Comfort and Abdelnour (2013). 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 (1977), Sturm and others (1997), Calonne and others (2011) and Riche and Schneebeli (2013). The most commonly used snow conductivity k _s (in units of W m⁻¹ °C⁻¹) in ice growth models is given as (Yen, 1981):

(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 (2019) 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, 2003). For example, the maximum snow thickness of 140 mm observed in a prairie pond by Andres and van der Vinne (2001) 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, 1962).

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, 1962). The surface ice growth caused by flooding was predicted by Leppäranta (1983). 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, 1983). 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 (2000), Shirasawa and others (2005), Cheng and others (2014) and Wang and others (2015).

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 (1998). 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, 2003, Comfort and Abdelnour, 2013). 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, 2000; Shirasawa and others, 2005; Ohata and others, 2016). β 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, 2000). 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, 2016).

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, 2000). 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, 1983; Saloranta, 2000), as well as the availability of the cracks or other water routes through the ice layer (Ashton, 2011). 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, 2000), as illustrated in Figure 1. However, it is important to mention that the ice thickness presented by Saloranta (2000) 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 (2000).

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, 2014). 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 (1983) and Saloranta (2000), 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 (2014).

A partly numerical approach has been introduced recently aiming to study the effect of snowfall on the initial ice growth stage (Toyota and others, 2020). 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, 1961, 1964), 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, 1992; Hori and others, 2006; Cheng and Liang, 2014). 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, 1994).

The sky emissivity depends on vertical temperature profile, humidity, vapour pressure, precipitated water and cloud occurrence (Niemelä and others, 2001). 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, 1981). 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, 1990). 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, 1998, 2002a, 2002b). 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, 1996). 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, 2006). 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, 1996) and the bulk ice density (Jin and others, 1994). 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, 2001; Toyota and others, 2020). The sensible and latent heat fluxes can be determined from equations proposed by Maykut (1978):

(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, 1978).

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, 1978). Consequently, the heat fluxes suggested by Maykut (1978) for an ice growth from 0 to 0.2 m are considerable. In addition, Toyota and others (2020) 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, 1987).

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, 1997; Ohata and others, 2016):

(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, 1982):

(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 (2016):

(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, 2016), 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, 1995). The seasonal variation of the oceanic heat flux was suggested to be important when estimating the annual thermodynamic ice growth (Wang and others, 2015). 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, 2017), 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, 2014). 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, 2018). 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 (2017).

3. Brash ice formation and growth

The lifespan of brash ice consists of three phases: formation, growth and melting. Furthermore, the brash ice formation phase may comprise two different events: first, the formation from breaking the level ice, or second the brash ice initiation from vessel passages starting at the beginning of ice formation when brash ice grows from slush.

The first type of brash ice formation occurs when a vessel initially navigates in thick unbroken level ice. The first ship passage breaks the level ice into sharp-edged ice floes. This initial breaking forms a vessel track that consists of ice blocks floating in the water. The channel initially has a width slightly larger than the vessel's beam and comprises many open water patches. The time between two consecutive ship passages allows the thermodynamic ice growth (consolidation) to occur at the open water where ice growth is similar to level ice growth. The consolidation (growth) of the broken ice pieces is initially driven by the temperature gradient of the submerged ice with water and continues until the temperature equilibrium is reached. The subsequent growth is driven by the temperature gradient between the ice/atmosphere and ice/water interfaces. The same consolidation phenomenon is observed in the ridge formation (Høyland and Liferov, 2005) where the ice initially consolidates due to the stored energy (the specific heat capacity of ice), followed by the thermodynamic growth driven by the temperature gradient between boundaries. However, contrary to ridge formation, the brash ice consolidation is a repeated process after each ship passage. This consolidation process freezes the ice pieces into a solid ice layer. The following vessel passages in the same track break the consolidated layer, displace and erode the ice pieces (Ettema and Huang, 1990). Consequently, the size and shape of the ice pieces are modified after several vessel passages into rounded but still floe-shaped ice pieces (Greisman, 1981; Sandkvist, 1986). Transiting ships displace the brash pieces below the adjacent level ice and overtime these form side ridges (Kitazawa and Ettema, 1985; Sandkvist, 1986; Ettema and others, 1998). The piece size of the ice floes reduces with frequent breaking. After several ship passages, new ice pieces are generated from the broken consolidated ice, and the existing ice floes have eroded, and eventually, the brash ice covers the whole water surface. This brash ice formation is common in fairways which are made later in the winter, and often as an alternative to the existing fairways with high brash ice accumulation or in regions with less shipping traffic. A schematic illustration of this process is shown in Figure 4, and a discussion of the brash ice formation and growth is given by Ettema and Huang (1990).

Fig. 4. Brash ice formation and accumulation scheme showing: (1) the first ship passage in level ice (LI) and the vessel's beam (VB), (2) the width of brash ice channel (BICh) and side ridge (SR) formation phase, (3) the end of brash ice formation phase where the channel is fully covered with ice pieces and (4) the brash ice and side ridge accumulation and thickening. Adapted from Ettema and Huang (1990).

The second type of brash ice formation begins with frequent navigation at the start of winter. When the air temperature drops below the water-freezing temperature, ice crystals form and suspend on the water surface. Over time, the initial ice crystals consolidate into a thin layer of ice. The initial ship passages in the early stage of ice forming agitate the ice with water and form a slush layer. The snowfall may add to the mass of the initial slush formed (Toyota and others, 2020; Zhaka and others, 2020). The heat loss due to cold air temperature and slush ice agitation initiate consolidation of the slush (Sandkvist, 1982). Initially, the slush formed at each ship passage does not cover the whole channel and ice growth is a mixture of open water freezing and slush consolidation. When the slush gets thicker, the whole channel becomes covered. However, there is no published information on when this initial phase ends, i.e. the point when the whole channel is ice covered. Further breaking and refreezing leads to the accumulation and growth of slush into brash pieces to the point when the ship channel is fully covered with a layer of brash ice pieces.

After both types of brash ice formation, the second phase of development consists of the thickening of the brash ice layer after the surface of the brash ice channel is fully covered. After each vessel passage, the brash ice is partly redistributed in the channel and partly displaced sideways. Ice that moves sideways builds up under the level ice edge, thus increasing the side ridge thickness. The lateral confinement from the adjacent side ridges and level ice force the brash ice to accumulate in the channel's centre (Greisman, 1981, Sandkvist, 1982). This accumulation and packing of ice blocks create voids in the brash ice layer, so-called macroporosity which is the fraction of water volume to ice volume. Theoretically, the macroporosity can be considered as the sum of the individual macropore volumes. The macropores can be either water-filled voids surrounded by the ice blocks after the brash ice displacement stops or air-filled which are commonly formed on the channel's surface (in the freeboard above the water level). Between two vessel passages, the consolidation occurs initially in the open water fraction in the brash ice surface (if present), and subsequently the freezing proceeds downwards by freezing the water-filled macropores. Macroporosity is to be distinguished from microporosity which usually means the ratio of air bubbles (and other additions) to the total volume of solid ice. Repeated vessel passages break the consolidated layer which generates an additional volume of ice and thickens the accumulated brash ice layer, while also changing the macroporosity.

Melting of the brash ice channels is the last phase of the annual brash ice lifespan. The melting phase can be faster in a frequently navigated fairway as a consequence of the turbulent heat fluxes caused by ship passages and an increase in heat transfer due to the water-ice agitation.

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 (1974). 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, 1974). 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 (1980, 1981). 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, 1983). 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, 1990; Riska and others, 2019).

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, 1990). 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, 1990). 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 (2019). 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 (2019).

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 (1974) and Sandkvist (1980, 1981) 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 (2019).

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, 1980). 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 (1980). Adapted from Sandkvist (1980, 1986).

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, 1985). A more recent study of a model scale brash ice channel observed porosity values between 40 and 60% (Bridges and others, 2020). 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, 2018). Higher values reported by Høyland (2007) reached 45%. However, the rubble ice forms naturally from compressing or shearing of level ice (Høyland, 2007), 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 (2021). 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 (2019) and Matala (2021). 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, 2019), this distribution trend is also supported by Tuovinen (1979), 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, 1982).

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, 2002). 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, 2005). 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, 2001; Høyland and Liferov, 2005; Chen and Høyland, 2016; Chen and others, 2020).

The influence of the rubble thickness in the initial ice consolidation was investigated by Høyland and others (2001). 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 (2016). 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 (2020). 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 (2016) and Chen and others (2020). 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, 2016). 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 (2016).

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, 2005). 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 (1963) based on Assur's (1960) 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, 2005).

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, 2016). 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, 2016) (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 (2016). Finally, a study on the consolidated layer along the cross section of the channel was obtained by Sandkvist (1980). 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 (2016).

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 (1980).

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.

4. Summary of knowledge gaps

This review investigated the level ice and brash ice growth models and parameters that influence the growth and accumulation of ice. Level ice analytical growth models derived from Stefan's approach were analysed for level ice and level ice covered with dry snow. The slush formation from snow flooding, its transformation to snow ice and its contribution to the total level ice thickness were all investigated. The influence of surface heat budget on the surface temperature and ice growth was evaluated including the net radiation, sensible and latent heat fluxes. The oceanic heat flux convected in the ice/water interface and its contribution to the annual ice growth was also analysed. It was found that snow physics, in particular, the slush-to-ice transformation process has a major contribution to the differences between measured and predicted ice thickness values. The main sources of model errors were related to accurately estimating the time of flooding, slush and snow-ice formation. In addition, shortwave radiation in spring and summer initiates the melting of snow followed by a snow-ice formation which retards the ice growth. The shortwave radiation absorbed in the ice causes ice ablation. The freeze-up of the water at the beginning of the winter season depends on the oceanic heat flux which is partly caused by the heat accumulated during summer from shortwave radiation. In addition to oceanic flux, the sensible and latent heat fluxes are observed to be an important parameter that contributes to the freeze-up of the open water and the initial growth of level ice.

In addition to thermodynamical and meteorological phenomena that influence level ice growth and brash ice consolidation, the brash ice formation and accumulation in ice-infested waters is affected by the mechanical process of frequent icebreaking. A description was given for the brash ice formation initiating from breaking the thick level ice or from initially navigating in very thin ice skim. This included the description of brash ice growth and the accumulation of brash ice caused by the frequent vessel navigation and the thermodynamic consolidation. The study comprised of a description of the brash ice growth models and the important parameters that contribute to brash ice development, including the lateral movement of brash ice and the side ridge formation, macroporosity, piece size distribution and consolidation. However, the existing brash ice growth models are based on few parameters such as the cumulative freezing degree-days, icebreaking frequency and macroporosity. Parameters included in the level ice models such as snow, shortwave radiation, sensible, latent and oceanic heat fluxes are not included in the existing brash ice models and their contribution to the brash ice development is not supported by measured data for verification. The effect of snow as an insulator between two consecutive passages and its effect on adding to the volume of brash ice after each ship passage has not been previously studied, and consequently, the influence of snow is neglected from the growth models. Considering the effect that these parameters have on the level ice development, it is suggested that i.e. the snow, and in particular the slush-to-ice transformation process can affect the brash ice consolidation in locations and winters when the amount of snowfall is significant. The latest model developed by Riska and others (2019) divided the brash ice model into three layers characterised by different porosity and thermal properties; the dry brash ice, the consolidated layer and the wet brash ice layer. A further improvement in the model should include the snow layer on the dry brash ice layer. Their model validation with brash ice thickness values from Sabetta port (Arctic) showed a higher accuracy compared to model validation with brash ice thickness data from Luleå port (Sub-Arctic). The highest deviation from the measurements in Luleå was observed in mid-March to April, which can be explained by the possible influence of the shortwave radiation that was omitted from the model. The effect of shortwave radiation and surface heat budget on the consolidation of the brash ice should also be further investigated and incorporated into future ice growth models.

The brash ice measurements are mostly obtained from full scale and model scale brash ice channels formed from breaking an initial level ice layer of thickness H ₀. Little is known about the brash ice's initial formation from ships' navigation channels in early winter, and the importance of this initial stage has not been previously addressed, i.e. the transformation of slush to brash ice, as well as the effect of snow on this process. Full-scale observation and measurement in the early winter can give insights into the importance of this initial formation phase. After each ship transit, the temperature of the ice blocks that submerge in the freezing water varies with the cumulative freezing temperatures between vessel passages. The temperature difference between the water and ice blocks causes a temperature gradient that is followed by ice formation until the temperature equilibrium between the two phases (water and ice) is reached. The time required to reach the temperature equilibrium after each vessel passage and the quantification of this porosity reduction require further investigation and research. Even though this phenomenon is addressed in previous models (Ettema and Huang, 1990; Riska and others, 2019), there is limited data to support the parametrisation of the phenomenon. Furthermore, after each vessel passage, a part of the brash ice is pushed sideways. Most of the brash ice forecast models assume that all ice remains in the channel which theoretically could lead to an overestimation of brash ice layer thickness. This is also concluded from the validation of brash ice growth models described previously by Ettema and Huang (1990) and Riska and others (2019). The volume fraction of the brash ice that is pushed from the channel at each vessel passage was previously included in the brash ice accumulation forecast by Ettema and Huang (1990); however, the accurate determination of this fraction requires systematic full-scale measurements.

To conclude, a better understanding is required for the initial brash ice formation phase, the amount and influence of brash ice lateral movement, the initial consolidation after each breaking event, the porosity variation during the brash ice lifespan, the snow and effect of radiation on brash ice formation and growth. These phenomena should be studied further, along with full-scale and model-scale observations and measurements, and be introduced in the brash ice growth models.