Snow on sea ice plays a critical role in the polar oceans’ energy balance, but also in satellite retrievals of sea-ice thickness among other variables. The density of snow on sea ice evolves over the winter season, generally increasing as grains become rounder and the snowpack settles due to the effect of overburden. It is, therefore, desirable to form a simple equation for the snow density as a function of the time-of-year. In order to investigate the role of snow in radar-derived estimates of Arctic sea-ice thickness, such an equation was put forward by Mallett and others (Reference Mallett, Lawrence, Stroeve, Landy and Tsamados2020, henceforth M20):

(1)\begin{equation} \qquad \qquad \qquad \qquad \qquad \qquad \rho_\mathrm{s} = 6.5t_\mathrm{m} + 274.51 \end{equation}

where $\rho_\mathrm{s}$ is the snow density in $\mathrm{kg\ m}^{-3}$, and $t_\mathrm{m}$ is the number of months since October. The equation has now been used in several publications (e.g. Dong and others, Reference Dong, Shi, Lin and Zeng2022, Reference Dong, Shi, Lin, Jia, Zeng and Wu2023; Jiang and others, Reference Jiang, Zhong, Xu and Jia2023; Shi and others, Reference Shi2023; Sievers and others, Reference Sievers, Rasmussen and Stenseng2023; Chen and others, Reference Chen2024; Fredensborg Hansen and others, Reference Fredensborg Hansen2024).

Equation (1) was computed as follows: a large dataset of snow depth and snow water equivalent (SWE) was compiled from in situ measurements at Soviet North Pole (NP) drifting stations by Warren and others (Reference Warren, Rigor, Untersteiner, Radionov, Bryazgin, Aleksandrov and Colony1999), and monthly quadratic fits were published for both variables. Following common practice in radar altimetry processing chains, M20 divided the quadratic fits for SWE by those for depth to produce spatial distributions for snow density. The spatial average of these density distributions in a subdomain of the Arctic Ocean was then computed, producing one mean snow density value for each winter month. These values were then regressed against the month number to generate Eqn (1) of this manuscript. The above method has several drawbacks; their impact and remediation are the subject of this communication.

The first limitation of the method described above concerns the original quadratic fits for SWE and depth themselves, the parameters of which were published by Warren and others Reference Warren, Rigor, Untersteiner, Radionov, Bryazgin, Aleksandrov and Colony(1999). In some months, the quadratic fits can produce negative values in the marginal seas of the Arctic and are not inherently ‘snow conserving’ (i.e. the mean value in the Arctic Ocean is not inherently the mean value of the underlying values, particularly since the spatial definition of the Arctic Ocean is not well-defined). Furthermore, it is sub-optimal to compute monthly spatial distributions for density by dividing those for SWE by those for depth: it would be better to compute the density distributions directly from the density measurements and their positions in the month concerned.

Further drawbacks exist in the averaging and regression process underpinning Eqn (1): the area over which M20 averaged the density distributions in each month goes beyond the area sampled by the NP station data. For example, the area considered by M20 includes the Laptev Sea, from which stations rarely collected data. It was also only performed in the months of October–April, when the source data from NP stations would potentially allow a function to apply beyond those months. Finally, $t_\mathrm{m}$ in Eqn (1) represents the integer number of months since October, indicating that the formula is not weighted for the variable lengths of the winter months. In a sense, it is linear in month number and thus not strictly linear in time.

All the methodological issues described above can be reduced (and some resolved), by directly regressing the mean densities calculated from the original transect data against the time-of-year at which they were generated. These data can be downloaded from the Joint US-Russian Sea Ice Atlas (Environmental Working Group, Reference Tanis and Smolyanitsky2000). Measurements were taken in bulk, by weighing a cylinder of $50\ \mathrm{cm}^{2}$ in cross section that had been pushed vertically down to the snow–ice interface (Colony and others, Reference Colony, Radionov and Tanis1998). After some data cleaning (see below), this regression yields:

(2)\begin{equation} \qquad \qquad \qquad \qquad \qquad \qquad \rho_\mathrm{s} = 0.35t_\mathrm{A_1} + 239.78 \end{equation}

where $t_\mathrm{A_1}$ represents the number of days since August 1st, and $\rho_\mathrm{s}$ remains the snow density in $\mathrm{kg\ m}^{-3}$ as in Eqn (1). Five out of 578 data points have been removed for quality-control reasons. These were recorded in the months of July and August: four of them are $ \gt 500\,\mathrm{kg\ m}^{-3}$ and one of them is $25\ \mathrm{kg\,m}^{-3}$ (this is likely a measurement error). These extreme values exist near the August 1st break-point of the analysis, and their inclusion makes the slope of the regression highly sensitive to the choice of this date. Because of their removal, it is inadvisable to generate snow densities from Eqn (2) in July and August. Despite this, it is clear that Eqn (2) can sensibly be used to produce values outside of the ‘cold season’ considered by the M20 calculation, for instance in September, May and June. Individual transect mean values in Fig. 1 are scattered about the regression line (Eqn (2)) with a root-mean-square error of $34.9\ \mathrm{kg\ m}^{-3}$. This is the typical error that a user of the function should expect in an ice environment similar to that from which the NP data were collected.

Figure 1. Transect-mean snow densities (n = 573; black scatter), with the M20 values shown as red lines. Linear regression through the scatter points shown in dark blue. Where possible, NP station transects were performed at 10-day intervals on the 10th, 20th and 30th of each month, generating a periodic distribution of scatter along the time-axis. Dates are shown on upper x-axis for non-leap-years.

Figure 1 also makes clear that the new regression slope is not very different from the M20 function in a quantitative sense. Density calculations in the publications cited above using M20 can therefore be trusted. So why make a new one? The first reason is that the new, simpler, more robust methodology can be better trusted in future to represent the underlying data, and in more months of the year. In addition, the new function also takes a more continuous input of days since August 1st rather than the month number, aiding its utility as described above.

This new densification function retains some key limitations. It still relies on data collected by Soviet NP drifting stations that operated on multiyear ice and overwhelmingly in the Central Arctic, East Siberian and Chukchi seas (see Figure 2 of Mallett and others, Reference Mallett2021, for trajectories of stations contributing measurements to this analysis). Snow in the multiyear ice environment may well have a different densification rate to that in the first-year ice environment due to its relative lack of salinity and the rougher underlying ice. Relatedly, the high latitude of the measurements means that the densification rate in Eqn (2) may not reflect that of lower latitudes where periods of diurnal cycling are more protracted and temperatures are often higher.