Investigating snow sinks on level sea ice: A case study in the western Arctic

Abstract

SnowModel-LG reconstructs snow depth and density over sea ice, explicitly resolving important snow sinks like blowing snow sublimation, static surface sublimation and melt, but not snow-ice formation. To examine snow sinks on level sea ice, we coupled SnowModel-LG with HIGHTSI, a 1-D thermodynamic sea-ice model, to create SMLG_HS. SMLG_HS simulations of snow depth and level ice thickness were evaluated against high-resolution airborne observations from the western Arctic, highlighting the importance of snow mass redistribution processes, i.e. snow’s tendency to leave level ice and accumulate over deformed ice due to wind-induced redistribution. Not accounting for snow mass redistribution, SMLG_HS overestimates snow depth on level ice, resulting in underestimation of level ice thickness and overestimation of snow-ice thickness. Our case study shows that snow depth on level ice needs to be reduced by 40% to simulate both snow depth and level ice thickness realistically in the western Arctic in April 2017. An independent analysis of snow volume distribution between level and deformed sea ice using airborne radar observations supported the model results and revealed a linear relationship that enables estimating the amount of snow remaining on level ice at the end of winter based on the amount of ice deformation.

1. Introduction

Arctic sea ice is going through unprecedented changes, decreasing both in extent (e.g. Stroeve and others, 2014) and in thickness (Maslanik and others, 2007; Kwok and others, 2009), and transitioning from a multi-year ice to a seasonal, first-year ice system (Meier and others, 2014). The role of snow in sea-ice mass balance is becoming predominant in many ways, because of the higher sensitivity of sea ice to its environmental conditions. The thermal resistance of snow cover significantly reduces the atmosphere-ocean heat fluxes, regulating sea-ice growth in winter (Maykut, 1978; Ledley, 1991). The high snow albedo reflects most of the solar radiation back to space, delaying sea ice from melting in spring (Perovich and others, 2017). The snow load may submerge thinner ice underneath the water level, creating negative freeboard conditions (Granskog and others, 2017; Merkouriadi and others, 2020). If sea water floods at the ice/snow interface and freezes there, snow ice is formed that is a mixture of frozen seawater and snow (e.g. Leppäranta, 1983) and increases the thickness of the sea ice. Snow ice is a common phenomenon in seas that are seasonally covered by ice (i.e. Baltic Sea, Sea of Okhotsk) and in large parts of the Antarctic sea ice (Massom and others, 2001), but it was not commonly observed in situ in drifting Arctic sea ice until the Norwegian Young Sea ICE (N-ICE2015) expedition (Granskog and others, 2017; Provost and others, 2017). Snow ice is a sink for snow, and it can positively contribute to the sea-ice mass balance (Merkouriadi and others, 2017; 2020), which has implications for remote sensing retrievals of sea-ice thickness. Therefore, it is essential to consider it for understanding sea-ice mass balance, both in contemporary times in peripheral seas and in future scenarios where sea ice may be thinner than present-day conditions.

Accounting for different snow processes is also relevant for remote sensing applications. Satellite altimetry is the most common method for monitoring sea-ice thickness, providing nearly full coverage of the Arctic Ocean (Laxon and others, 2003; Markus and others, 2017; Landy and others, 2022). Information on the snow load on sea ice is crucial for accurate altimetry retrievals of sea-ice thickness, because radar and laser altimeters measure the elevation of the ice or snow surface from the water surface, i.e. ice or snow freeboard. Snow depth and density are required to convert freeboard to sea-ice thickness information (e.g. Laxon and others, 2003). According to Giles and others (2007), uncertainties in snow depth and density contribute 48% and 14%, respectively, to the total error of sea-ice thickness retrievals from radar altimetry. A more recent study by Landy and others (2020) estimated these uncertainties at 11% for snow depth and 16% for density. Similarly, snow depth and density uncertainties were found to contribute 70% and 30–35%, respectively, to the total error of sea-ice thickness retrievals from laser altimetry (Zygmuntowska and others, 2014).

Snow depth and density estimates used in satellite altimetry applications are often derived from snow climatologies or modified versions of snow climatologies from historical observations. The most widely used snow-on-sea-ice climatology is compiled from a snow depth and density dataset collected mostly over multi-year ice in 1954–91 (Warren and others, 1999). In a changing Arctic sea-ice system, snow conditions are expected to change as well (Blanchard-Wrigglesworth and others, 2015; Webster and others, 2021), and these changes are not captured by the Warren and others (1999) climatology. In addition to the long-term changes, climatological data are not suitable for representing the spatiotemporal variability of the snow conditions in the Arctic, which are evidently strong (Warren and others, 1999; Webster and others, 2024). To account for spatiotemporal variability, efforts have focused on reanalysis-based snow depth and density reconstructions (e.g. Kwok and Cunningham, 2008; Blanchard-Wrigglesworth and others, 2018; Petty and others, 2018), i.e. pan-Arctic simulations of snow depth and density evolution on sea ice. A recent contribution was SnowModel-LG, a state-of-the-art Lagrangian snow evolution model (Liston and others, 2020a). Compared to other reanalysis-based products, SnowModel-LG implemented Lagrangian parcel tracking and included an improved representation of snow evolution physics. It has been bias-corrected and validated against a wide observation framework in all seasons and yielded good agreement, especially with in situ measurements (Stroeve and others, 2020).

SnowModel-LG explicitly resolves many snow mass sources and sinks, such as blowing snow sublimation, static-surface sublimation and melt, by performing a snow mass-budget calculation in each time step (Liston and others, 2020a). However, SnowModel-LG, similarly to all the abovementioned, operationally used, pan-Arctic snow models, is not coupled to a sea-ice model. Therefore, it does not account for snow sinks caused by snow-ice formation. Moreover, being configured over ice parcels of kilometer-scale (25 km × 25 km), it does not resolve wind-driven snow mass redistribution. The latter describes the tendency of snow to accumulate on the lee side of pressure ridges and other roughness elements (e.g. Liston and others, 2018) as a result of snow redistribution by the wind. This process results in uneven snow load over a sea-ice floe (i.e. reduced snow over level ice areas and increased snow over deformed ice) (Sturm and others, 2002; Webster and others, 2015; Itkin and others, 2023). Because in this study we are examining level ice only, we will be referring to the sub-parcel snow mass redistribution process as a snow sink. In other words, the term ‘snow sinks’ here refers to all processes that reduce snow mass on level ice, including snow sublimation, melting, snow-ice formation and snow mass redistribution.

To investigate snow sinks on level Arctic sea ice, we coupled SnowModel-LG with the High-Resolution Thermodynamic Sea Ice model (HIGHTSI) (Launiainen and Cheng, 1998) to produce SMLG_HS. In SMLG_HS, snow-ice forms when the ice surface is depressed below the water surface (negative freeboard). SMLG_HS outputs of snow depth, snow-ice and sea-ice thickness from 1 August 2007 until 31 July 2021 were evaluated against airborne observations in the western Arctic to examine and to mitigate the biases introduced when sub-parcel snow mass redistribution processes are ignored.

2. Materials and methods

2.1. SnowModel-LG

SnowModel is a collection of snow distribution and snow evolution modeling tools, applicable to any environment experiencing snow, including sea-ice applications (Liston and Elder, 2006a; Liston and others, 2018). SnowModel-LG is adapted for snow depth and density reconstruction over sea ice (Liston and others, 2020a). It is implemented in a Lagrangian framework to simulate snow properties on drifting sea-ice parcels. SnowModel-LG accounts for physical snow processes such as sublimation from static surfaces and blowing snow, snow melt, evolution of snow density and temperature profiles, energy and mass transfers within the snowpack, and superimposed ice formation in a multi-layer configuration. The ice parcels are 1-D, and they do not interact with each other.

At each time step (3 hours here), SnowModel-LG performs a mass-budget calculation, where snow water equivalent (SWE) depth (m) is defined by snow mass gains, losses and ice parcel dynamics,

(1)\begin{equation} \frac{\mathrm{d}\mathrm{SWE}}{\mathrm{d}t} = \frac{1}{\rho_\mathrm{w}} \left[ \left( P_\mathrm{r} + P_\mathrm{s} \right) - \left( S_\mathrm{ss} + S_\mathrm{bs} + M \right) + D \right] \end{equation}

where t (s) is time; $\rho_\mathrm{w}=1000\,\mathrm{kg\,m}^{-3}$ is the water density; $P_\mathrm{r}$ (kg m⁻² s⁻¹) and $P_\mathrm{s}$ (kg m⁻² s⁻¹) are the water-equivalent rainfall and snowfall fluxes, respectively; $S_\mathrm{ss}$ (kg m⁻² s⁻¹) and $S_\mathrm{bs}$ (kg m⁻² s⁻¹) are the water-equivalent sublimation from static-surface and blowing-snow processes, respectively; M (kg m⁻² s⁻¹) is the melt-related mass loss; and D (kg m⁻² s⁻¹) represents the mass losses and gains from sea-ice dynamics processes (i.e. parcels being created and lost with ice motion, divergence and convergence).

Snow depth $h_\mathrm{s}$ (m) is related to SWE through the ratio of snow $(\rho_\mathrm{s})$ and water $(\rho_\mathrm{w})$ densities,

(2)\begin{equation} \mathrm{SWE} = \frac{\rho_\mathrm{s}}{\rho_\mathrm{w}} h_\mathrm{s}. \end{equation}

Therefore, the evolution of snow depths and densities are calculated by

(3)\begin{equation} \frac{\mathrm{d}\left(\rho_\mathrm{s} h_\mathrm{s} \right)}{\mathrm{d}t} = \left( P_\mathrm{r} + P_\mathrm{s} \right) - \left( S_\mathrm{ss} + S_\mathrm{bs} + M \right) + D. \end{equation}

In SnowModel-LG, snow density evolves and changes in response to compaction (weight of the above snow layers), wind force, freezing of liquid water and vapor flux through the snowpack. Additional information on the components and the configuration of SnowModel-LG are summarized and provided in great detail in Liston and others (2020a). The model configuration in this study is identical to the one used in Liston and others (2020a), only here we have coupled it to a sea-ice model (Sections 2.2–2.3). According to Stroeve and others (2020), SnowModel-LG performed well in capturing the spatial and seasonal variation of snow distributions, when evaluated against several Arctic datasets, including NASA Operation IceBridge (OIB), ice mass-balance buoys, snow buoys, MagnaProbes and ruler measurements.

In the simulations presented herein, Lagrangian parcel tracking began on 1 August 2007. At the start of the first simulation year, the model assumes no snow atop the sea ice, which is well supported by in situ observations from the contemporary period (Radionov and others, 1997; Chapman-Dutton and Webster, 2024; Webster and others, 2024); the following years carry available snow from 31 July to 1 August. Essential inputs are atmospheric reanalysis estimates of near-surface air temperature, relative humidity, precipitation, wind speed and direction, and sea-ice motion and concentration products, described in detail in Section 2.4.

2.2. HIGHTSI

HIGHTSI is a 1-D thermodynamic sea-ice model designed to simulate the evolution of snow and sea-ice thickness and temperature profiles (Launiainen and Cheng, 1998) by solving the heat conduction equation for multiple ice and snow layers. The sea-ice thermal conductivity is parameterized following Pringle and others (2007). HIGHTSI simulates snow-ice formation following Saloranta (2000).

HIGHTSI has been widely used in process studies and validated extensively against observations (Cheng and others, 2008b; 2013; Wang and others, 2015; Merkouriadi and others, 2017; 2020). In this study, we used a model configuration that is derived from validation studies on Arctic sea ice. The model’s vertical resolution has been found to be critical for its performance in the Arctic (Cheng and others, 2008a). Here, we used 20 layers in the ice which is considered optimal for capturing internal thermodynamic processes (Cheng and others, 2008a; 2008b; 2013; Wang and others, 2015). Detailed information on model parameterizations is given in Table S1 in the supplementary material.

Merkouriadi and others (2020) implemented HIGHTSI in a Lagrangian framework to examine pan-Arctic snow-ice distributions. In the study presented herein, HIGHTSI was modified further, so that snow depth and bulk density evolution were simulated by SnowModel-LG in a 25-layer configuration. We did this because SnowModel-LG provides a more advanced representation of snow physics compared to HIGHTSI’s snow configuration. Additionally, we wanted to explore the effects of snow sinks using a publicly available snow product such as SnowModel-LG.

2.3. SMLG_HS

We performed two separate snow-on-sea-ice simulations. First, we simulated snow depth and density with SnowModel-LG (i.e. Liston and others, 2020a). Second, SnowModel-LG’s snow depth and density evolution were coupled with HIGHTSI’s snow-ice and thermodynamic ice growth representations. Hereafter, we will refer to the original SnowModel-LG output as SMLG and to the coupled output as SMLS_HS.

For the SMLG_HS runs, snow density was simulated following Appendix C of Liston and others (2020a) and stored as a bulk density value. To represent the typical snow stratigraphy of snow on Arctic sea ice (i.e. high-density wind slab layer at the top, low-density depth hoar layer at the bottom), the vertical density profile was parameterized as being a linear fit between densities that are 20% greater than the bulk snow density of SMLG at the top of the snowpack and 20% less at the bottom of the snowpack. These percentages are consistent with snow-pit measurements made during the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) expedition in 2019–20 (Macfarlane and others, 2023). This approach was chosen to provide a best-possible fit to available snow density observations and to account for changes in snow density in response to snow-ice formation. When snow ice was formed, the corresponding snow-depth amount was removed from the lower density bottom layers of the snowpack, and the new bottom density was calculated based on the the linear fit of depth and density of the remaining snow. Additional model specifications are presented in the supplementary material (Table S1).

2.4. Input datasets

Daily ice concentrations (15–100%) by the NASA team algorithm (DiGirolamo and others, 2022) were used to define whether an ice parcel existed and whether snow could accumulate on that parcel. Ice motion vectors from Tschudi and others (2019); (2020) gridded over 25 km spatial resolution were used as Lagrangian ice parcel tracks. NASA’s Modern Era Retrospective Analysis for Research and Application Version 2 (MERRA-2; Global Modeling And Assimilation Office (GMAO), 2015a; 2015b; Gelaro and others, 2017) was used as atmospheric forcing to SMLG_HS. Specifically, SMLG_HS was forced with 10 m wind speed and direction, 2 m air temperature and relative humidity, and total water-equivalent precipitation from MERRA-2. During these simulations, MicroMet (Liston and Elder, 2006b) provided the required liquid and solid precipitation, and the downwelling shortwave and longwave radiation following Liston and others (2020a).

We applied the same bias-correction in MERRA-2 reanalysis as in Liston and others (2020a), where snow depth observations from NASA OIB (2009–16) were used to scale the precipitation inputs. In Liston and others (2020a), 8-year averages of precipitation scaling factors were calculated and they were applied over all ice parcels and through the whole simulation period, making the results of MERRA-2 and the European Centre for Medium-Range Weather Forecasts (ECMWF) ReAnalysis-5th Generation (ERA5; Hersbach and others, 2020) model runs similar. This is why we only use MERRA-2 reanalysis in this study. Scaling factor was 1.37 for MERRA-2, indicating the need to increase the precipitation inputs in order to match the OIB observations (Liston and others, 2020a). The same scaling factor was used in this study for the results to be comparable with the publicly available SMLG snow depth and density dataset (Liston and others, 2020b).

For the ocean boundary forcing, at the ice/ocean interface, we used ocean heat flux from the Ocean Reanalysis System 5 (ORAS5) provided at the ECMWF (Zuo and others, 2019). ORAS5 resolution is eddy-permitting (0.25^∘ latitude and longitude) horizontally and 1 m vertically. ORAS5 includes five ensemble members and covers the period from 1979 onward. In our study, we used the ensemble mean, providing one unique value on a 1^∘ grid for each simulation day.

2.5. Model configuration and outputs

The simulations began on 1 August 2007 and ran through 31 July 2021. Temporal resolution was 3 hours to capture diurnal variations, and the parcel-specific outputs (e.g. snow depth, snow bulk density, sea-ice thickness and snow-ice thickness) were saved at the end of each day. Ice parcel trajectories were linearly interpolated from weekly to daily time steps. On 1 August of each year (except in the first year), the multi-year ice thicknesses were calculated from the sea-ice thickness distribution on 31 July. The initial ice thickness conditions on 1 August 2007 were defined by performing a 1-year simulation with a domain-wide initial condition of 1 m, and then using the ice thickness distribution at the end of the first simulation year as the initial condition for the beginning of the 14 year-simulation (i.e. the model ran the first year twice and assumed the 31 July 2008 ice thickness distribution equaled the 1 August 2007 distribution). In addition, any snow remaining at 00:00 UTC on 1 August (the last time step on 31 July) was used as the initial condition for the following simulation year that started at 03:00 UTC on 1 August (these are the standard model spin-up procedures as implemented in Liston and others (2020a)).

The daily simulation outputs for each parcel (approximately 61 000 parcels each year) were gridded to the 25 km × 25 km Equal-Area Scalable Earth (EASE) grid, provided by the National Snow and Ice Data Center (NSIDC). The location of each parcel was used to calculate the overlap between that parcel and the EASE grid cell, i.e. the fractional area of the EASE grid cell that was occupied by the parcel. The fractional area was then multiplied by the sea-ice concentration of the parcel, and the result was used to weigh the parcels’ contribution to each EASE grid cell. This procedure of area- and concentration-weighted averages within the EASE grid cells conserved the examined parameters, similar to Merkouriadi and others (2020); Liston and others (2020a).

2.6. Evaluation exercise

To evaluate SMLG_HS snow depth and sea-ice thickness, we compared them against a total of >100 airborne surveys from the Alfred Wegener Institute’s (AWI) IceBird and NASA OIB campaigns over the western Arctic in late winter 2009–19 (Fig. 1). Summarizing descriptions of the respective datasets are given in the Sections 2.6.1 and 2.6.2. We averaged the airborne measurements over the same model EASE grid when >50 values were present in a grid cell.

Figure 1. Spatial and annual coverage of the 11 AWI IceBird survey flights in 2017 and 2019 (Table 1) and the 99 NASA Operation IceBridge (OIB) survey flights in 2009–19 (Table A1 in Appendix A). The background shows the average March–April monthly sea-ice concentration in 2009–19.

2.6.1. AWI IceBird

The AWI IceBird program carried out 11 survey flights over the western Arctic Ocean in April 2017 and 2019, monitoring the regional sea-ice conditions in very high resolution (Table 1). The nominal measurement spacing along-track is 5–6 m. Snow depth data were derived from an airborne snow radar similar to OIB using the Peakiness retrieval algorithm (Jutila and others, 2021a; 2021b; 2022b). Sea-ice thickness was derived by subtracting snow depth from the total (sea ice + snow) thickness data measured simultaneously with a towed electromagnetic induction sounding instrument (Jutila and others, 2022a; 2024a; 2024b). We distinguished measurements over level ice by using the flag in the data product that implements a sea-ice thickness gradient threshold of 4 cm within an along-track distance of 1 m over continuous sections of at least 100 m long.

Table 1. Statistics of the 11 AWI IceBird survey flights over the western Arctic Ocean in 2017 and 2019 (Fig. 1) used in this study, where L is the total length of the survey flight, ${\bar{h}}_{\mathrm{s,level}}$ is the average snow depth on level ice, $\bar{h}_{\mathrm{s,deformed}}$ is the average snow depth on deformed ice, $\bar{h}_{\mathrm{s,all}}$ is the average snow depth of the entire survey flight including all ice types, $\frac{\bar{h}_{\mathrm{s,level}}}{\bar{h}_{\mathrm{s,deformed}}}$ is the fraction of the average snow depth on level ice to the average snow depth on deformed ice, $f_{\mathrm{level}}$ is the level ice fraction of the survey flight, $f_{V_{\mathrm{s,level}}}$ is the fraction of snow volume on level ice, $f_{\mathrm{MYI}}$ is the fraction of multi-year ice (MYI) and $f_{\mathrm{NaN}}$ is the fraction of missing snow depth data

^aNot applicable; no sea-ice thickness measurements.

^bUsed in the sensitivity experiment (Section 2.7).

2.6.2. NASA OIB

Annual NASA OIB campaigns over the western Arctic Ocean took place in March–April 2009–19 (MacGregor and others, 2021) and comprise 99 survey flights in total (Table A1 in Appendix A). We used the data products of Kurtz and others (2015); (2016) where the snow depth data were derived from airborne snow radars using the retrieval algorithms described in Kurtz and Farrell (2011); Kurtz and others (2013). The data are averaged in the along-track direction to a 40 m length scale. We did not use the OIB sea-ice thickness to evaluate modeled sea-ice thickness, because it is not directly measured but converted from freeboard and snow depth measurements assuming hydrostatic equilibrium. However, we did use it together with surface roughness data included in the product to guide a level ice identification similar to the IceBird data (Jutila and others, 2022a). To compensate for the increased uncertainty of sea-ice thickness and the approximately 8–10 times coarser along-track resolution compared to the IceBird data, we applied here more strict conditions including a sea-ice thickness gradient threshold of 2 cm m⁻¹ over continuous sections of at least 200 m long as well as ensuring a corresponding surface roughness value of <0.3 m. We determined the numerical values of these conditions through manual iteration and visually inspecting along-track sea-ice transect profiles. Using snow freeboard, snow depth and sea-ice thickness data, we produced figures of snow and sea-ice layer heights with respect to the local sea level (i.e. zero freeboard) in 10 km long along-track segments for five OIB flights from different years and different regions to verify that only truly level ice was identified.

2.7. Sensitivity experiment

Because sub-parcel snow mass redistribution is not considered in the model, we hypothesized that SMLG_HS would overestimate snow depth on level ice and consequently underestimate level ice thickness and overestimate snow-ice thickness. To test our hypothesis, we performed a modeling sensitivity experiment, where we decreased snow depth in SMLG_HS by 10% intervals and compared the simulated snow depth and sea-ice thickness to observations. We conducted in total seven SMLG_HS runs, decreasing the snow depth by 10%, 20%, 30%, 40%, 50%, 60% and 70%, respectively. As observations, we used a subset of IceBird flights in April 2017 that have the highest fractions of level ice and that extend far enough from the vicinity of coastlines, i.e. the edges of the simulation domain (Table 1). The data collected in April 2019 were not included in the experiment due to the smaller fractions of level ice, the presence of deformed multi-year ice and/or the closeness of coastlines. We derived a snow depth fraction that resulted in best fitting of both snow depth and level ice thickness simulations (including snow-ice formation) to the observations. We argue that this snow depth decrease represents the sub-parcel snow mass redistribution sink, because all other snow sinks are accounted for by the model.

As an independent evaluation, we investigated the snow mass redistribution between level and deformed ice also along the 100+ airborne surveys in 2009–19 by calculating the fraction of snow volume on level ice for each flight:

(4a)\begin{equation} V_{\mathrm{s,tot}} = V_{\mathrm{s,level}} + V_{\mathrm{s,deformed}} = \bar{h}_{\mathrm{s,level}} \times f_{\mathrm{level}} + \bar{h}_{\mathrm{s,deformed}} \times \left(1 - f_{\mathrm{level}}\right), \end{equation}

(4b)\begin{equation} f_{V_{\mathrm{s,level}}} = \frac{V_{\mathrm{s,level}}}{V_{\mathrm{s,tot}}} = \frac{\bar{h}_{\mathrm{s,level}} \times f_{\mathrm{level}}}{V_{\mathrm{s,tot}}}, \end{equation}

where $V_{\mathrm{s,tot}}$ is the total snow volume, $V_{\mathrm{s,level}}$ is the snow volume on level ice, $V_{\mathrm{s,deformed}}$ is the snow volume on deformed ice, $\bar{h}_{\mathrm{s,level}}$ is the average snow depth on level ice, $\bar{h}_{\mathrm{s,deformed}}$ is the average snow depth on deformed ice, $f_{\mathrm{level}}$ is the level ice fraction and $f_{V_{\mathrm{s,level}}}$ is the fraction of snow volume on level ice.

3. Results

In the evaluation exercise, we compared SMLG_HS simulations of snow depth and sea-ice thickness against independent airborne observations from the IceBird and OIB campaigns, and we were able to examine snow depth on level ice separately. The results of the evaluation exercise confirmed our hypothesis. They indicated that SMLG_HS overestimated snow depth over level ice on average by 0.06–0.07 m with a root-mean-square error (RMSE) of 0.10–0.11 m, but with an absolute error up to 0.45 m (Fig. 2a–d and Fig. 3a–b). For comparison, in the SMLG, the maximum absolute error was even higher, 0.60 m. Therefore, SMLG_HS underestimated level ice thickness on average by 0.45 m with an RMSE of 0.62 m, but with an absolute error up to 1.76 m (Fig. 2e–h and Fig. 3e). This result was consistent in all IceBird flights examined in the evaluation exercise.

Figure 2. Panels (a)–(d) show the evaluation of modeled snow depth from SMLG and SMLG_HS against airborne radar-derived snow depth measurements from the AWI IceBird survey flight on 8 April 2017. Red color refers to the original SMLG and black color to the new, coupled SMLG_HS. Panels (e)–(h) show the evaluation of thermodynamically grown (TD-grown) sea ice and snow ice modeled with SMLG_HS against airborne sea-ice thickness measurements over level ice from the same flight. The red square in panels (d) and (h) show the extent of panels (b), (c), (f) and (g). Red color refers to only thermodynamically grown (TD-grown) sea ice, black color indicates the sum of TD-grown sea ice and snow ice, i.e. total sea-ice thickness. In panels (a) and (e), the size of the data point reflects the relative number of airborne measurements in the grid cell. Upper and lower right corners of each panel show the statistics of the corresponding year: Pearson correlation coefficient r, p-value in parenthesis, root-mean-square error (RMSE) and lastly mean bias in parenthesis.

Figure 3. Evaluation of the simulations compared against gridded airborne measurements. Panels (a)–(d) with white background show the modeled snow depth against radar-derived snow depth. The upper panels (a)–(b) show only measurements over level ice and the lower panels (c)–(d) show measurements over all ice types. The left-side panels (a) and (c) show the NASA Operation IceBridge (OIB) flights in 2009–19 and the middle panels (b) and (d) show the AWI IceBird flights in 2017 and 2019. Red color refers to the original SMLG and black color to the new, coupled SMLG_HS. The upper right panel (e) with gray background shows the modeled sea-ice thickness compared against gridded airborne sea-ice thickness measurements over level ice from the AWI IceBird campaigns in 2017 and 2019. Red color refers to only thermodynamically grown (TD-grown) sea ice, black color indicates the sum of TD-grown sea ice and snow ice, i.e. total sea-ice thickness. The size of the data point reflects the relative number of airborne measurements in the grid cell. Upper and lower right corners of each panel show the statistics of the corresponding year: Pearson correlation coefficient r, p-value in parenthesis, root-mean-square error (RMSE) and lastly mean bias in parenthesis.

When we did not distinguish between level and deformed ice and we evaluated SMLG_HS simulations against the total snow depth observations instead (over all ice types), SMLG_HS demonstrated better fit to the snow depth observations from both IceBird and OIB flights (Fig. 2a–d), with reduced RMSEs and biases compared to SMLG (Fig. 3c–d). However, based on a non-parametric Wilcoxon signed-rank test, the differences between the models were not statistically significant. In any case, this indicates that total snow-on-sea-ice amounts given by SMLG_HS are realistic, but they do not account for the sub-grid spatial variations of snow depth (25 km × 25 km). Without considering the sub-grid snow distribution, SMLG_HS overestimated snow depth on level ice resulting in thinner level ice thickness that is more prone to snow-ice formation. The question now becomes: how much snow is removed from the level ice due to snow mass redistribution?

We examined two different approaches to address the question above and to assess the sub-parcel snow mass redistribution: (1) by conducting a modeling sensitivity experiment with a subset of IceBird flights and (2) by performing an analysis of snow volume distribution between level and deformed sea ice based on all IceBird and OIB flight transects. The results of the modeling sensitivity experiment revealed that snow depth on level ice should be reduced by 40% to simulate level ice thickness realistically (Fig. 4) and, at the same time, to maintain snow depth and sea-ice thickness within their respective measurement uncertainties of 0.05 m and 0.12 m for the western Arctic in April 2017. This reduced the mean bias in level sea-ice thickness by 53% from 0.30 m to 0.14 m. The analysis of snow volume distribution along all the OIB and IceBird flight transects in 2009–19 revealed a relationship between the fraction of level ice ( $f_{\mathrm{level}}$) along a sea-ice transect and the fraction of snow volume on level ice ( $f_{V_{\mathrm{s,level}}}$), demonstrating the effect of the sub-parcel snow mass redistribution (Fig. 5). This relationship is linear for fractions of level ice up to 0.5, and it can be given by

(5)\begin{equation} f_{V_{\mathrm{s,level}}} = \left(0.68 \pm 0.05\right) \times f_{\mathrm{level}}, \end{equation}

where ±0.05 represents the 95% confidence interval of the slope.

Figure 4. Results of the sensitivity experiment showing (a) snow depth over all (level and deformed) ice types, (b) snow depth over level ice only, (c) sea-ice thickness over level ice and (d) location of the three IceBird flights (red lines; Table 1) together with the sea-ice type in April 2017 at the time of the flights. The control simulation with unmodified snow depth (SMLG_HS ctrl) is shown as red circles and the simulation with snow depth reduced by 40% (SMLG_HS 0p6) as black circles. The size of the data point reflects the relative number of airborne measurements in the grid cell. While 38% of the total data are from the level ice, the total number of the grid cells (N = 57) is not reduced. Upper and lower right corners of panels (a)–(c) show the statistics of the datasets: the number above is the Pearson correlation coefficient r with p-value in parenthesis, while below are the root-mean-square error and lastly mean bias in parenthesis. OW stands for open water, FYI for first-year ice, SYI for second-year ice (i.e. sea ice that has survived one melt season) and MYI for multi-year ice (i.e. sea ice that has survived at least two or more melt seasons).

Figure 5. The relationship between the fraction of level ice and the fraction of snow volume on level ice demonstrating the effect of snow mass redistribution (gray hatching). The NASA OIB survey flights are marked with black circles and their linear fit with a black dashed line, whereas the AWI IceBird ones are shown with red crosses and a red dashed line. The solid black line shows the linear fit of all airborne data and the gray shading is its 95% confidence interval. The blue stars show the corresponding end-of-winter values in March–April 2020 from the MOSAiC expedition ground-based transect by Itkin and others (2023).

4. Discussion

We performed a modeling study to investigate snow sinks on level ice in the western Arctic. Specifically, we examined snow sinks caused by snow and ice interactions, such as snow-ice formation and snow mass redistribution. Snow loss in leads was not considered in this study, because observations from the MOSAiC expedition demonstrated that this is likely an insignificant snow sink in winter, due to quick refreezing of the leads (Clemens-Sewall and others, 2023). We coupled SMLG snow depth and density evolution with HIGHTSI thermodynamic sea-ice and snow-ice growth to create SMLG_HS. Being in fact a 1-D model, SMLG_HS should be considered for level ice only. It does not account for dynamic ice thickening, nor for sub-parcel snow mass redistribution processes, i.e. the preference of snow to accumulate over ice deformations (Liston and others, 2018). It also assumes that negative freeboard will lead to snow-ice formation. Therefore, it is expected to overestimate snow depth on level ice. Being a very effective insulator, this additional snow decelerates level ice growth, resulting in underestimation of level ice thickness and overestimation of snow-ice thickness. This hypothesis was confirmed when we compared SMLG_HS simulations to airborne observations of snow depth over level ice and level ice thicknesses. SMLG_HS matched the overall snow depth observations from airborne radars better compared to SMLG, with reduced RMSEs and biases. However, based on a statistical test (non-parametric Wilcoxon signed-rank test), the differences between the snow depths in the two models were not statistically significant.

AWI IceBird data are ideal for evaluating SMLG_HS, because they offer simultaneous snow depth and sea-ice thickness observations over hundreds of kilometers of transects in high resolution, with a possibility to examine level and deformed ice conditions separately. However, IceBird campaigns that provide a concrete dataset of both snow depth and sea-ice thickness observations are limited to the western Arctic and in April 2017 and 2019 only. In 2019, IceBird flew over multi-year ice that was heavily deformed, resulting in small fractions of level ice along the flight tracks. The limited level ice observations would impose a risk of unreliable conclusions; therefore, we focused our analysis on flights with the largest level ice fraction in 2017. Moreover, the monitored region was occasionally close to the coast where parcel trajectory data are unavailable, rendering these regions outside the simulation domain and the sensitivity experiment.

The modeling sensitivity experiment showed that reducing snow depth by 40% produced the best agreement between snow depth (on level ice) and level ice thickness in the western Arctic in April 2017 and reduced the mean bias in sea-ice thickness by 53%. The analysis of the snow volume distribution between level and deformed sea ice using observations from the IceBird and OIB transects was in good agreement with the model results when considering their respective limitations. The sensitivity experiment relies on a 1-D thermodynamic model that does not account for lateral conduction of heat, a factor that becomes significant when snow depth varies spatially (Clemens-Sewall and others, 2024; Zampieri and others, 2024). Regarding the airborne approach, it is not possible to account for snow sinks in snow-ice formation. Omitting snow-ice formation, that mostly occurs over level ice, would result in underestimation of the snow mass redistribution.

We argue that the snow depth decrease on level ice represents the sub-parcel snow mass redistribution process; however, this mechanism is not yet fully understood. The deformation rate of a sea-ice floe, together with the atmospheric conditions (e.g. wind, warm intrusions) and the properties of snow cover (density, wetness, sintering level and snow-surface shear strength) are expected to affect the snow redistribution, i.e. the amount of snow removed from the level to deformed ice. Ice and snow conditions are not uniform across the Arctic Ocean, but they vary regionally and temporally. Therefore, a 40% reduction of snow depth on level ice is empirical and more data are needed across the Arctic and the different seasons to study the spatiotemporal variability of snow mass redistribution. In another, yet more local example by Itkin and others (2023), data from the MOSAiC expedition indicated that 31% of level ice contained only 18% of the snow volume at the end of spring (see the blue stars in Fig. 5). In the Surface Heat Budget of the Arctic Ocean (SHEBA) study in 1997–98, snowdrifts associated with ridges occupied between 3% and 6% of the total study area. The drift sections had mean depths that were on average 30% higher than the surrounding snow (Sturm and others, 2002).

Although the snow depth reduction suggested by the sensitivity experiment cannot be generalized across the entire Arctic and across different years, as an illustrative attempt, we compared snow-ice formation results from the SMLG_HS simulation spanning the years 2007–21, with and without a 40% decrease in snow depth. The 14 year-average snow-ice thickness on the day of maximum snow-on-sea-ice volume is shown in Figure 6. Even with a 40% decrease in snow depth, snow ice still has the potential to form and is characterized by strong seasonal and regional variations. However, 40% less snow on level ice would greatly limit snow-ice formation in the central and western Arctic. This process would be primarily restricted to the Atlantic sector of the Arctic, particularly along Greenland’s east coast and north of Svalbard underneath the North Atlantic storm track, where the N-ICE2015 campaign was conducted. Snow-ice formation has also been observed with autonomous sea-ice mass-balance buoys similar to Provost and others (2017) (Text S2 and Fig. S1 in the supplementary material) and in fully coupled climate models (Webster and others, 2021) in these regions in the contemporary period. Understanding the importance of sub-parcel snow mass redistribution will guide the development of necessary modeling tools that capture snow sinks properly.

Figure 6. Snow-ice thickness, 14 year-average over the day of maximum snow-on-sea-ice volume in 2007–21, from (a) the control run (SMLG_HS ctrl), (b) the run with snow depth reduced by 40% (SMLG_HS 0p6) and (c) the difference between the two simulations (reduced minus control).

5. Conclusions

We showed that a 1-D sea-ice model and snow thermodynamic model approach would overestimate snow sink in snow-ice formation. Even though the total snow depth (over both deformed and level ice) matched well with both OIB and IceBird observations, not accounting for snow redistribution from level to deformed ice resulted in overestimation of snow depth over level ice. As expected, this additional snow decelerated thermodynamic ice growth in the model, resulting in thinner level ice that is more prone to snow-ice formation. Based on the evaluation of our simulations against IceBird data in April 2017, fitting both snow on level ice and level ice thickness simulations to the IceBird observations, snow depth in SMLG_HS should be reduced by 40%. We argue that in our 2017 case study in the western Arctic, 40% reduction in snow depth over level ice represented the sub-parcel snow mass redistribution process. Based on the analysis of >100 airborne survey flights spanning a full decade, the fraction of snow volume on level ice in spring is linearly related to the level ice fraction, and it is given by $f_{V_{\mathrm{s,level}}} = \left(0.68 \pm 0.05\right) \times f_{\mathrm{level}}$. This linear relationship indicates that the amount of snow remaining on level sea ice in the end of winter is proportional to the amount of ice deformation.

When snow models do not account for snow sinks caused by snow and sea-ice interactions, such as snow-ice formation or sub-parcel snow mass redistribution processes, they overestimate snow depth on level ice. Uneven snow-on-sea-ice load within a sub-grid area will result in biases in altimetry retrievals of sea-ice thickness by overestimating level ice and underestimating deformed ice thickness. Regarding sea-ice modeling applications, spatial variability in snow depth will impact sea-ice thermodynamic growth in winter, affecting both vertical and horizontal heat fluxes, and will influence melt pond formation in summer (Thielke and others, 2023). Therefore, snow-on-sea-ice reconstructions should be used with caution depending on the application requirements. This study emphasizes the need to account for sub-grid scale heterogeneity in snow and sea-ice interactions to improve the representation of snow in remote sensing and model studies. It also highlights the crucial need for additional independent but simultaneous observations of snow depth and sea-ice thickness, together with information on snow properties, to understand the mechanism behind snow mass changes due to coupled physical processes.