Articles | Volume 18, issue 4
Research article
09 Apr 2024
Research article |  | 09 Apr 2024

Using Icepack to reproduce ice mass balance buoy observations in landfast ice: improvements from the mushy-layer thermodynamics

Mathieu Plante, Jean-François Lemieux, L. Bruno Tremblay, Adrienne Tivy, Joey Angnatok, François Roy, Gregory Smith, Frédéric Dupont, and Adrian K. Turner

Icepack (v1.1.0) – the column thermodynamics model of the Community Ice CodE (CICE) version 6 – is used to assess how changing the thermodynamics from the Bitz and Lipscomb (1999) physics (hereafter BL99) to the mushy-layer physics impacts the model performance in reproducing in situ landfast ice observations from two ice mass balance (IMB) buoys co-deployed in the landfast ice close to Nain (Labrador) in February 2017. To this end, a new automated surface retrieval algorithm is used to determine the in situ ice thickness, snow depth, basal ice congelation and snow-ice formation from the measured vertical temperature profiles. Icepack simulations are run to reproduce these observations using each thermodynamics scheme, with a particular interest in how the different physics influence the representation of snow-ice formation and ice congelation. Results show that the BL99 parameterization represents well the ice congelation but underrepresents the snow-ice contribution to the ice mass balance. In particular, defining snow-ice formation based on the hydrostatic balance alone does not reproduce the negative freeboards observed for several days in the IMB observations, resulting in an earlier snow-flooding onset, a positive ice thickness bias and reduced snow depth variations. We find that the mushy-layer thermodynamics with default parameters significantly degrades the model performance, overestimating both the congelation growth and snow-ice formation. The simulated thermodynamics response to flooding, however, better represents the observations, and the best results are obtained when allowing for negative freeboards in the mushy-layer physics. We find that the mushy-layer thermodynamics produces a larger variability in congelation rates at the ice bottom interface, alternating between periods of exceedingly fast growth and periods of unrealistic basal melt. This pattern is related to persistent brine dilution in the lowest ice layer by the congelation and brine drainage parameterizations. We also show that the mushy-layer congelation parameterization produces significant frazil formation, which is not expected in a landfast ice context. This behavior is attributed to the congelation parameterization not fully accounting for the conductive heat flux imbalance at the ice–ocean boundary. We propose a modification of the mushy-layer congelation scheme that largely reduces the frazil formation and allows for better tuning of the congelation rates to match the observations. Our results demonstrate that the mushy-layer physics and its parameters can be tuned to closely match the in situ observations, although more observations are needed to better constrain them.

1 Introduction

The sea ice and oceanography of the Canadian Arctic are largely modulated by the formation of landfast ice in fjords, along the coasts and in narrow channels. Each winter, this landlocked sea ice transforms the Canadian Arctic Archipelago (CAA) into a seasonal continent of stationary sea ice (Melling2002; Galley et al.2012), effectively insulating the seawater from the cold atmosphere and barring the transport of ice through the CAA passages (Howell et al.2013; Kwok2006). The landfast ice edge represents a seasonal boundary where the air–ocean exchanges and ice dynamics processes are concentrated, in particular by the opening of semipermanent polynyi under divergent surface forcing conditions (Melling et al.2001; Dumont et al.2010). These flaw polynyi in turn drive the regional meteorology (Barber et al.2001; Gultepe et al.2003; Lüpkes et al.2008; Raddatz et al.2011) and ocean circulation (Dumont et al.2010), producing sediment-rich waters that are key to the Arctic marine ecosystem (Stirling1980, 1997; Carmack and Macdonald2002; Tremblay et al.2002). As changes in the landfast ice cover are expected to alter these processes, its monitoring, representation in forecast models and inclusion in climate projections are a concern not only for the study of the Arctic climate but also for a wide range of socio-economical aspects such as on-ice transport safety, food security and navigation planning (Gearheard et al.2006; Eicken et al.2011; Cooley et al.2020).

In dynamic sea ice models, the physics of landfast ice is represented using a combination of thermodynamic relations governing the ice growth and melt (i.e., a column thermodynamics model; Maykut and Untersteiner1971; Semtner1976; Bitz and Lipscomb1999; Huwald et al.2005; Turner et al.2013) and of dynamic parameterizations governing its stability against external forces (i.e., a rheological model; Hibler1979; Hunke and Dukowicz1997; Tremblay and Mysak1997; Wilchinsky and Feltham2004; Rampal et al.2016). While these components are mostly treated (and developed) independently, they remain deeply interconnected and the formation of landfast ice usually results from their combined action. In many areas, for instance, the landfast ice is held by the grounding of ice keels on the ocean floor, which involves prior ridging (dynamics) of sufficiently thick ice (thermodynamics). In the absence of ice grounding, landfast ice can form during periods of calm and cold weather (Divine et al.2004; Kirillov et al.2021), during which leads freeze to a sufficient ice thickness for the unconsolidated ice floes to coalesce together (thermodynamics), allowing the formation of ice arches between pining points that resist subsequent surface forcings (dynamics;  Dammann et al.2019; Liu et al.2022). In sea ice models, this interplay between thermodynamic and dynamic factors is represented by ice thickness dependencies in the dynamic parameters, such as the seabed stress term (Lemieux et al.2015) or the material strength parameters (Dumont et al.2009; Lemieux et al.2016; Plante et al.2020; Liu et al.2022). The accurate representation of landfast ice extent, trends and variability in sea ice models therefore requires not only the permitting dynamics (i.e., ice grounding and tensile strength) but also thermodynamics that reproduce well the landfast ice growth and melt.

In the Environment and Climate Change Canada (ECCC) ice–ocean forecasting systems (e.g., RIOPSv2; Smith et al.2021), the implementation of the aforementioned landfast ice dynamics was shown to greatly improve the representation of landfast ice in hindcast (free-run) simulations (Lemieux et al.2016). The timings of landfast ice formation and breakup, however, remain difficult to reproduce, often offset by a couple of weeks with respect to those recorded in operation ice charts (Lemieux et al.2016). While this could be improved by modifications to the ice grounding mechanics (e.g., Dupont et al.2022) or by changes to the ice strength formulation (Ungermann et al.2017), it is also possible that the discrepancy is associated with a misrepresentation of the landfast ice thermodynamics, which in the ECCC systems is based on the model of Bitz and Lipscomb (1999) (hereafter BL99). Thermodynamics models have grown in sophistication over the years, in particular with the representation of mushy-layer physics (Feltham et al.2006), brine dynamics (Notz and Worster2009; Turner et al.2013) and melt ponds (Flocco et al.2010; Holland et al.2012; Hunke et al.2013). These developments are implemented in the Los Alamos Community Ice CodE version 5 (CICE5) and were shown to have competing effects on the overall pan-Arctic ice mass balance, both in long-term global simulations (Turner and Hunke2015) and in coupled climate simulations (Community Earth System Model version 2;  Bailey et al.2020; DuVivier et al.2021). The use of the mushy-layer physics was in particular shown to produce larger amounts of frazil and snow ice, together increasing the overall ice thickness. Whether this increase is also seen in the landfast ice context (without sensitivities to the offshore sea ice dynamics) remains to be determined.

In recent years, the deployment of ice mass balance (hereafter IMB but used here as a general term not referring to specific designs) buoys in both the Arctic and Antarctic provided in situ observations of the thermodynamics in the sea ice interior by measuring the internal sea ice temperature at a high vertical (cm) and temporal (h) resolution (Richter-Menge et al.2006; Jackson et al.2013; Planck et al.2019). The snow depth and ice thickness conditions are inferred from the recorded vertical temperature profiles, traditionally by visual inspection (Tian et al.2017; Provost et al.2017) but more recently by using automated algorithms (Liao et al.2019; Cheng et al.2020; Richter et al.2023). These measurements give new insights into thermodynamic processes that are otherwise not detectable by traditional ice thickness measurements, ice core analyses or remote sensing, such as snow-ice formation (Provost et al.2017; Rösel et al.2018), heat fluxes within and between the material interfaces (Trodahl et al.2000; West et al.2020), brine convection, and mushy-layer properties (Wongpan et al.2018). IMB buoys have also been used to assess the performance of thermodynamics models in the context of 1D simulations (Caixin et al.2015; Tian et al.2017; Duarte et al.2020). The mushy-layer physics in CICE version 5, for instance, has been tested against IMBs deployed in the pack ice (first year and multiyear) during the N-ICE2015 expedition north of Svalbard, where it was shown to adequately represent the observed sea ice growth but also to overrepresent snow flooding under large snow depth conditions (Duarte et al.2020).

In this study, we investigate how updating the model thermodynamics from the BL99 physics to the mushy-layer parameterization impacts the simulated sea ice mass balance in a landfast ice context, away from the pack ice dynamics. This assessment is based on the in situ observations from two IMB buoys that were deployed in a landfast ice channel well sheltered from offshore dynamics, close to Nain (Nunatsiavut, Labrador). A particular interest is placed on the ice growth from congelation and snow-ice formation, which is determined from the recorded internal temperature profiles using a novel surface retrieval algorithm building on the work of Liao et al. (2019) and Cheng et al. (2020). Multiple Icepack (v1.1.0) model simulations (Hunke et al.2023) are run to reproduce these observations using the BL99 physics or the mushy-layer physics to determine the effect of the brine physics on the model's performance. In particular, we find that the use of the mushy-layer physics with default parameters significantly degrades the model performance despite the improved representation of flooding and brine processes. The basal ice growth in mushy simulations is overrepresented, includes a significant contribution from frazil production and exhibits unexpected periods of basal melt. The snow-ice formation is also overrepresented due to early snow flooding when observations are under negative freeboard conditions. We show that these discrepancies are largely resolved by simple modifications and tunings of the congelation and snow-ice parameterizations. The contributions of this paper include a modified mushy-layer congelation parameterization not conducive to frazil formation and a parameterized dependency of the snow flooding rate on negative freeboard values.

This paper is organized as follows. A description of the buoys and surface forcing data used in the analysis is provided in Sect. 2. The Icepack model physics are briefly presented in Sect. 3, first describing the BL99 physics currently used in the ECCC forecast systems and then the differences when using the mushy-layer thermodynamics. Our modifications to the congelation and snow-ice parameterizations are also included in Sect. 3. The methods are detailed in Sect. 4, including the surface retrieval algorithm, the numerical simulation setup and the model performance diagnostics. Results from the in situ observations and Icepackv1.1.0 simulations are presented in Sect. 5. Discussions of the model performance and conclusions are summarized in Sect. 6.

Figure 1Location of the two snow ice mass balance apparatus (SIMBA) buoys on the Labrador coast (a) and in a landfast ice channel close to the Nain community (b). The buoys are located at ∼56.42° N, 61.7° W (SIMBA1) and ∼56.43° N, 61.50° W (SIMBA2), 12 km from each other and ∼50 km from the nearest landfast ice edge. Images are corrected reflectance imagery taken from MODIS Worldview (, last access: 7 September 2023).

2 Data

2.1 Ice mass balance buoy observations

Two Scottish Association for Marine Science (SAMS) snow ice mass balance apparatus (hereafter SIMBA) buoys were deployed in winter 2017 as part of an ongoing collaboration with the Nunatsiavut Research Centre (NRC), with the goal of serving the Nain community with the deployment of scientific instruments in the local landfast ice. The buoys were thus not deployed as part of a wider scientific field observation campaign: the deployment dates and locations were chosen with NRC collaborators based on their sea ice monitoring interests. The first buoy (SIMBA1) was deployed on 23 February 2017 at ∼56.42° N, 61.7° W in a landfast channel close to the southern coast of Satosoak Island (see Fig. 1) and recovered 2 months later on 18 April. The second buoy (SIMBA2) was deployed during the same season on 24 February at ∼56.43° N, 61.50° W, ∼12 km east of SIMBA1 in the same fjord close to Palungitak Island and recovered 3 months later on 31 May. To our knowledge, this was the first time IMB buoys were deployed in this area.

Figure 2Schematics of the deployed SIMBA buoy thermistor strings through the snow, snow-ice and sea ice layers (a) and the vertical temperature profiles they measure, with the sensor positions identified by the surface retrieval algorithm (b). Note the section of the thermistor string (thermistor plateau, red lines) laid flat on the bare ice surface at deployment but later embedded within the ice layer after flooding. The marked sensor locations are, from top to bottom as follows: the air–snow interface (Za−s), the snow interior (Zs), the snow–ice interface above the snow ice (Zs−i), the snow–ice interface at deployment (Zice0), the lower end of the thermistor plateau (Zp), a point in the ice interior close to the ice–ocean interface (Zc) and the ice–ocean interface (Zi−o).


The SIMBA buoys consist of a 5 m long thermistor string with temperature sensors (Maxim DS28EA00, with a 0.0625 °C resolution and 0.0625 °C accuracy) placed every 2 cm (Jackson et al.2013). The thermistor strings are deployed vertically through a 5 cm hole such that the sensors measure the vertical temperature profile from the atmosphere above the snow layer down to the seawater below the ice (Fig. 2a). At deployment, a section of the thermistor string is laid flat on the ice surface to mark the initial snow–ice interface in the data (see red arrows and dashed lines in Fig. 2a–b). The sensors within this thermistor string section are thus all at the same depth and show nearly identical temperature readings, making this segment easily identifiable in the vertical temperature profiles. The hole is then refilled with slush and the snow cover is carefully restored to its original depth. The vertical temperature profiles are measured with a 6 h time resolution and are transmitted remotely via Iridium satellite along with the recorded air temperature, atmospheric pressure and GPS location.

The SIMBA buoys also perform daily heat cycle measurements, which consist in recording the temperature change associated with a 1 min and 2 min heating from a resistor component besides each temperature sensor (Jackson et al.2013). This change in temperature can be used to infer the heat capacity and conductivity of the medium surrounding the sensors, and it is used in this study to visually locate the material interfaces and validate the accuracy of our surface retrieval algorithm.

2.2 GDPS atmospheric forcing

Data from the ECCC Global Deterministic Prediction System (GDPS; Buehner et al.2015; Smith et al.2018) are used to compute the atmospheric fluxes driving our thermodynamic simulations at the air–snow interface. The GDPS was previously shown to be equally representative of observations as more commonly used reanalysis data (Smith et al.2014), and it offers an accurate estimate of the atmospheric conditions in our study region with limited surface observations.

The GDPS is a coupled atmosphere–ice–ocean forecasting system using the Global Environment Multiscale (GEM) model for the atmosphere (Côté et al.1998b, a), the Los Alamos multicategory Community Ice CodE (CICE) model version 4 for the sea ice (Hunke et al.2010) and the Nucleus for European Modelling of the Ocean (NEMO) model for the ocean (Madec et al.1998; Madec and the NEMO team2008). This system produces 10 d forecasts with 3 h outputs of the atmosphere, ice and ocean, initialized each day at 00:00 UTC with fields from a data assimilation system (e.g., a four-dimensional ensemble-variational data assimilation scheme for the atmosphere; for details see Buehner et al.2013, 2015). Here, we use the archived surface fields from the 006–027 h forecasts (i.e., after a 6 h spin-up) to drive the atmospheric fluxes in our model. At these very short lead times, only limited deviations from the initial analysis fields are expected (Smith et al.2014). The GDPS variables used in our analysis include surface winds, air temperature, humidity, shortwave and longwave radiation, and precipitation, all taken at the grid point location closest to the buoy deployment.

3 Model

Column thermodynamics simulations are produced using Icepackv1.1.0, which is the thermodynamics package from CICE6. This package corresponds to a collection of thermodynamics parameterizations that can be chosen by the user. In this analysis, we use Icepack with two different thermodynamics schemes: simulations are first ran using the BL99 thermodynamics, available in CICE version 4 and employed in the ECCC systems, and then repeated using the mushy-layer thermodynamics, available from CICE version 5 onward. All simulations share the same forcing (atmosphere and ocean fluxes) and snow model, but the mushy-layer thermodynamics includes improvements in the representation of brine processes and modifications to the sea ice congelation and snow-ice formation parameterizations (Turner and Hunke2015; Bailey et al.2020).

3.1 Standard BL99 thermodynamics

3.1.1 Surface thermodynamic balance

The thermodynamic growth and melt of sea ice are governed by the net energy balance at the top and bottom ice (or snow) surfaces. At the top interface, the atmospheric fluxes are calculated from the GDPS data, and the net heat flux F0 (positive downward) at the top interface is written as

(1) F 0 = F s + F l + F LW + ( 1 - α ) ( 1 - i 0 ) F SW ,

where Fs is the sensible heat flux, Fl is the latent heat flux, FLW is the net longwave flux, α is the surface shortwave albedo, i0 is the fraction of shortwave penetration into the ice or snow surface, and FSW is the incoming shortwave flux. In all simulations, the shortwave albedo and penetration are defined by the Community Climate System Model version 3 (CCSM3; Collins et al.2006).

Due to the absence of ocean salinity and ocean currents observations at the buoy locations, no forcing data are used in our simulations to represent the oceanographic conditions. The ice–ocean fluxes are represented using the mixed-layer parameterization included in Icepack, which determines the sea surface temperature (SST) and heat exchanges between the sea ice and the ocean based on a fixed mixed-layer depth, sea surface salinity (SSS) and skin friction velocity. Here, we set the SSS to 33 psu (a value coherent with our measured ocean surface temperature of -1.85 °C), the mixed-layer depth to 20 m (default value) and the skin friction velocity to 0.005 m s−1 (the set minimum in Icepack). The SST is prognostic but initialized at the freezing point (as calculated from the liquidus).

The net heat exchange Fbot between the ice and the ocean is given by

(2) F bot = - ρ w c w c h u * ( T w - T f ) ,

where ρw (1026 kg m−3) is the seawater density, cw is the seawater specific heat capacity (4.218 kJ kg−1 K−1), ch (0.006) is a heat transfer coefficient, u* is the ocean friction velocity (0.005 m s−1), and Tw is the water temperature at the SST, and Tf is the ice bottom temperature at the freezing point. Note that when the SST is at the freezing point, Tw=Tf and Fbot=0.

3.1.2 Enthalpy, temperature and salinity profiles

The vertical temperature profiles are computed with boundary conditions set from the surface energy balance described above. The temperature in the snow and ice interior layers is solved to satisfy a prognostic temperature equation, which treats sea ice as a single-phase solid but represents brine via salinity dependencies in the heat conductivity and specific capacity definitions (see Bitz and Lipscomb1999, for details).

The top surface temperature Tsf is determined by the conductive flux needed to balance the net heat flux F0:

(3) F 0 = F ct = K sf 2 ( T sf - T t ) Δ h t ,

where Fct is the top interface conductive flux, Ksf is the conductivity at the air–snow (or air–ice) interface, and Tt and Δht are the internal temperature and layer thickness of the top snow or ice layer. If F0>0, then Tsf is capped to the melting temperature and the remaining imbalance is used to melt snow or ice. The ice bottom temperature Tb is set to the freezing point of surface seawater (Tf).

The internal temperatures in each of the snow or ice layers are governed by the following prognostic equation:

(4) ρ i c i T i t = z ( K i T i z ) - z ( I pen ( z ) ) ,

where ρi is the ice or snow density (917 kg m−3 for sea ice or ρs=330 kg m−3 for snow), ci(T,S) is the specific heat of sea ice or snow, Ti is the internal temperature in the ice or snow layer, Ki is the thermal conductivity based on the bubbly parameterization (Pringle et al.2007), and Ipen(z) is the flux of penetrating solar radiation at depth z according to Beer's law.

The enthalpy q(T,S) of any interface or layer can be retrieved from the solved temperatures as follows:

(5) q ( T , S ) = - ρ c 0 ( T m - T ) + L 0 1 - T m T - c w T m ,

where S is the sea ice bulk salinity (fixed and based on observed vertical salinity profiles; Bitz and Lipscomb1999), c0 (2.106 kJ kg−1 K−1) is the specific heat of fresh ice at 0 °C, Tm(S) is the melting temperature of sea ice as determined by a salinity-dependent liquidus relation, L0 (334 kJ kg−1) is the latent heat of fusion of fresh ice at 0 °C and cw is the specific heat capacity of brine.

3.1.3 Ice congelation

The amount of ice congelation or melt at the ice bottom is given by the imbalance between Fbot and the conductive heat flux adjacent to the ice base (Fcb), according to

(6) - q ( T , S ) h c t = ( F bot - F cb ) ,

where q is the enthalpy at the ice bottom interface as given from Eq. (5). Fcb is defined as

(7) F cb = K b 2 ( T n - T b ) Δ h n ,

where hc represents the thickness of ice formed by congelation at the ice–ocean interface, Kb and Tb (set to Tf) are the conductivity and temperature at the ice–ocean interface, and Tn and Δhn are the temperature and thickness of the lowest ice layer.

3.1.4 Snow-ice formation

The formation of snow ice is represented by converting a fraction of the snow layer to sea ice whenever the hydrostatic balance pushes the snow–ice interface below the water line. This conversion is mass-conserving and instantaneous. The threshold for snow-ice formation is based on Archimedes' law:

(8) h s > ( ρ w - ρ i ) h i ρ s ,

where hs is the snow thickness. The change in snow and ice thicknesses (δhs and δhi) associated with snow-ice formation is written as


where h* is the is the snow thickness in excess of the threshold for snow-ice formation (Eq. 8) before the snow-ice conversion.

3.2 Mushy-layer thermodynamics

3.2.1 Enthalpy, temperature and salinity profiles

In the mushy-layer thermodynamics, sea ice is assumed to be a mixed-phase layer composed of both fresh ice and liquid brine inclusions, with proportions that are determined by prognostic temperature and salinity relations (Feltham et al.2006; Turner et al.2013). The boundary conditions at the top and bottom interface are the same as in the BL99 parameterization, but the internal temperatures in the snow and ice layers are governed by a prognostic equation for enthalpy:

(11) q t = z ( K i T i z ) + w q br z - z ( I pen ( z ) ) ,

where qbr is the enthalpy of the brine and w is the Darcy velocity of the brine. The enthalpy q is defined in terms of the brine fraction and temperature as

(12) q = ϕ q br + ( 1 - ϕ ) q i = ϕ ρ w c w T + ( 1 - ϕ ) ( ρ i c i T - ρ i L 0 ) ,

where qi is the enthalpy of fresh ice and ϕ is the liquid fraction defined as

(13) ϕ = S S br ,

where Sbr is the salinity of the brine as defined by an observation-based liquidus relation (Turner et al.2013). Together, Eqs. (11) and (12) differ from the BL99 thermodynamics only from the additional heat advection from brine flow and the mixed-phase enthalpy definition.

The prognostic salinity equation in each ice layer includes dependencies on brine processes such as gravity drainage and melt pond flushing (Notz and Worster2009; Turner et al.2013). The brine drainage component (Turner et al.2013) is written as

(14) S t + v z S br z = S t | slow ,

where vz is the vertical velocity of the ocean water percolating upward through the ice layer in response to the brine drainage (rapid drainage mode). The right-hand side represents a slow mode of brine drainage that varies with the surface temperature according to the following equation:

(15) S t | slow = - ω ( S - ϕ c S br ) ( T b - T sf ) h i if  T b > T sf , 0 otherwise,

where ω is a tuning parameter set by the user (-5.0×10-9 m s−1 is the default value) determining the strength of the slow drainage and ϕc is a critical liquid fraction for the slow drainage to occur, also set by the user (0.05 is the default value in Icepack). More details can be found in Turner et al. (2013).

3.2.2 Standard mushy-layer congelation

In mushy-layer physics, there is no sharp interface between solid ice and ocean water, but rather a downward transition within the mushy medium towards a 100 % liquid fraction. As such, ice congelation is not made by forming a layer of solid ice but by moving the ice–ocean boundary downward at a rate defined by the conductive heat flux imbalance and then by integrating the corresponding amount of seawater in the bottom ice layer. The solidification of the seawater is thus only treated in subsequent time steps when implicitly solving for the temperature profiles, during which the liquid fraction is adjusted to satisfy the liquidus relation.

Specifically, the congelation rate (i.e., the migration of the ice–ocean boundary) is first defined based on the energy needed to form a mushy layer with an initial liquid fraction ϕinit=0.85 (default value) as

(16) h c t = F bot - F cb L ρ i ( 1 - ϕ init ) .

Then, the enthalpy and salinity of the lowest ice layer are updated to account for the added seawater spanned by the moving boundary, using


where the subscript n refers to the lowest ice layer, qw is the enthalpy of seawater at the freezing point and Sc is the bulk salinity of the integrated seawater (i.e., SSS).

Note, however, that Eq. (17) does not include the phase change of the solid ice fraction (1−ϕinit) assumed in Eq. (16) and thus does not fully account for conductive heat imbalance FbotFcb. This leads to a leftover being taken from the ocean, resulting either in a cooling of the SST, or, if the SST is already at the freezing point, in frazil formation. Defining Focn as the leftover heat flux imbalance taken from the ocean after congelation, the rate of frazil formation in the latter case is defined as

(19) h f t = F ocn q f ,

where hf is the ice thickness from frazil formation and qf is the enthalpy of the frazil as defined from Eq. (12) using a liquid fraction of 0.75 (smaller than ϕinit for congelation) and temperature corresponding to the liquidus for a brine salinity of Sbr=SSS-3.

The total growth at the ice base in mushy-layer simulations is thus obtained by combining the congelation growth/melt with the frazil (i.e., hc+f=hc+hf; see Appendix A for more details).

3.2.3 Modified mushy-layer congelation

To improve our mushy simulation results, we propose a modification to the congelation parameterization that better accounts for the heat flux imbalance at the ice base (i.e., reducing the associated ocean cooling or frazil formation). Our modifications are twofold. First, we assume that the solid ice formation is simultaneous with the moving boundary, such that the congelation mush layer with liquid fraction ϕinit is explicitly incorporated into the lowest ice layer (instead of only seawater in the standard parameterization described above). Second, we define the congelation rate as a function of the energy needed to decrease the enthalpy of the original seawater to that of the new congelation mush. Together, these modifications ensure that the enthalpy of the added congelation layer corresponds with the conductive heat imbalance at the ice–ocean interface. More details can be found in Appendix B.

Specifically, the mushy-layer congelation rate (i.e., the migration of the ice–ocean boundary) is now defined as

(20) h c t = F bot - F cb q w - q m ,

where qm is the enthalpy of the integrated congelation mush layer as defined by Eq. (12), with a liquid fraction ϕinit and at the freezing point temperature. The enthalpy of the lowest ice layer is updated by integrating the congelation mushy congelation layer spanned by the moving boundary:

(21) ( h n q n ) t = q m h c t .

The salinity update in our scheme remains given by Eq. (18) but using Sc=ϕinitSbr, where Sbr= SSS.

3.2.4 Snow-ice formation

In the mushy-layer scheme, the snow-ice formation remains based on the hydrostatic balance (Eq. 8), but the conversion of snow to ice is no longer mass-conserving (in stand-alone simulations). Instead, snow flooding is parameterized by adding seawater to a fraction of the snow layer, thus assuming that seawater either is advected laterally or percolates through the ice layer. The change in snow and ice thickness associated with snow-ice formation is given by

(22) δ h i = - δ h s = m fb ρ w + ρ s - ρ snice ,

where mfb (=hiρi+hsρs-hiρw) is the combined mass of snow and ice in excess of the hydrostatic equilibrium prior to the snow-ice formation and ρsnice is the density of the newly formed snow ice. The snow-ice density and liquid fraction ϕsnice are defined by assuming that seawater has filled the porosity of the snow layer:


In this analysis, we also test the use of additional criteria for snow flooding. In these specific simulations, the flooding onset either is allowed only after the observed flooding onset date (i.e., set manually, as in Duarte et al.2020) or if the ice layers are sufficiently permeable (i.e., if the smallest liquid fraction in all ice layers is larger than a liquid fraction threshold ϕmin). Otherwise, negative freeboards can develop. To avoid large and sudden snow flooding once the criteria are met, we include a simple linear dependence of the flooding rate on the negative freeboard:

(25) δ h i = - δ h s = γ m fb ρ w + ρ s - ρ snice ,

where γ is a free parameter set here to 0.01 to match the observations.

4 Methods

4.1 Snow depth and ice thickness retrieval

The in situ snow depth, ice thickness, congelation growth and snow-ice formation are determined using a new automated surface retrieval algorithm. Our algorithm is similar to that of Liao et al. (2019) and Cheng et al. (2020) with a few adaptations that aim to reduce its sensitivity to large diurnal cycles and to improve its performance in near-isothermal conditions. As in Cheng et al. (2020), it detects the material interfaces based on the vertical gradients in the temperature profiles and is built to detect snow flooding (i.e., upward ice growth at the snow–ice interface), which was suspected at our deployment sites. Similar vertical-gradient-based algorithms were recently shown to be most appropriate compared to other methods for the automated retrieval of ice thickness from IMB data (Gough et al.2012; Richter et al.2023).

The ice thickness and snow depths are determined from the position of three material interfaces on the SIMBA temperature profiles: the top of the snow layer (the air–snow interface, shown as Za−s in Fig. 2), the snow–ice interface (Zs−i) and the bottom ice–ocean interface (Zi−o). Since a segment of the thermistor string is laid flat (horizontal) on the ice surface at deployment, the algorithm also needs to identify the first (Zice0 in Fig. 2) and last (Zp) sensors of this “thermistor plateau”, which becomes embedded in the ice after flooding events (see Fig. 2b for the flooded ice case). These locations are first detected for each individual profile (at a 6 h interval) and then smoothed using a 24 h running mean to remove any sensitivity to the diurnal cycles.

The ice thickness hi (including snow ice), snow depth hs, and snow-ice thickness hsi are calculated from the five identified positions, according to


The changes in ice thickness can thus be associated with an upward displacement of the snow–ice interface (defining the snow-ice contribution to the mass balance) or a downward displacement of the ice bottom interface (defining the congelation contribution to the mass balance).

The surface retrieval algorithm is based on the following assumptions:

  1. The temperature profiles are piecewise linear.

  2. The snow–ice interface does not move downward along the thermistor string (i.e., no vertical slip between the buoys and the ice and no surface melting).

  3. The minimum temperature along the thermistor string is located above the snow layer.

  4. The vertical profiles are isothermal in the ocean.

These assumptions are similar to those from Liao et al. (2019), Zuo et al. (2018) and Cheng et al. (2020), and they relate to the dependency of the algorithm on the difference in heat conductivity (i.e., vertical temperature gradient) in the snow and ice layers. Heat-conductivity-based surface retrieval algorithms are thus, by construction, not suited for near-isothermal conditions (e.g., during thaw), in which case other observations (e.g., from sonar data or the SIMBA heat cycles) are needed to determine the ice mass balance. In all, the algorithm described below is similar in principle to that of Cheng et al. (2020) and only differs in the detection criteria for each interface.

4.1.1 Temperature gradient and curvature

The vertical temperature gradient β and curvature γ are first calculated at each sensor location and for the entire data record using a centered finite-difference scheme. The vertical temperature gradient at the kth sensor location is defined as

(29) β k = T k z T k + 1 - T k - 1 2 Δ z ,

where Tk represents the temperature reading of the kth sensor and Δz is the spacing between two sensors (here 2 cm). The curvature at point k is defined as

(30) γ k = 2 T k z 2 T k + 1 - 2 T k + T k - 1 Δ z 2 .

4.1.2 Initial ice surface and thermistor plateau

For each buoy, the thermistor plateau is set at deployment and remains fixed over the entire record. The initial ice surface Zice0 (with temperature Tice0) and lower end of the thermistor plateau serve as reference points for the algorithm.

The position Zice0 is identified by the minimum curvature (min(γk)) below the maximum vertical temperature gradient in the profiles (assumed to be inside the snow layer, Fig. 2). The other end of the thermistor plateau Zp is identified by the closest local maxima in curvature below Zice0. To remove sensitivity to sporadic variations in the detected interfaces (±2 cm), the reference locations are defined as the statistical mode of Zice0 and Zp over the first 7 d of records.

4.1.3 Ice–ocean interface

For each profile, the position of the ice–ocean interface is determined using a minimization approach to find the sensor location best matching the corresponding change in the vertical temperature slope. That is, to each tentative ice bottom position Zl, where l represents a specific sensor location k=l close to the expected ice bottom, we assign a theoretical piecewise linear vertical temperature profile, defined as

(31) T k th = T c + ( z k - Z c ) β ice if  Z c > z k > Z l , T w if  z k < Z l ,

where Tkth is the theoretical temperature at sensor location zk, Tc is the temperature observed at a position Zc in the ice interior (here defined as TcTw+r(Tice0-Tw), where r=1/3 is an arbitrary scaling factor), βice is an ice temperature gradient initial guess and Tw is the observed ocean temperature. The initial guess βice is defined as

(32) β ice = T w - T c Z l - Z c .

The position of the ice bottom interface Zi−o (and associated βice) is then defined from the position Zl with the theoretical profile that minimizes the following error function:

(33) err = k = l - 10 l + 10 T k th - T k obs 2 ,

where Tkobs is the observed temperature at sensor position k.

Note that this detection method differs significantly from the temperature selection method of Liao et al. (2019) and Cheng et al. (2020), with the benefit of not depending on the sensor type and precision.

4.1.4 Air–snow interface

The air–snow interface position Za−s is found by identifying the maximum vertical temperature curvature γk below the sensor with the coldest temperature reading (assumed to be in the air) and above the initial ice surface Zice0. The temperature gradient directly below Za−s must also be smaller than a threshold for snow detection, set to 0.1 °C cm−1. Note that this threshold is smaller than in Liao et al. (2019) but is only used to discriminate curvatures associated with noise in the data. The temperature gradient in the snow layer is then defined as

(34) β snow = T ice 0 - T a - s Z ice 0 - Z a - s ,

where Ta−s is the temperature reading at Za−s.

4.1.5 Snow–ice interface

The presence of snow ice above the initial ice surface is detected by comparing the temperature gradient directly above the initial ice surface Zice0 with βsnow and βice. That is, sensors above the original ice surface are associated with snow ice if the local temperature gradient satisfies

(35) β k < β ice + r si ( β snow - β ice ) ,

where rsi (0.2) is a ratio between 0 and 1. If such a gradient is present above Zice0, then the new ice surface position (Zs−i) is updated to the lowest point where βk<βsi but only if detected for at least 4 consecutive days.

Note that while arbitrary, the ratio rsi for snow-ice detection ensures that the snow-ice conductivity is closer to that of sea ice, while also filtering fluctuations due to changing temperature conditions. The snow-ice detection is the only component of the algorithm that depends on the other detected interfaces.

4.2 Freeboard computation

The ice freeboard hfb is the elevation of the snow–ice interface above the water line. A negative freeboard value indicates that the snow–ice interface is below the water line with the ice in hydrostatic imbalance. In both the observations and simulations, we compute the freeboard based on the hydrostatic balance using the same material densities defined in Icepack (see Sect. 3):

(36) h fb = h i - ρ s h s + ρ i h i ρ w .

Based on the propagation of uncertainty and assuming an error of 2 cm for the snow and ice thicknesses and of 33 kg m−3 for the snow density (King et al.2020), these freeboard estimates have a precision of ∼1.0 cm.

Table 1Parameters and performance in mean integrated error (MIE) for all SIMBA1 simulations.

Download Print Version | Download XLSX

4.3 Experiment setup

Multiple Icepack simulations are run with the BL99 physics or the mushy-layer physics to reproduce each of the SIMBA observations. All simulations are initialized using the ice thickness, snow depth and internal ice temperature (at locations corresponding to the center of the snow and ice layers) recorded by the buoys. The initialization values are taken on 1 March, a few days after the SIMBA deployment to ensure that the deployment holes are completely refrozen. The simulations are run with seven ice layers, one snow layer and a time step of 1 h (outputs only every 3 h) from 1 March until well past the buoy recovery date. Results are only shown for the period corresponding with observation records.

The parameter specifications for each simulations are listed in Tables 1 and 2 for the SIMBA1 and SIMBA2 experiments, respectively. Control simulations are first run using the default Icepack BL99 physics and mushy-layer parameterizations. A series of experiments are then conducted to investigate the influence of individual mushy physics components on the model performance by (1) removing the snow-ice parameterization, (2) adding a minimum liquid fraction (ϕmin) criteria for snow flooding, (3) varying the strength of the brine drainage and (4) using the modified congelation scheme with varying initial liquid fraction ϕinit.

Table 2Parameters and performance in mean integrated error (MIE) for all SIMBA2 simulations.

Download Print Version | Download XLSX

4.4 Model evaluation

The performance of each simulation is quantified using the mean integrated error (MIE) of the ice thickness, snow depth, cumulative congelation and snow-ice formation. For each variable, the MIE is calculated first by linearly interpolating the SIMBA and simulation data into an hourly time series. The MIE is then defined as

(37) MIE = τ = 1 n ( X sim τ - X obs τ ) n ,

where n is the number of valid data points in the time series and Xsimτ and Xobsτ are the simulated and observed variable values at time τ in the interpolated time series.

Figure 3Time series of (a) air temperature from the GDPS (black) and recorded by the SIMBA buoys (SIMBA1 in blue, SIMBA2 in green), (b) precipitations from the Nain ECCC weather station (black) and from the GDPS (blue), (c) recorded temperatures along the SIMBA1 thermistor string with the detected material interfaces (air–snow interface in blue, ice interfaces in black and thermistor string plateau in red); (d) same as (c) but for the SIMBA2 buoy. The purple arrow points to the warming at the snow–ice interface, indicating flooding.


5 Results

5.1 In situ landfast ice thermodynamics

5.1.1 Observed temperature and weather conditions

The late-winter conditions along the Labrador coast are characterized by increasingly large diurnal cycles in air temperature, with longer (synoptic-scale) events of colder or warmer weather (Fig. 3a). The 2 m air temperatures calculated from the GDPS data correspond well with the air temperatures recorded in situ but are generally colder (MIE of −0.78 and −0.58 °C compared to the SIMBA1 and SIMBA2 records, respectively). These biases are mostly associated with differences in the short-term temperature peaks; the buoys record larger maxima in air temperatures than represented in the GDPS data.

Several precipitation events occurred during the observational periods. The snow precipitation events from the GDPS data correspond well with the precipitations recorded at a nearby weather station (Nain Airport; Fig. 3b). The precipitation phases were not documented in the airport records, and all events were snowfalls in the GDPS data. In particular, two events with heavy snowfalls are recorded on 9–11 March and on 6–10 April, which also correspond to periods of warmer weather during which temperatures slightly exceeded the freezing point.

The vertical temperatures recorded along the two SIMBA thermistor strings are coherent with these patterns (Fig. 3c–d). Short-term variations in air temperature are rapidly damped in the snow layer even though heat from longer periods of warm weather reaches, and has a larger impact on, the ice interior. The downward propagation of the surface heat is often followed by a slower cooling once colder conditions return. Despite the similar air temperature patterns, the SIMBA2 recorded significantly warmer ice temperatures than SIMBA1, with sharp warming events (see purple arrow in Fig. 3d) that suggest a snow-flooding onset (Provost et al.2017).

Figure 4Rates of temperature changes (a, b), vertical temperature gradients (c, d) and change in temperature recorded after 2 min of heating during the daily heating cycles (e, f) at each sensor as measured for SIMBA1 (a, c, e) and SIMBA2 (b, d, f). Colored lines indicate the detected material interfaces (air–snow interface in blue, ice interfaces in black and thermistor string plateau in red). The purple arrow points to the warming at the snow–ice interface, indicating flooding.


5.1.2 Surface retrieval algorithm validation

The surface retrieval algorithm is able to identify the snow and ice interfaces in most of the records (Fig. 3c–d). The algorithm fails during the two warm spells when negligible vertical temperature gradients or temperature inversions are present within the snow and ice layers (i.e., the piecewise linear assumption does not hold). The surface retrieval algorithm is also generally not successful during the melt season (beyond 16 April) for the same reason, except on occasional colder days.

While we do not have independent data to validate the retrieved snow and ice thicknesses, we find that the selected interfaces are coherent with the interfaces detectable by visual inspection in the temperature profiles (Fig. 3c–d). The detected snow interfaces correspond well with the layer within which most of the variability associated with diurnal cycles or synoptic systems is damped (Fig. 4a–b) and where large vertical temperature gradients are present (Fig. 4c–d). The algorithm detects an upward migration of the snow–ice interface (i.e., snow flooding; see the upward displacement of the top black line, representing the snow–ice interface, above its original position) that also corresponds well with the warm temperatures recorded above the initial ice surface. In particular, the onset of flooding at the SIMBA2 site (on 26 March) coincides with a sudden warming event observed at the snow–ice interface, propagating upward in the snow layer despite a cooling in the surface air temperature above (see profiles at the purple arrow in Figs. 3d and 4b). This signal is expected when the snow flooding is caused by upward percolation or lateral advection of seawater (Provost et al.2017), since the warm seawater increases the snow–ice interface temperature, and the heat later diffuses upward. In contrast, flooding by liquid precipitation or snow melting would show the entire snow layer at the freezing point. This could be the case at the SIMBA1 site, where flooding is only detected late in the observational record (on 25 April) when surface air temperatures above freezing are regularly present.

The top and bottom ice interfaces show good agreement with those seen in the recorded warming of sensors during the SIMBA heating cycles (Fig. 4e–f). The detected snow layers are also coherent with the thermistor string sections measuring the largest heating (smallest conductivity), although this is more difficult to assess with certainty due to the large variations within this layer, likely resulting from a vertically varying snow density. Note that the SIMBA2 heat cycle records (Fig. 4f) present a rather smooth vertical gradient over 2–4 cm within the thermistor plateau, supposedly sitting on the snow–ice interface. We speculate that this is due to the thermistor plateau not being exactly horizontal on the ice surface, indicating a (∼1–2 cm) thickness uncertainty related to this deployment method for a marking of the initial ice surface. This positional uncertainty remains for the entire record but is no longer visible once the thermistor plateau is flooded.

Figure 5(a) Ice (blue lines), snow (green lines) and freeboard (orange lines) thicknesses from the SIMBA data. (b) Contribution of snow ice (blue lines) and congelation ice (orange lines) to the ice mass balance inferred from the SIMBA data.


5.1.3 In situ landfast ice mass balance

The SIMBA observations show large snow depths (∼20–40 cm) over relatively thin ice (∼75–100 cm) from the beginning of the records, and the measured freeboard occasionally dips to negative values (Fig. 5a). Both sites present significant snow depth increases during each warm event with a subsequent reduction likely resulting from snow compaction and redistribution by winds. The snow depths are generally larger at the SIMBA2 site (by ∼5–10 cm), with a large but short-lived maximum of 50 cm likely resulting from snow accretion and subsequent removal by the winds around the buoy.

The local ice mass balance at the two sites is largely influenced by the snow layer thickness and its insulating effect on the sea ice below. The thinner snow cover at the SIMBA1 site results in colder internal ice temperatures, larger congelation rates at the ice base and less snow flooding (Fig. 5b). With an initial ice thickness and snow depth of 80 and 26 cm, respectively (on 1 March), the SIMBA1 freeboard reaches negative values after each snowfall event: −1.8 cm on 13 March and −1.6 cm on 14 April. Snow flooding is only detected from 25 April onward. The ice thickness reached its maximum (100 cm) on 1 May, for a total ice growth of 20 cm, from which 16 cm is associated with (downward) congelation at the ice–ocean interface and 4 cm is associated with (upward) snow-ice formation. In comparison, the SIMBA2 buoy initially recorded a 30 cm snow depth and 76 cm ice thickness (on 1 March), already corresponding to a negative freeboard (−1.6 cm). Snowfalls during the first warm event brings the freeboard to a minimum of −6.4 cm on 15 March. Snow flooding is detected from 25 March onward, coinciding with a large (∼10 cm) reduction in the snow depth. By 6 April, the ice thickness reached a maximum of 97 cm for a total ice growth of 21 cm, 13 cm of which is attributed to snow-ice formation and 8 cm to congelation.

Figure 6Simulated internal temperatures interpolated into 2 cm intervals from the BL99 (a, b) and mushy simulations (c, d), initialized from the SIMBA1 (a, c) and SIMBA2 (b, d) data. Solid lines indicate the simulated material interfaces (air–snow interface in blue, ice interfaces in black).


Figure 7Ice mass balance in the BL99 simulations against the SIMBA1 (a, c) and SIMBA2 (b, d) observations. Top panels (a, b) represent ice thickness (blue lines), snow depth (green lines) and freeboard (yellow lines) values, with the observations in black. Bottom panels (c, d) represent cumulative ice growth from ice bottom (yellow lines) and snow-ice formation (blue lines), with the observations in black. Thin lines indicate results from the BL99 simulation ran without using the snow-ice parameterization.


5.2 BL99 simulations

The BL99 thermodynamics represents generally well the observed internal temperature profiles but with a larger downward heat conduction in the ice interior during periods of warm weather compared to the observations (Fig. 6). The simulated snow thicknesses present large discrepancies with observations (MIE of +1.88 and −3.07 cm for SIMBA1 and SIMBA2, respectively), with a tendency to remain mostly constant except for the occasional increases associated with precipitation events. This is mostly due to the simple snow model not accounting for snow compaction and redistribution. The simulated ice thickness is in general accord with observations (Fig. 7a–b) with a small positive bias +0.25 cm (MIE) for SIMBA1 and +2.13 for SIMBA2 (see Tables 1 and 2). Despite the positive MIE values, the ice thickness at the time of the observed maximum is smaller than the observed values, at 97.4 and 90.8 cm for the SIMBA1 and SIMBA2 simulations, respectively (−2.6 and −6.2 cm underestimations). Most of the ice growth is attributed to ice congelation at the ice bottom (14.9 and 10.5 cm), showing slowly decreasing ice growth rates from ∼0.3 cm d−1 to near zero in May. The volume of snow ice is largely underestimated at 2.4 and 4.3 cm (−1.6 and −8.7 cm underestimations) despite the fact that conditions for snow-ice formation are met from the very start of the simulation (Fig. 7c–d).

The ice thickness and snow depth discrepancies are partly attributed to the misrepresentation of snow-ice formation, specifically to the hydrostatic-balance-based snow-flooding onset: the initialized snow depths are sufficient enough to depress the ice surface near (or already below in the SIMBA2 case) the water line, with any subsequent snow precipitation leading to a portion of the snow cover being immediately transformed into snow ice (see orange lines representing the freeboard values in Fig. 7a–b and blue lines representing the snow-ice volumes in Fig. 7c–d). This leads to the ice thickness temporarily exceeding the observations early in the simulations up to the observed snow-flooding onset, after which the thickness bias turns negative (in the SIMBA2 case) due to the small snow-ice volume.

Figure 8Ice mass balance in the mushy-layer simulations against the SIMBA1 (a, c) and SIMBA2 (b, d) observations. Top panels (a, b) represent ice thickness (blue lines), snow depth (green lines), and freeboard (yellow lines) values, with the observations in black. Bottom panels (c, d) represent cumulative ice growth from ice bottom (yellow lines, including frazil) and snow-ice formation (blue lines), with the observations in black. Thin lines indicate results from the mushy simulations run without using the snow-ice parameterization.


5.3 Mushy simulations

Compared to the BL99 simulations, the mushy-layer physics produces warmer sea ice temperatures (see Fig. 6c–d) and faster ice growth at both interfaces (i.e., upward snow-ice formation and downward bottom ice growth; see Fig. 8). These differences are present despite the fact that the simulated snow depths are very similar in both simulations.

The ice thickness in the mushy simulations reached 103.8 and 97.5 cm for SIMBA1 and SIMBA2, respectively, at the time of the observed maximum, corresponding to 3.8 and 0.5 cm overestimations. The ice growth presents larger variations than observations due to a combination of spurious snow-ice formation and variable basal ice growth. The spurious snow-ice formation is similar to but larger than in the BL99 simulations, yielding large ice thickness discrepancies during the period with observed negative freeboards. The total volume of snow ice is, however, closer to observations, with 8.0 and 13.0 cm for the SIMBA1 and SIMBA2 simulations, respectively (+4.0 and 0.0 cm deviations from observations). The basal growth variability could be considered an asset when compared with the slowly varying congelation rates in the BL99 simulations, but it is largely overestimated and effectively degrades the model performance. In particular, the simulated basal growth features periods of weak basal melt that are not coherent with observations (Fig. 8c–d). Furthermore, a quarter of the basal growth is attributed to the frazil formation during periods of rapid congelation (see the Appendix A). This frazil formation occurs despite the landfast ice conditions with a 100 % sea ice concentration and uniform thickness and is related to the treatment of the ice–ocean boundary in the mushy-layer congelation parameterization (see Sect. 3.2.2 and Appendix for details).

All control simulations (BL99 and mushy) thus present discrepancies early in the simulations due to the hydrostatic-balance criteria not accounting for negative freeboards. This difficulty lead Duarte et al. (2020) to manually activate/deactivate the snow-ice parameterization (i.e., by adding/removing the condition of the hydrostatic equilibrium according to the observations) to adequately reproduce in situ conditions. In our experiment, deactivating the snow-ice parameterization in all simulations effectively allows the snow depth to increase during precipitation events (see thin green lines in Figs. 7 and 8) and reduces the ice thickness discrepancy up to the flooding onset (thin blue lines in Figs. 7 and 8). This, however, leads to an underestimation of the ice thickness by the end of the simulations due to the missing snow-ice contribution in the ice mass balance.

The removal of snow-ice formation causes different responses in the BL99 and mushy-layer thermodynamics. Using the BL99 physics, simulations without snow ice show smaller congelation rates due to the increased insulation from the larger snow depths. Using the mushy-layer physics, simulations without snow ice show larger congelation rates despite the larger snow depths, since the ice interior is colder without the influx of warm ocean water associated with flooding, resulting in larger conductive heat fluxes at the ice base.

Figure 9Time series of the simulated vertical temperature profiles interpolated in 2 cm intervals to reproduce the SIMBA2 records using the mushy-layer physics with different criteria for snow flooding. Thick lines indicate the material interfaces (air–snow in blue, ice interfaces in black and the initial snow–ice interface in dashed red). (a) Without snow flooding, (b) using ϕ=0.005 as a snow-flooding onset criteria and (c) manually setting the snow-flooding onset on 26 March to match the observations.


We note that while the mushy-layer simulations quantitatively represent a degradation of the model performance (see larger MIE values in Tables 1 and 2), this is largely due to the snow-flooding onset discrepancy combined with the wider-ranging effects of the flooding on the ice thickness growth, interior ice salinity and temperatures. These effects, however, are physically meaningful and correspond well with previously recorded thermodynamics of snow flooding (see Provost et al.2017, for instance). This suggests that the model performance could be largely improved by a simple tuning of the snow-flooding onset and rates. For instance, we find that adding a simple minimum porosity criterion (ϕmin=0.005) to the snow-ice parameterization and setting the flooding rate inversely proportional to hfb largely improves the SIMBA2 simulations by delaying the snow-ice formation by several days (Fig. 9b). The model in particular presents very small MIE values for snow-ice formation when the flooding onset is set manually to the observed date (Fig. 9c and Table 2).

Figure 10Time series of the (a) temperature, (b) bulk salinity and (c) brine salinity in the upper ice layer in mushy simulations with different criteria for snow-ice formation (blue lines represent no flooding).


5.4 Basal ice temperature, brine salinity and congelation

The inclusion of prognostic salinity and brine parameterizations in the mushy-layer physics adds model sensitivities relating to the liquidus relationship. That is, for a given salinity, the liquidus relation interconnects changes in temperature with changes in brine salinity and liquid fraction. As such, updating the brine salinity in explicit parameterizations, such as the snow-flooding parameterizations or ice congelation parameterizations, later affects the layer temperature solved implicitly in subsequent time steps. For instance, the seawater added in the upper ice layer in the snow-ice parameterization increases the layer bulk salinity but also dilutes the brine salinity towards SSS values. This effectively warms the layer according to the liquidus balance (Fig. 10). The layer temperature then slowly returns to colder values as the brine pockets refreeze, concentrating the brine salinity to its original value (see curves converging back to values from the negative freeboard simulations in Fig. 10a).

Figure 11Time series of the temperature (a, b), bulk salinity (c, d), brine salinity (e, f) and brine desalination rate from the slow-brine-drainage parameterization (g, h, represented in psu d−1), all in the lowest ice layer. Left panels showy mushy-layer simulations using different brine drainage strength parameters ω (using the standard congelation). Right panels show mushy-layer simulations using the different initial congelation liquid fraction ϕinit and using the modified mushy-congelation scheme.


Figure 12Ice mass balance in mushy-layer simulations tuned to best represent the observations from the SIMBA1 (a, c) and SIMBA2 (b, d) buoys (observations are in black). The snow-flooding onset is set manually according to the observed flooding onset dates, and the simulations use the modified congelation scheme with ϕinit=0.45. Top panels (a, b) show ice thickness (blue lines), snow depth (green lines) and freeboard (yellow lines) values. Bottom panels (c, d) show cumulative ice growth from ice bottom (yellow lines, including frazil) and snow-ice formation (blue lines).


Similarly, the alternating periods of sea ice congelation and melt in the standard mushy-layer simulations are attributed to a similar brine–temperature feedback in the lowest ice layer: any process reducing the brine salinity yields an increase in the layer temperature Tn. This reduces the conductive flux at the ice base (see Eq. 7, with Tf constant at the freezing point) and thus the available energy for congelation. Specifically, there are two explicit parameterizations inducing brine salinity changes at the ice base in the mushy-layer thermodynamics: the brine drainage parameterizations (reducing the brine salinity) and the ice congelation parameterizations (diluting the brine towards the SSS). These parameterizations together bring the brine salinity close to SSS values early in the simulations (see blue curve in Fig. 11c for the control simulation). Later, brine drainage under cold weather further dilutes the brine to values below the SSS (and thus, Tn>Tf), causing a reversal of the conductive flux and sea ice melt. This pattern can be suppressed by reducing the strength of the brine drainage (reducing the parameter ω, Fig. 11, left panels), although it also consequently yields too large congelation rates.

The basal ice growth can be improved by modifying the congelation parameterization to reduce the associated salinity increase (Fig. 11b, d, f and h). To do so, we repeat the experiments using a modified congelation scheme in which a mush layer with a liquid fraction ϕinit and Sbr= SSS is incorporated in the lowest ice layer during congelation (see Sect. 3.2.3 and Appendix B). Using this scheme, reducing the liquid fraction of congelation ice (ϕinit) results in smaller congelation rates and salinity in the lowest ice layer. This in turn reduces the strength of the brine drainage (a lower salinity in Eq. 15), diminishing the variations in congelation while bringing the congelation rates closer to observations.

6 Discussion and conclusions

In this study, the thermodynamic growth of landfast ice in the vicinity of Nain (Labrador) is investigated from two Scottish Association for Marine Science (SAMS) snow ice mass balance apparatus (SIMBA) buoys deployed in winter 2017. The observed thermodynamics are reproduced using Icepack (v1.1.0), which is the column thermodynamics package of the Community Ice CodE (CICE) version 6, with two different physical schemes: the Bitz and Lipscomb (1999) physics that represent the thermodynamics currently used in the Environment and Climate Change Canada (ECCC) ice–ocean forecasting systems and the mushy-layer thermodynamics (Feltham et al.2006; Notz and Worster2009; Turner et al.2013) that include new physics available in CICE6. The performance of Icepack in reproducing the IMB observations is assessed with particular attention to the improvements associated with the use of the mushy-layer physics. The contributions of this paper include a new automated surface retrieval algorithm to infer the ice and snow thicknesses from thermistor string records, a modified mushy-layer congelation scheme less conducive to frazil formation and modifications to the snow-flooding parameterization to allow for negative freeboards and slow snow flooding rates.

The in situ observations presented in this analysis are in line with a number of negative freeboard measurements reported in recent years in the Arctic (Rösel et al.2018; Provost et al.2017; Duarte et al.2020), which are likely to become more frequent as the sea ice thins and precipitation increases in the transition to a seasonal ice cover (Merkouriadi et al.2020). Snow flooding remains nonetheless relatively infrequent: our in situ snow flooding observations were associated with anomalous 2017 snow conditions that have not yet reoccurred in subsequent (2018–2023) landfast ice observation campaigns. The frequency at which snow flooding contributes to the ice mass balance in landfast ice areas, in Nain but also more widely along the Canadian Arctic, remains to be determined. Note, however, that as snow-ice formation occurs more easily over thin ice (Granskog et al.2017), it is likely contributing to the ice growth early in the season and in new leads. This could be better assessed with IMB buoys deployed in open water prior to the freeze-up. Such a deployment was attempted in 2022 in Nain, but buoy icing, floe drifting and wave battering prevented the measurement of a continuous time series during the freeze-up period.

The large discrepancies between the observed and simulated snow-flooding onset in the analysis join the results of Duarte et al. (2020) in demonstrating that the use of the hydrostatic balance alone is insufficient to define snow flooding and to capture the more complex processes observed in situ (Eicken et al.1995; Maksym and Jeffries2000; Provost et al.2017). Our results show that while this conclusion also applies to the BL99 parameterization, the snow flooding exerts a much wider-ranging thermodynamic response under the mushy-layer physics as the flooding increases the temperature, salinity and liquid fraction in the upper ice layers. It better represents the observed thermodynamics and is an improvement compared to the BL99 physics, as indicated by the smaller MIE values when the flooding onset is corrected according to the observations (see Table 2). One advantage of the mushy-layer physics is that it contains the necessary ingredients to improve the snow-flooding parameterization with additional porosity conditions for the percolation of seawater through the brine channels.

In our analysis, no porosity criterion was found to reproduce the observed snow-flooding onset date. This could indicate the influence of nearby sea ice dynamics, although in our case, the deployed IMBs were located in a well-sheltered landfast channel dozens of kilometers away from the landfast ice edge. Moreover, the slow rate of snow-ice formation corresponds well with percolation through the porous sea ice medium (i.e., as opposed to the sudden flooding expected when floodwater is advected laterally from neighboring deformation sites; Provost et al.2017). One difficulty in reproducing the snow-flooding onset with porosity criteria is that they do not account for a percolation associated with the larger-scale porosity (e.g., from thermal cracking) unrelated to the smaller-scale mushy-layer characteristics. At the kilometer scale of most dynamic sea ice models, the volume of snow ice will likely not be uniform over a grid-cell area. This is made evident in our results by the different in situ flooding onset recorded by our two neighboring SIMBAs. Most likely, the snow-ice volume will be spatially distributed according to the ability of the floodwater to penetrate the snow layer, ultimately depending on the ice topography (ice thickness distribution), local snow conditions and the ice heterogeneity (i.e., the presence and average distance between cracks). The snow-ice volume at this scale would thus likely be better represented by a subgrid parameterization relating the snow conversion to a spatial probability for water penetration.

Our results further demonstrate that the mushy-layer physics leads to a much larger salinity and temperature variability at the ice bottom, with significant sensitivity to new free parameters (e.g., ω, ϕinit). This highly impacts the simulated ice congelation rates and, using the default Icepack parameters, yields a degradation of the model performance despite the improved representation of brine processes. This performance is, however, mostly associated with the treatment of the brine salinity in the explicit congelation parameterization producing overinflated congelation rates, erroneous melt and significant frazil formation. The frazil formation in particular is not expected in our sheltered landfast context, but its overrepresentation is coherent with previous studies reporting large frazil volumes in pan-Arctic or Antarctic simulations using the mushy-layer thermodynamics (Turner and Hunke2015; Bailey et al.2020; DuVivier et al.2021).

The mushy-layer thermodynamics can, nonetheless, outperform the BL99 simulations using a simple modification to the congelation parameterization with a tuning of the initial congelation liquid fraction. Note, however, that the modified congelation is not salt-conserving (i.e., similar to the frazil formation) and should be treated accordingly in the context of coupled simulations. The best model performance was obtained when the mushy-layer physics was used with the modified congelation parameterization, a reduced initial congelation liquid fraction and a manual snow-flooding onset (Fig. 12). We note, however, that this represents significant tuning towards our SIMBA observations, which are not representative of typical high-Arctic conditions. As such, this tuning exercise is not meant to determine specific parameter values to be used in Icepack. It nonetheless demonstrates the need to better constrain the mushy-layer parameters, which could be made in future work with larger sets of in situ observations including salinity measurements.

Finally, we note that the increased sensitivity to physical processes in the mushy-layer thermodynamics is likely to positively affect the landfast ice dynamics. For instance, the larger congelation rates simulated under conditions of colder air may allow for faster sea ice consolidation (increasing the effective ice strength) and ease the formation of ice arches in narrow passages. The impact of snow flooding, precipitation and surface melt on the ice interior via brine dynamics is also likely to increase the preconditioning from large-scale ice strength heterogeneity early in the thaw season, which also could affect the landfast ice variability and the timing of its breakup. The mushy-layer thermodynamics thus presents itself as a useful, if not necessary, step towards improving the coupling between the thermodynamic and dynamic sea ice model components.

Appendix A: Frazil formation from the standard mushy-layer congelation parameterization

In Icepack, the mushy-layer congelation parameterization is composed of two components: the downward migration of the ice–ocean boundary, based on the computed congelation rate, and the integration of mass (seawater and/or fresh ice) in the bottom sea ice layer. In mushy-layer physics, the ice–ocean interface is defined by the position where the mush medium reaches a 100 % liquid fraction. Accordingly, the standard congelation parameterization assumes that its downward migration precedes the freezing of seawater, such that seawater without fresh ice is being added in the mush medium. Later, solidification in the bottom ice layer occurs in subsequent time steps when solving for the internal temperature profiles via the liquidus relation. This, however, implies that the enthalpy used to define the congelation rate (i.e., Eq. 16, based on initial congelation liquid fraction ϕinit) differs from the energy actually being integrated in the sea ice by congelation.

Figure A1Total basal ice growth (hc+f, solid lines) and contributions from congelation (hc, dashed lines) and frazil formation (hf, dot-dashed lines) in mushy simulations using the standard congelation scheme, with different initial congelation liquid fractions: default ϕinit=0.85 in blue and ϕinit=0.65 in green.


Specifically, defining EA=Fbot-Fcb as the energy available for congelation (i.e., the conductive heat imbalance at the ice–ocean interface) and EU as the energy used during congelation (i.e., from Eq. 17), the fraction r of the available energy used by the mushy-layer congelation can be written as

(A1) r = E U E A = - q w h c t F bot - F cb ,

where hc is the congelation ice thickness and qw=cwρwTf is the enthalpy of seawater at the freezing point. Using Eq. (16), this reduces to

(A2) r = q w - L 0 ρ i ( 1 - ϕ init ) = - c 1 T 1 - ϕ init ,

where c1=ρwcw/L0ρi∼0.014 and T is in degrees Celsius. This demonstrates that unless ϕinit is close to 1, the congelation only accounts for a fraction r≪1 of the heat flux imbalance (e.g., using ϕinit=0.85 and T-1.8°C, we find r=0.17), and the remaining energy is taken from the ocean. If the ocean is at the freezing point and there is no heat transfer from below the mixed layer, this energy transfer leads to frazil formation. The rate of frazil formation can be estimated using Eq. (19) with Focn=(1-r)UA, as

(A3) h f t = 1 - r - q f E A .

The total basal ice growth in the standard mushy-layer physics is then obtained by adding Eqs. (16) and (A3):

(A4) h t = h c t + h f t = r - q w E A + 1 - r - q f E A ,

where Eq. (A2) has been used to rewrite the growth in terms of the available energy and the fraction r. This indicates that increasing the fraction r (e.g., by decreasing ϕinit; see Eq. A2) decreases the congelation rate but also increases the amount of frazil formation by a proportional amount. Varying ϕinit thus results in a similar total basal growth (Fig. A1).

Appendix B: Modified mushy-layer congelation parameterization

Here, we propose a modified mushy-layer congelation scheme that aims to reduce the amount of frazil formation. The modifications are twofold: (1) the congelation rate is defined by the energy needed to bring seawater enthalpy to that of the integrated mushy layer and (2) the mushy layer integrated in the lowest ice layer has a liquid fraction ϕinit. This implies that some solidification occurs simultaneously as the ice–ocean interface migrates downward.

Specifically, instead of Eqs. (16) and (17), we use Eqs. (20) and (21). The energy integrated in the bottom ice layer in this case corresponds to

(B1) E U = - q m h c t ,

and the fraction r of the available energy used by the modified mushy-layer congelation is (instead of Eq. A2) represented as

(B2) r = - q w h c t F bot - F cb = q m q m - q w .

Given that qmqw (qm and qw being negative), we have r∼1, and the volume of frazil associated with congelation is negligible.

Figure B1Total basal ice growth (hc+f, solid lines) and contributions from congelation (hc, dashed lines) and frazil formation (hf, dot-dashed lines) in mushy simulations using the standard (blue) and modified (green) congelation schemes, both with default ϕinit=0.85. Using the modified congelation scheme, the total basal growth and congelation lines are superposed as the frazil formation is zero.


Figure B2Total basal ice growth (hc+f, solid lines) and contributions from congelation (hc, dashed lines) and frazil formation (hf, dot-dashed lines) in mushy simulations using the modified congelation scheme, with a different initial congelation liquid fraction. The total basal growth and congelation lines are superposed as the frazil formation is zero.


Using ϕinit=0.85, the modified congelation schemes produce total basal growth rates similar to the ones simulated by the standard parameterization, but all of the growth is attributed to the congelation as there is only a negligible frazil formation (Fig. B1). The sensitivity to the parameter ϕinit is, however, increased (Fig. B2), as the changes in ice congelation rates are no longer balanced by changes in frazil formation. This allows for better tuning with the observations.

Note that in this analysis, we define Sbr= SSS to satisfy the liquidus at the boundary where T=Tf. This implies some salt rejection associated with congelation, and it should be treated accordingly when coupling with an ocean model. To keep the congelation salt-conserving, the brine salinity of the integrated mush layer could be set to Sbr=SSS/ϕinit. Note, however, that this does not satisfy the liquidus at boundary and would thus affect the simulated temperature in the lowest ice layer.

Code and data availability

All codes and data (model and analysis) are available on Zenodo as follows. The Icepack version 1 releases are available on (Hunke et al.2023). The modifications and code used to produce the simulations are available on (Plante et al.2024b). The simulation diagnostics, buoy data, surface retrieval algorithm and all data analysis code are available on (Plante et al.2024a).

Author contributions

MP adapted Icepack v1.1.0 to run in stand-alone mode using GDPS atmospheric forcing and produced the simulations with assistance from JFL, FR and FD. AT and JA deployed the ice mass balance buoys. MP coded the surface retrieval algorithm with contributions from GS. MP and AKT coded the modified congelation scheme. MP, JFL, LBT, FR, FD and GS analyzed and discussed the results. MP wrote the paper with edits from JFL, LBT, FR, GS, AT, AKT and FD.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We thank the Nunatsiavut Research Centre for assistance and support for the ice mass balance buoy deployment and retrieval. We also thank Elizabeth Hunke, David Bailey, David Clemens-Sewall and Andrew Roberts for useful discussions on the mushy-layer sea ice congelation scheme.

Review statement

This paper was edited by John Yackel and reviewed by two anonymous referees.


Bailey, D. A., Holland, M. M., DuVivier, A. K., Hunke, E. C., and Turner, A. K.: Impact of a New Sea Ice Thermodynamic Formulation in the CESM2 Sea Ice Component, J. Adv. Model. Earth Sy., 12, e2020MS002154,, 2020. a, b, c

Barber, D., Hanesiak, J., Chan, W., and Piwowar, J.: Sea‐ice and meteorological conditions in Northern Baffin Bay and the North Water polynya between 1979 and 1996, Atmos. Ocean, 39, 343–359,, 2001. a

Bitz, C. M. and Lipscomb, W. H.: An energy-conserving thermodynamic model of sea ice, J. Geophys. Res.-Oceans, 104, 15669–15677,, 1999. a, b, c, d, e, f

Buehner, M., Morneau, J., and Charette, C.: Four-dimensional ensemble-variational data assimilation for global deterministic weather prediction, Nonlin. Processes Geophys., 20, 669–682,, 2013. a

Buehner, M., McTaggart-Cowan, R., Beaulne, A., Charette, C., Garand, L., Heilliette, S., Lapalme, E., Laroche, S., Macpherson, S. R., Morneau, J., and Zadra, A.: Implementation of Deterministic Weather Forecasting Systems Based on Ensemble–Variational Data Assimilation at Environment Canada. Part I: The Global System, Mon. Weather Rev., 143, 2532–2559,, 2015. a, b

Caixin, W., Bin, C., Keguang, W., Sebastian, G., and Olga, P.: Modelling snow ice and superimposed ice on landfast sea ice in Kongsfjorden, Svalbard, Polar Res., 34,, 2015. a

Carmack, E. C. and Macdonald, R.: Oceanography of the Canadian Shelf of the Beaufort Sea: A Setting for Marine Life, Arctic, 55, 29–45, 2002. a

Cheng, Y., Cheng, B., Zheng, F., Vihma, T., Kontu, A., Yang, Q., and Liao, Z.: Air/snow, snow/ice and ice/water interfaces detection from high-resolution vertical temperature profiles measured by ice mass-balance buoys on an Arctic lake, Ann. Glaciol., 61, 309–319,, 2020. a, b, c, d, e, f, g

Collins, W. D., Bitz, C. M., Blackmon, M. L., Bonan, G. B., Bretherton, C. S., Carton, J. A., Chang, P., Doney, S. C., Hack, J. J., Henderson, T. B., Kiehl, J. T., Large, W. G., McKenna, D. S., Santer, B. D., and Smith, R. D.: The Community Climate System Model Version 3 (CCSM3), J. Climate, 19, 2122–2143,, 2006. a

Cooley, S. W., Ryan, J. C., Smith, L. C., Horvat, C., Pearson, B., Dale, B., and Lynch, A. H.: Coldest Canadian Arctic communities face greatest reductions in shorefast sea ice, Nat. Clim. Change, 10, 533–538,, 2020. a

Côté, J., Desmarais, J.-G., Gravel, S., Méthot, A., Patoine, A., Roch, M., and Staniforth, A.: The Operational CMC–MRB Global Environmental Multiscale (GEM) Model. Part II: Results, Mon. Weather Rev., 126, 1397–1418,<1397:TOCMGE>2.0.CO;2, 1998a. a

Côté, J., Gravel, S., Méthot, A., Patoine, A., Roch, M., and Staniforth, A.: The Operational CMC–MRB Global Environmental Multiscale (GEM) Model. Part I: Design Considerations and Formulation, Mon. Weather Rev., 126, 1373–1395,<1373:TOCMGE>2.0.CO;2, 1998b. a

Dammann, D. O., Eriksson, L. E. B., Mahoney, A. R., Eicken, H., and Meyer, F. J.: Mapping pan-Arctic landfast sea ice stability using Sentinel-1 interferometry, The Cryosphere, 13, 557–577,, 2019. a

Divine, D. V., Korsnes, R., and Makshtas, A. P.: Temporal and spatial variation of shore-fast ice in the Kara Sea, Cont. Shelf Res., 24, 1717–1736,, 2004. a

Duarte, P., Sundfjord, A., Meyer, A., Hudson, S. R., Spreen, G., and Smedsrud, L. H.: Warm Atlantic Water Explains Observed Sea Ice Melt Rates North of Svalbard, J. Geophys. Res.-Oceans, 125, e2019JC015662,, 2020. a, b, c, d, e, f

Dumont, D., Gratton, Y., and Arbetter, T. E.: Modeling the Dynamics of the North Water Polynya Ice Bridge, J. Phys. Oceanogr., 39, 1448–1461,, 2009. a

Dumont, D., Gratton, Y., and Arbetter, T. E.: Modeling Wind-Driven Circulation and Landfast Ice-Edge Processes during Polynya Events in Northern Baffin Bay, J. Phys. Oceanogr., 40, 1356–1372,, 2010. a, b

Dupont, F., Dumont, D., Lemieux, J.-F., Dumas-Lefebvre, E., and Caya, A.: A probabilistic seabed–ice keel interaction model, The Cryosphere, 16, 1963–1977,, 2022. a

DuVivier, A. K., Holland, M. M., Landrum, L., Singh, H. A., Bailey, D. A., and Maroon, E. A.: Impacts of Sea Ice Mushy Thermodynamics in the Antarctic on the Coupled Earth System, Geophys. Res. Lett., 48, e2021GL094287,, 2021. a, b

Eicken, H., Fischer, H., and Lemke, P.: Effects of the snow cover on Antarctic sea ice and potential modulation of its response to climate change, Ann. Glaciol., 21, 369–376,, 1995. a

Eicken, H., Jones, J., Meyer, F., Mahoney, A., Druckenmiller, M. L., Rohith, M., and Kambhamettu, C.: Environmental Security in Arctic Ice-Covered Seas: From Strategy to Tactics of Hazard Identification and Emergency Response, Mar. Technol. Soc. J., 45, 37–48,, 2011. a

Feltham, D. L., Untersteiner, N., Wettlaufer, J. S., and Worster, M. G.: Sea ice is a mushy layer, Geophys. Res. Lett., 33, L14501,, 2006. a, b, c

Flocco, D., Feltham, D. L., and Turner, A. K.: Incorporation of a physically based melt pond scheme into the sea ice component of a climate model, J. Geophys. Res.-Oceans, 115, C08012,, 2010. a

Galley, R. J., Else, B. G. T., Howell, S. E. L., Lukovich, J. V., and Barber, D. G.: Landfast Sea Ice Conditions in the Canadian Arctic: 1983–2009, Arctic, 65, 133–144, 2012. a

Gearheard, S., Matumeak, W., Angutikjuaq, I., Maslanik, J., Huntington, H. P., Leavitt, J., Kagak, D. M., Tigullaraq, G., and Barry, R. G.: “It's Not that Simple”: A Collaborative Comparison of Sea Ice Environments, Their Uses, Observed Changes, and Adaptations in Barrow, Alaska, USA, and Clyde River, Nunavut, Canada, AMBIO, 35, 203–211,[203:INTSAC]2.0.CO;2, 2006. a

Gough, A. J., Mahoney, A. R., Langhorne, P. J., Williams, M. J., Robinson, N. J., and Haskell, T. G.: Signatures of supercooling: McMurdo Sound platelet ice, J. Glaciol., 58, 38–50,, 2012. a

Granskog, M. A., Rösel, A., Dodd, P. A., Divine, D., Gerland, S., Martma, T., and Leng, M. J.: Snow contribution to first-year and second-year Arctic sea ice mass balance north of Svalbard, J. Geophys. Res.-Oceans, 122, 2539–2549,, 2017. a

Gultepe, I., Isaac, G. A., Williams, A., Marcotte, D., and Strawbridge, K. B.: Turbulent heat fluxes over leads and polynyas, and their effects on arctic clouds during FIRE.ACE: Aircraft observations for April 1998, Atmos. Ocean, 41, 15–34,, 2003. a

Hibler, W. D.: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815–846, 1979. a

Holland, M. M., Bailey, D. A., Briegleb, B. P., Light, B., and Hunke, E.: Improved Sea Ice Shortwave Radiation Physics in CCSM4: The Impact of Melt Ponds and Aerosols on Arctic Sea Ice, J. Climate, 25, 1413–1430,, 2012. a

Howell, S. E. L., Wohlleben, T., Dabboor, M., Derksen, C., Komarov, A., and Pizzolato, L.: Recent changes in the exchange of sea ice between the Arctic Ocean and the Canadian Arctic Archipelago, J. Geophys. Res.- Oceans, 118, 3595–3607,, 2013. a

Hunke, E. C. and Dukowicz, J. K.: An Elastic–Viscous–Plastic Model for Sea Ice Dynamics, J. Phys. Oceanogr., 27, 1849–1867,<1849:AEVPMF>2.0.CO;2, 1997. a

Hunke, E. C., Lipscomb, W. H., Turner, A. K., Jeffery, N., and Elliott, S.: Cice: the los alamos sea ice model documentation and software user’s manual version 4.1 la-cc-06-012, T-3 Fluid Dynamics Group, Los Alamos National Laboratory, 675, 500, 2010. a

Hunke, E. C., Hebert, D. A., and Lecomte, O.: Level-ice melt ponds in the Los Alamos sea ice model, CICE, Ocean Model., 71, 26–42,, 2013. a

Hunke, E. C., Allard, R., Bailey, D. A., Blain, P., Craig, A., Dupont, F., DuVivier, A. K., Grumbine, R., Hebert, D. A., Holland, M. M., Jeffery, N., Lemieux, J.-F., Osinski, R., Rasmussen, T., Ribergaard, M., Roach, L., Roberts, A., Turner, M., and Winton, M.: CICE-Consortium/Icepack: Icepack 1.4.0 (1.4.0), Zenodo [code],, 2023. a, b

Huwald, H., Tremblay, L.-B., and Blatter, H.: A multilayer sigma-coordinate thermodynamic sea ice model: Validation against Surface Heat Budget of the Arctic Ocean (SHEBA)/Sea Ice Model Intercomparison Project Part 2 (SIMIP2) data, J. Geophys. Res.-Oceans, 110, C05010,, 2005. a

Jackson, K., Wilkinson, J., Maksym, T., Meldrum, D., Beckers, J., Haas, C., and Mackenzie, D.: A Novel and Low-Cost Sea Ice Mass Balance Buoy, J. Atmos. Ocean. Tech., 30, 2676–2688,, 2013. a, b, c

King, J., Howell, S., Brady, M., Toose, P., Derksen, C., Haas, C., and Beckers, J.: Local-scale variability of snow density on Arctic sea ice, The Cryosphere, 14, 4323–4339,, 2020. a

Kirillov, S., Babb, D. G., Komarov, A. S., Dmitrenko, I., Ehn, J. K., Worden, E., Candlish, L., Rysgaard, S., and Barber, D. G.: On the Physical Settings of Ice Bridge Formation in Nares Strait, J. Geophys. Res.-Oceans, 126, e2021JC017331,, 2021. a

Kwok, R.: Exchange of sea ice between the Arctic Ocean and the Canadian Arctic Archipelago, Geophys.Res. Lett., 33, L16501,, 2006. a

Lemieux, J.-F., Tremblay, L. B., Dupont, F., Plante, M., Smith, G. C., and Dumont, D.: A basal stress parameterization for modeling landfast ice, J. Geophys. Res.-Oceans, 120, 3157–3173,, 2015. a

Lemieux, J.-F., Dupont, F., Blain, P., Roy, F., Smith, G. C., and Flato, G. M.: Improving the simulation of landfast ice by combining tensile strength and a parameterization for grounded ridges, J. Geophys. Res.-Oceans, 121, 7354–7368,, 2016. a, b, c

Liao, Z., Cheng, B., Zhao, J., Vihma, T., Jackson, K., Yang, Q., Yang, Y., Zhang, L., Li, Z., Qiu, Y., and Cheng, X.: Snow depth and ice thickness derived from SIMBA ice mass balance buoy data using an automated algorithm, Int. J. Digit. Earth, 12, 962–979,, 2019. a, b, c, d, e, f

Liu, Y., Losch, M., Hutter, N., and Mu, L.: A New Parameterization of Coastal Drag to Simulate Landfast Ice in Deep Marginal Seas in the Arctic, J. Geophys. Res.-Oceans, 127, e2022JC018413,, 2022. a, b

Lüpkes, C., Vihma, T., Birnbaum, G., and Wacker, U.: Influence of leads in sea ice on the temperature of the atmospheric boundary layer during polar night, Geophys. Res. Lett., 35, L03805,, 2008. a

Madec, G. and the NEMO team: NEMO ocean engine, Note du Pole de modélisation 27, Institut Pierre-Simon Laplace (IPSL), France, ISSN 1288-1619, 2008. a

Madec, G., Delecluse, P., Imbard, M., and Levy, C.: OPA 8 Ocean General Circulation Model – Reference Manual, Tech. rep., LODYC/IPSL Note 11, 1998. a

Maksym, T. and Jeffries, M. O.: A one-dimensional percolation model of flooding and snow ice formation on Antarctic sea ice, J. Geophys. Res.-Oceans, 105, 26313–26331,, 2000. a

Maykut, G. A. and Untersteiner, N.: Some results from a time-dependent thermodynamic model of sea ice, J. Geophys. Res., 76, 1550–1575,, 1971. a

Melling, H.: Sea ice of the northern Canadian Arctic Archipelago, J. Geophys. Res.-Oceans, 107, 2-1–2-21,, 2002. a

Melling, H., Gratton, Y., and Ingram, G.: Ocean circulation within the North Water polynya of Baffin Bay, Atmos. Ocean, 39, 301–325,, 2001. a

Merkouriadi, I., Liston, G. E., Graham, R. M., and Granskog, M. A.: Quantifying the Potential for Snow-Ice Formation in the Arctic Ocean, Geophys. Res. Lett., 47, e2019GL085020,, 2020. a

Notz, D. and Worster, M. G.: Desalination processes of sea ice revisited, J. Geophys. Res.-Oceans, 114, C05006,, 2009. a, b, c

Planck, C. J., Whitlock, J., Polashenski, C., and Perovich, D.: The evolution of the seasonal ice mass balance buoy, Cold Reg. Sci. Technol., 165, 102792,, 2019. a

Plante, M., Tremblay, B., Losch, M., and Lemieux, J.-F.: Landfast sea ice material properties derived from ice bridge simulations using the Maxwell elasto-brittle rheology, The Cryosphere, 14, 2137–2157,, 2020. a

Plante, M., Lemieux, J.-F., Tremblay, L. B., Tivy, A., Angnatok, J., Roy, F., Smith, G., Dupont, F., and Turner, A. K.: Using Icepack to reproduce Ice Mass Balance buoy observations in landfast ice – data and codes, Zenodo [code, data set],, 2024a. a

Plante, M., Allard, R., Angnatok, J., Bailey, D. A., Blain, P., Craig, A., Dupont, F., DuVivier, A., Grumbine, R., Hebert, D., Holland, M., Hunke, E., Jeffery, N., Lemieux, J.-F., Osinski, R., Rassmussen, T., Ribergaard, M., Roach, L., Roberts, A., Roy, F., Turner, A. K., Turner, M., Winton, M., Tremblay, L. B., Tivy, A., and Smith, G.: mathieuslplante/Icepack: Modified Icepack v.1.1.0 code used in Plante et al. (2024), Zenodo [code],, 2024b. a

Pringle, D. J., Eicken, H., Trodahl, H. J., and Backstrom, L. G. E.: Thermal conductivity of landfast Antarctic and Arctic sea ice, J. Geophys. Res.-Oceans, 112, C04017,, 2007. a

Provost, C., Sennéchael, N., Miguet, J., Itkin, P., Rösel, A., Koenig, Z., Villacieros-Robineau, N., and Granskog, M. A.: Observations of flooding and snow-ice formation in a thinner Arctic sea-ice regime during the N-ICE2015 campaign: Influence of basal ice melt and storms, J. Geophys. Res.-Oceans, 122, 7115–7134,, 2017. a, b, c, d, e, f, g, h

Raddatz, R. L., Asplin, M. G., Candlish, L., and Barber, D. G.: General Characteristics of the Atmospheric Boundary Layer Over a Flaw Lead Polynya Region in Winter and Spring, Bound.-Lay. Meteorol., 138, 321–335,, 2011. a

Rampal, P., Bouillon, S., Ólason, E., and Morlighem, M.: neXtSIM: a new Lagrangian sea ice model, The Cryosphere, 10, 1055–1073,, 2016. a

Richter, M. E., Leonard, G. H., Smith, I. J., Langhorne, P. J., Mahoney, A. R., and Parry, M.: Accuracy and precision when deriving sea-ice thickness from thermistor strings: a comparison of methods, J. Glaciol., 69, 879–898,, 2023. a, b

Richter-Menge, J. A., Perovich, D. K., Elder, B. C., Claffey, K., Rigor, I., and Ortmeyer, M.: Ice mass-balance buoys: a tool for measuring and attributing changes in the thickness of the Arctic sea-ice cover, Ann. Glaciol., 44, 205–210,, 2006. a

Rösel, A., Itkin, P., King, J., Divine, D., Wang, C., Granskog, M. A., Krumpen, T., and Gerland, S.: Thin Sea Ice, Thick Snow, and Widespread Negative Freeboard Observed During N-ICE2015 North of Svalbard, J. Geophys. Res.-Oceans, 123, 1156–1176,, 2018. a, b

Semtner, A. J.: A Model for the Thermodynamic Growth of Sea Ice in Numerical Investigations of Climate, J. Phys. Oceanogr., 6, 379–389,<0379:AMFTTG>2.0.CO;2, 1976. a

Smith, G. C., Roy, F., Mann, P., Dupont, F., Brasnett, B., Lemieux, J.-F., Laroche, S., and Bélair, S.: A new atmospheric dataset for forcing ice–ocean models: Evaluation of reforecasts using the Canadian global deterministic prediction system, Q. J. Roy. Meteor. Soc., 140, 881–894,, 2014. a, b

Smith, G. C., Bélanger, J.-M., Roy, F., Pellerin, P., Ritchie, H., Onu, K., Roch, M., Zadra, A., Colan, D. S., Winter, B., Fontecilla, J.-S., and Deacu, D.: Impact of Coupling with an Ice–Ocean Model on Global Medium-Range NWP Forecast Skill, Mon. Weather Rev., 146, 1157–1180,, 2018. a

Smith, G. C., Liu, Y., Benkiran, M., Chikhar, K., Surcel Colan, D., Gauthier, A.-A., Testut, C.-E., Dupont, F., Lei, J., Roy, F., Lemieux, J.-F., and Davidson, F.: The Regional Ice Ocean Prediction System v2: a pan-Canadian ocean analysis system using an online tidal harmonic analysis, Geosci. Model Dev., 14, 1445–1467,, 2021. a

Stirling, I.: The Biological Importance of Polynyas in the Canadian Arctic, Arctic, 33, 303–315, 1980. a

Stirling, I.: The importance of polynyas, ice edges, and leads to marine mammals and birds, J. Mar. Syst., 10, 9–21,, 1997. a

Tian, Z., Cheng, B., Zhao, J., Vihma, T., Zhang, W., Li, Z., and Zhang, Z.: Observed and modelled snow and ice thickness in the Arctic Ocean with CHINARE buoy data, Acta Oceanol. Sin., 36, 66–75,, 2017. a, b

Tremblay, J.-E., Gratton, Y., Carmack, E. C., Payne, C. D., and Price, N. M.: Impact of the large-scale Arctic circulation and the North Water Polynya on nutrient inventories in Baffin Bay, J. Geophys. Res.-Oceans, 107, 26-1–26-14,, 2002. a

Tremblay, L.-B. and Mysak, L. A.: Modeling Sea Ice as a Granular Material, Including the Dilatancy Effect, J. Phys. Oceanogr., 27, 2342–2360,<2342:MSIAAG>2.0.CO;2, 1997. a

Trodahl, H. J., McGuinness, M. J., Langhorne, P. J., Collins, K., Pantoja, A. E., Smith, I. J., and Haskell, T. G.: Heat transport in McMurdo Sound first-year fast ice, J. Geophys. Res.-Oceans, 105, 11347–11358,, 2000.  a

Turner, A. K. and Hunke, E. C.: Impacts of a mushy-layer thermodynamic approach in global sea-ice simulations using the CICE sea-ice model, J. Geophys. Res.-Oceans, 120, 1253–1275,, 2015. a, b, c

Turner, A. K., Hunke, E. C., and Bitz, C. M.: Two modes of sea-ice gravity drainage: A parameterization for large-scale modeling, J. Geophys. Res.-Oceans, 118, 2279–2294,, 2013. a, b, c, d, e, f, g, h

Ungermann, M., Tremblay, L. B., Martin, T., and Losch, M.: Impact of the ice strength formulation on the performance of a sea ice thickness distribution model in the Arctic, J. Geophys. Res.-Oceans, 122, 2090–2107,, 2017. a

West, A., Collins, M., and Blockley, E.: Using Arctic ice mass balance buoys for evaluation of modelled ice energy fluxes, Geosci. Model Dev., 13, 4845–4868,, 2020. a

Wilchinsky, A. V. and Feltham, D. L.: A continuum anisotropic model of sea-ice dynamics, P. Roy. Soc. Lond. A Mat., 460, 2105–2140,, 2004. a

Wongpan, P., Hughes, K. G., Langhorne, P. J., and Smith, I. J.: Brine Convection, Temperature Fluctuations, and Permeability in Winter Antarctic Land-Fast Sea Ice, J. Geophys. Res.-Oceans, 123, 216–230,, 2018. a

Zuo, G., Dou, Y., and Lei, R.: Discrimination Algorithm and Procedure of Snow Depth and Sea Ice Thickness Determination Using Measurements of the Vertical Ice Temperature Profile by the Ice-Tethered Buoys, Sensors, 18, 4162,, 2018. a

Short summary
We use a sea ice model to reproduce ice growth observations from two buoys deployed on coastal sea ice and analyze the improvements brought by new physics that represent the presence of saline liquid water in the ice interior. We find that the new physics with default parameters degrade the model performance, with overly rapid ice growth and overly early snow flooding on top of the ice. The performance is largely improved by simple modifications to the ice growth and snow-flooding algorithms.