Retrieval and parametrisation of sea-ice bulk density from airborne multi-sensor measurements

. Knowledge of sea-ice thickness and volume depends on freeboard observations from satellite altimeters and in turn on information of snow mass and sea-ice density required for the freeboard-to-thickness conversion. These parameters, especially sea-ice density, are usually based on climatologies constructed from in situ observations made in the 1980s and earlier while contemporary and representative measurements are lacking. Our aim with this paper is to derive updated sea-ice bulk density estimates suitable for the present Arctic sea-ice cover and a range of ice types to reduce uncertainties in sea-ice 5 thickness remote sensing. Our sea-ice density measurements are based on over 3000 km of high-resolution collocated airborne sea-ice and snow thickness and freeboard measurements in the western Arctic Ocean in 2017 and 2019. Sea-ice bulk density is derived assuming isostatic equilibrium for different ice types. Our results show higher average bulk densities for both ﬁrst-year ice (FYI) and especially multi-year ice (MYI) compared to previous studies. In addition, we ﬁnd a small difference between deformed and possibly unconsolidated FYI and younger MYI. We ﬁnd a negative-exponential relationship between sea-ice bulk 10 density and sea-ice freeboard and apply this parametrisation to one winter of monthly gridded CryoSat-2 sea-ice freeboard data. We discuss the suitability and the impact of the derived FYI and MYI bulk densities for sea-ice thickness retrievals and the uncertainty related to the indirect method of measuring sea-ice bulk density. The results suggest that retrieval algorithms be adapted to changes in sea-ice density and highlight the need of future studies to evaluate the impact of density parametrisation on the full sea-ice thickness data record.

waterline ice is saturated by sea water and has a density of 900-940 kg m −3 less dependent on its age (Timco and Frederking, 1996;Timco and Weeks, 2010;Pustogvar and Kulyakhtin, 2016). Especially with satellite altimetry applications in mind, 60 Alexandrov et al. (2010) derived a density of 916.7 ± 35.7 kg m −3 for FYI using the drill-hole data set of airborne Sever expeditions concentrated in the shelf seas of the Eurasian Russian Arctic in the 1980s and a density of 882 ± 23 kg m −3 for MYI as a weighted average of the layers above and below the waterline using values from literature. The values by Alexandrov et al. (2010) (hereafter A10) are the most commonly used in sea-ice thickness retrieval algorithms (Sallila et al., 2019). It has to be noted that majority of the density measurements originate from Arctic sites, which is the case also in our study. Different 65 properties and processes of Antarctic sea ice could lead to different density values.
Accurate and representative measurements of sea-ice density using traditional techniques are temporally and spatially limited. Most of them require coring or cutting out pieces of ice, such as the mass/volume, displacement (submersion), or specific gravity techniques, making them susceptible for inaccuracies through brine drainage and imprecise volume of the samples (Timco and Frederking, 1996). This can be avoided by carefully recording sea-ice thickness and freeboard, e.g., with drill-hole 70 measurements, in addition to the snow depth atop and calculating the sea-ice bulk density assuming isostatic equilibrium and densities of snow and sea water. However, significant error may be introduced locally where sea ice is not isostatically compensated due to lateral stresses, e.g., close to pressure ridges (Timco and Frederking, 1996). Previous parametrisations of sea-ice density include the effective freeboard approach (snow depth converted to ice thickness using their density ratio) by Ackley et al. (1976) using drill-hole measurements from 400 m of profile lines on MYI in the Beaufort Sea and a parametrisation 75 based on ice floe thickness by Kovacs (1997) utilising 17 FYI and 4 MYI sea-ice cores from the Beaufort Sea. Neither of these parametrisations have been widely used. Moreover, the multi-phase nature of sea ice is an ongoing challenge for modelling approaches (Hunke et al., 2011). There is a definite need for evaluating sea-ice density because there is no density climatology available representing the current state of sea ice, nor is it possible to observe density by satellites from space.
Simultaneous, collocated, and preferably single-platform measurements of the key parameters of the entire sea-ice-snow 80 layer covering a wide range of ice types and conditions on regional scales are required to decrease the uncertainties related to the conversion of freeboard to sea-ice thickness. Since 2017, a unique sensor configuration on the Alfred Wegener Institute's (AWI) IceBird winter campaigns combines airborne laser, radar, and electromagnetic induction sounding instruments making it now possible to measure them on a single platform. In this paper, we present high-resolution data of simultaneous airborne sea-ice thickness, freeboard, and snow depth over late-winter Arctic sea ice from the AWI IceBird campaigns in 2017 and 85 2019. Observing the locations of the air-snow, snow-ice, and ice-water interfaces in the sea-ice system along survey tracks allows us to estimate sea-ice bulk density that also serves as a consistency check between the sea-ice thickness, freeboard, and snow depth measurements. Our objective is to derive an updated parametrisation of sea-ice bulk density based on a single variable observable from space. The resulting parametrisation should be suitable for the present Arctic sea-ice cover including the densities of deformed sea ice which, if unconsolidated, can deviate even more strongly from the density of solid ice. 2 Data and methods

Aircraft campaigns
The AWI IceBird program (see reference list for a link to webpage) is a series of airborne campaigns carried out using the institute's two Basler BT-67 research aircraft Polar5 and 6 (Alfred-Wegener-Institut Helmholtz-Zentrum für Polar-und Meeresforschung, 2016) to measure Arctic sea ice and its change since 2009. The campaigns operate from airfields extending from 95 Longyearbyen in Svalbard to Utqiagvik (Barrow) in Alaska and coincide closely with the Arctic sea ice summer minimum (August) and winter maximum (April). The primary scientific instrumentation on the aircraft includes an electromagnetic (EM) induction sounding instrument (EM-Bird) to measure total (i.e., ice+snow) thickness, an airborne laser scanner (ALS) for surface topography and freeboard measurements, a microwave radar to measure snow depth, and an infrared radiation pyrometer to record surface temperature (Fig. 1). We describe each instrument in the following sections. The low altitude of 100 200 ft (≈ 60 m) and slow speed of 110 kn (≈ 60 m s −1 ) during the nominal surveys are beneficial for high-resolution data acquisition.
In this study, we used data collected during the IceBird winter campaigns in early April of 2017 and 2019 (Table 1; Fig. 2) utilising the unique data set of simultaneous total thickness, snow freeboard, and snow depth measurements. From 2017, we used measurements from four survey flights that took place over the Beaufort and Chukchi Seas as part of the Polar Airborne 105 Measurements and Arctic Regional Climate Model Simulation Project (PAMARCMiP; Haas et al., 2010;Herber et al., 2012).
From 2019, we included five survey flights that covered regions in the Lincoln Sea and the Arctic Ocean in addition to an overlap with the measurements in the Beaufort Sea in 2017.

Sea-ice thickness
We measured total (ice+snow) thickness of sea ice (h tot ) using the towed EM induction sounding instrument, the EM-Bird, 110 suspended below the aircraft 10-20 m above the sea ice surface (Haas et al., 2009) as illustrated in Fig. 1. The EM-Bird utilises the contrast between resistive snow and ice layers and conductive sea water by transmitting an EM field that induces eddy currents only in the latter (Kovacs and Morey, 1991;Haas et al., 1997). The EM-Bird measures the phase and amplitude of a secondary EM field induced by those eddy currents in relation to the primary field. The phase and amplitude of the secondary field depend on the distance between the instrument and the ice-water interface and decrease negative-exponentially with 115 increasing distance. Subtracting the instrument height above the surface, measured by an integrated laser altimeter, from the distance to the sea water gives the total thickness (Haas et al., 2009(Haas et al., , 2021. Sea-ice thickness (h i ) was derived by subtracting snow depth from total thickness (see Sect. 2.4). The EM-Bird sampling rate was 10 Hz, which translated to 5-6 m point spacing at the nominal survey speed. Approximately every 15-20 minutes brief ascents to more than 100 m were carried out to monitor the sensor drift during post-processing (Haas et al., 2009). Comparison to drill-hole measurements over level ice 120 have indicated an accuracy of 0.1 m, whereas ridge peak thicknesses are generally underestimated by up to 50 % as a result of mass-conserving averaging effects within the approximately 40 m diameter footprint of the instrument (Pfaffling et al., 2007;Haas et al., 2009). For our analysis, we disregarded measurements of total thickness (i) less than the instrument accuracy of 0.1 m, (ii) where total thickness was less than the mean snow freeboard or snow depth, and (iii) where the surface temperature was above −5 • C within the footprint of the instrument to avoid open or newly frozen leads with total thickness below the 125 accuracy of the EM-Bird (see Sect. 2.4 and 2.5.1).

Freeboard
The near-infrared (1064 nm), line-scanning Riegl VQ-580 airborne laser scanner (ALS) measured ellipsoidal elevations of ice surfaces with a 60 • field of view resulting in a swath width approximately equal to the aircraft's altitude above ground (nominally ≈ 60 m). We obtained freeboard from the ALS data by subtracting the local sea-surface height from the ice surface of the sea-ice cover or newly frozen leads with negligible freeboard and we manually selected the corresponding elevations.
The spacing of leads was diverse and depended on the ice regime varying between less than 10 km (FYI) and potentially more than 30 km (densely packed MYI). We subtracted the mean sea surface (DTU15 MSS; Andersen et al., 2016) from the surface elevations to remove large scale variations and reduce interpolation errors before applying a spline interpolation between 135 the sea-surface height tie points along the flight track. Data before the first lead detection and after the last lead detection during a survey flight were discarded to avoid extrapolation errors. Subtracting the interpolated sea-surface height from the ice elevations results in snow freeboard (h f s ), as the elevation measurement of the ALS always includes the snow layer ( Fig. 1).
Sea-ice freeboard (h f i ), i.e., the location of the snow-ice interface in relation to the local sea level, was derived by subtracting snow depth from snow freeboard (see Sect. 2.4). We then interpolated the obtained point cloud data of snow freeboard onto 140 a regular grid with a 0.25 m resolution. Freeboard uncertainties are dominated by the accuracy of the interpolation of the instantaneous sea surface anomaly that depends on the abundance of leads. Therefore, especially regions with packed MYI and low lead density are associated with high uncertainties (Ricker et al., 2016) and we manually masked out such areas.
Supervising the along-track interpolation, we estimated an overall uncertainty of 0.1 m. For each total thickness measurement, we calculated the corresponding mean snow freeboard within the EM-Bird footprint (Fig. 3).

Snow depth
We used an ultra-wideband ( .

155
To calibrate the raw Snow Radar data, we used a workflow described in Jutila et al. (2022) including coherent noise removal and system impulse deconvolution. Using an open-source python package pySnowRadar (King et al., 2020a), we detected airsnow and snow-ice interfaces with the algorithm by Jutila et al. (2022) that is based on a pulse peakiness approach by Ricker et al. (2014) used in satellite altimetry. Taking into account the decreased wave propagation speed in snow by assuming a snow density, the distance between the identified interfaces determined the snow depth estimate. In post-processing, we filtered 160 out values which were acquired during the EM-Bird calibration manoeuvres with a simple altitude threshold of 100 m and when the absolute roll or pitch of the aircraft exceeded 5 • . Additionally, for each snow depth estimate we calculated a surface topography estimate (h topo ) as the difference between the 95th and 5th percentile of the ALS surface elevation data within the radar footprint to disregard potentially erroneous snow depth estimates over heavily deformed sea ice using a threshold value of 0.5 m (Jutila et al., 2022).

165
A validation exercise over level, landfast FYI yielded a negligible mean bias of 0.86 cm, which was below the radar resolution and within the accuracy of the ground truth data, and a root-mean-square error (RMSE) of 6.9 cm for the radar-derived snow depth estimates (Jutila et al., 2022). For each total thickness measurement, we calculated the corresponding mean snow depth requiring at least five valid snow depth estimates within the EM-Bird footprint (corresponds to approximately 50 % of the values; Fig. 3). The averaging reduced the error by a factor of √ N , where N is the number of averaged estimates. We 170 disregarded averaged snow depth values where the surface temperature within the EM-Bird footprint was above −5 • C and where the mean snow freeboard was less than the mean snow depth (negative sea-ice freeboard) to avoid potentially erroneous snow depth retrievals due to changes in the dielectric properties of snow induced by liquid water (Barber et al., 1995;Kurtz and Farrell, 2011;Kurtz et al., 2013;Rösel et al., 2021). However, it has to be noted that the studies of Kurtz and Farrell (2011), Kurtz et al. (2013), and Rösel et al. (2021) use different radar versions and retrieval algorithms compared to this study, which 175 may limit the direct applicability of their results.

Surface temperature
Surface temperature was acquired using the Heitronics infrared radiation pyrometer KT19.85II that recorded the 9.6-11.5 µm spectral band response of the surface with a sampling rate of 50 Hz resulting in an approximately 1 m sample spacing and 3.1 m 180 diameter footprint at the nominal survey speed and altitude. The manufacturer reported an accuracy of ±0.5 • C + 0.7 % of the temperature difference between the target and the instrument housing. We used the measurements to filter out total thickness and snow depth measurements where the surface temperature was above −5 • C (see Sect. 2.2 & 2.4).

Sea-ice type
Information of sea-ice type is required for accurate classification the sampled ice. However, no remote sensing product or 185 modelling output is able to match the resolution of the airborne survey data. Therefore, we used a custom sea-ice classification scheme.
We started with identifying level and deformed ice following the approach of Rabenstein et al. (2010). The filter is based on the observation that level ice is mostly flat and extends over long distances. We identified data points that fulfilled those characteristics using two criteria. First, we calculated the along-track total thickness gradient using a three-point Lagrangian 190 interpolator. We applied a threshold gradient of 4 cm within an along-track distance of 1 m, below which the ice was classified as level following Rabenstein et al. (2010). Second, this condition must be met continuously for at least 100 m of the profile length. Choosing the value of 100 m, which represents approximately twice the footprint size of the EM-bird, makes sure that the conditions were met over two completely independent EM total thickness measurements. If these criteria were not fulfilled, the ice was deemed deformed.
We then chose the nearest neighbour data point from the coinciding EASE-Grid Sea Ice Age (Version 4) product from the National Snow & Ice Data Center (NSIDC; Tschudi et al., 2019) providing weekly sea-ice age estimates in 12.5 km resolution.
Where the NSIDC Sea Ice Age data were not available (landfast ice and close to coasts), we manually assigned the ice type to FYI or MYI (old ice, including SYI) according to the Canadian Ice Service (CIS) regional and weekly ice charts (Canadian Ice Service, 2009). Finally, we defined the ice type as (i) first-year ice (FYI), if the ice was younger than 1 year according 200 to NSIDC/CIS or the observed ice thickness was below 2 m regardless of its age; (ii) second-year ice (SYI), if the ice had a thickness of 2 m or more and its age was 1-2 years; and (iii) multi-year ice (MYI), if the ice had a thickness of 2 m or more and was at least 2 years old. To account for the spatial and temporal limitations of the NSIDC Sea Ice Age product and the drift of sea ice, we adjusted the ice type classification from FYI to MYI for any ice that indicated an age less than 1 year but was level and thicker than 2 m or, after along-track averaging over a length scale (see Sect. 2.6), the lower quartile (25th quantile) 205 of the averaged ice thickness values within the chosen length scale was above 2 m.
To support the analysis of the sampled sea ice and to evaluate the indicated sea-ice age in the NSIDC product, we investigated the sea-ice age, pathways, and origin using the Lagrangian drift analysis system ICETrack (Krumpen, 2018;Krumpen et al., 2020). We split the surveyed sea ice into 25 km along-track segments and tracked them backwards in time in daily increments utilising a publicly available low-resolution satellite sea ice motion product from the Ocean and Sea Ice Satellite Application  et al., 2007), along the backward trajectory dropped below 25 %, we assumed that the ice was formed in that specific location.
To quantify uncertainties of sea-ice trajectories, Krumpen et al. (2019) reconstructed the pathways of 57 drifting buoys. The authors showed that the deviation between actual and virtual tracks was rather small, 36 ± 20 km after 200 days, and considered 215 to be in an acceptable range.

Sea-ice bulk density
Simultaneous measurements of sea-ice thickness, snow depth, and freeboard enable us to calculate sea-ice bulk density using the so-called "freeboard and ice thickness technique" (Timco and Frederking, 1996). Archimedes' principle dictates where ρ is density for ice, snow, and sea water denoted with subscripts i, s, and w, respectively. The terms can be rearranged to solve for ρ i : To solve Eq.
(3), we need to assume values only for the densities of sea water and snow, but their impact on the uncertainty of sea-ice bulk density is small (see Eq. (4) below). Here we took sea-water density and its uncertainty according to Wadhams et al. (1992) as ρ w = 1024 kg m −3 and σ ρw = 0.5 kg m −3 , respectively. For snow density in April, when measurements were carried out, we chose ρ s = 300 kg m −3 following Warren et al. (1999) and for the respective uncertainty σ ρs = 34 kg m −3

230
from King et al. (2020b). These values and the uncertainties of the measured variables are summarised in Table 2. Assuming that the individual uncertainties are uncorrelated, we can derive uncertainty for sea-ice bulk density (σ ρi ) using Gaussian error propagation: Following the uncertainty source analysis of Giles et al. (2007) and using the values summarised in Table 2, the largest contributors to the uncertainty of sea-ice bulk density are, in descending order of magnitude, snow freeboard, snow depth, and total thickness. We disregarded density values with uncertainty exceeding 100 kg m −3 from further analysis.
While the assumption of isostatic equilibrium may not necessarily be valid locally, e.g., close to pressure ridges, it holds true 240 when averaging over a sufficient length scale. We varied the averaging length in 10 m increments and found that the mean bulk densities and the standard deviations for different ice types of the surveys did not change significantly beyond a length scale of about 200 m, if measurements of each ice type were abundant. Here, we computed sea-ice bulk density estimates at two length scales representing the sensor resolution of the CryoSat-2 satellite as well as the typical resolution of gridded sea-ice thickness data. In the first case, we approximated the scale of full-resolution altimeter footprint by the diameter of a circle with the same 245 area as the 300×1650 m pulse-Doppler-limited footprint of CryoSat-2 in the synthetic aperture radar (SAR) acquisition mode.
The diameter of that circle is equal to approximately 800 m. In the second case, we chose the typical satellite product grid cell size of 25 km. We assumed that the sea-ice and snow layers are in isostatic equilibrium at both scales. We calculated an along-track weighted average using the squares of individual uncertainty values as weights: where k is the index and N equals the number of values within the length scale to be averaged. The resulting uncertainty was determined with We calculated inverse-uncertainty weighted averages also for snow depth since its uncertainty varied spatially. For all other variables we calculated an arithmetic mean.

Results
This study included a total of 3410 airborne survey kilometres split approximately equally between the years 2017 and 2019 (Table 1). The abundance of different sea-ice types varied between the years and the individual surveys (see the percentages in Table 3). Surveyed sea ice in 2017 was solely FYI. This was a result of the ice-free conditions in the Beaufort and Chukchi Seas in the previous summer, additionally influenced by the collapse of the semipermanent Beaufort high pressure system,     Table 3 and shown in Fig. 6. On average, FYI bulk density was higher than the value derived by Alexandrov et al. (2010) (A10), also for individual surveys. FYI bulk density in 2017 was slightly higher than in 2019, and combined they resulted in an average density of 928.5 ± 16.4 kg m −3 . SYI and MYI bulk densities differed only a little from each other, but were 23-30 kg m −3 lower than FYI bulk density. Similar to FYI, the bulk densities of old ice were in the upper 280 range of or even beyond A10. Figure 7 shows the spatial distribution of the derived sea-ice bulk density when averaged over a typical satellite product grid cell size of 25 km. As expected, lower density values were generally encountered with increasing ice age ( Fig. 7e and f). The lowest sea-ice bulk density values were located 50-150 km northwest off the edge of the landfast Table 3. Summary of the 800 m along-track averaged sea-ice bulk density (inverse-uncertainty weighted mean ± one standard deviation) according to survey and ice type. The mean values correspond to the red crosses in Fig. 6. The percentage of each sea-ice type encountered on each survey is given in parentheses. ice and the coast of Ellesmere Island and Nansen Sound (Fig. 7e). Moreover, there was noticeable along-track variability of bulk density also within a single ice type, even FYI where also the highest bulk density values were observed.

Parametrisation of sea-ice bulk density
We explored the possibility to parametrise sea-ice bulk density using one of the measured sea-ice parameters. Out of the full parameter space, sea-ice freeboard showed the best correlation with the estimated sea-ice bulk density, r = −0.62 (p 0.001), indicating a significant, linear anti-correlation as expected according to Eq. (2) (ρ i ∝ −h f i ). However, the dependence on other sea-ice properties and the fact that keels of ridges with high freeboard contain voids filled with sea water introduced 290 non-linearity to this relationship. To avoid underestimating the density values near the lower and upper ends of the observed freeboard range and to increase the goodness of fit, we therefore fitted an exponential function to the full, along-track averaged sea-ice bulk density data set (N = 2246, Fig. 8). Least-squares fitting yielded a relationship of where

Discussion
We measure sea-ice bulk density indirectly based on (i) the direct measurements of sea-ice thickness, snow depth, and freeboard and (ii) the isostatic balance between the masses of snow, sea ice, and displaced sea water only assuming the densities of the snow layer and sea water. However, we measure the thickness of the sea ice layer and not its mass. Instead, we use sea-ice bulk density to relate sea-ice thickness to the displaced mass of sea water inherently by assuming a constant density of the 305 entire sea-ice layer. However, in reality, the material composition of sea ice is not constant throughout the vertical column which complicates the attribution of bulk sea ice density values to physical properties such as porosity. Above the waterline, the density is lower than that of pure ice due to air incorporated in the pore spaces. This feature is more pronounced in MYI.
Below the waterline, brine and sea water saturate the sea ice and increase the density above the pure ice density. Despite the indirect measurement method, we are able to detect a difference in FYI bulk density between 2017 and 2019 that can be linked 310 to the high sea-ice deformation in 2017. The effect of deformed and unconsolidated sea-ice is often overlooked and needs the attention of the scientific community. Dedicated sea-ice porosity studies and extensive field measurement programs, such as the recent Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC), will be able to shed more light on their effect and development.

315
Compared to A10, the average sea-ice bulk density estimates derived in this paper are larger by 11.8 kg m −3 (≈ 1.3 %) for FYI and by 20.4 kg m −3 (≈ 2.3 %) for MYI but still within the A10 uncertainties, albeit close to the upper limit. In general, our ice-type averaged bulk density estimates fall within the range of previous studies (Timco and Frederking, 1996). Reasons for the comparably high estimates are twofold. First, the A10 FYI density is representative only to level, undeformed ice whereas our estimate includes also deformed FYI. Alexandrov et al. (2010) used drill-hole measurements that were carried out   Timco and Frederking (1996) instead, we find a MYI density of 909 ± 28 kg m −3 that is closer to our estimate.

Uncertainties and limitations of the derived sea-ice bulk density 330
The effect of the uncertainties in the measured parameters and assumed sea-water and snow density values on the sea-ice bulk density was studied using Gaussian error propagation in Eq. (4). Single point measurements typically resulted in seaice bulk density uncertainties of approximately 70 kg m −3 and 35 kg m −3 for FYI and MYI, respectively. Other sources of error arise from the different length scales illustrated in Figs. 3 and 5. We are not able to resolve snow depth fully within the total thickness measurement given the comparably large footprint of the EM-Bird but we calculate it as an average of 335 the snow depth measurements along a chord of the circular EM-Bird footprint. Due to the cross-track movement of the EM-Bird under the aircraft, the ground locations and the number of the snow depth measurements within the total thickness measurement vary slightly along the survey track but generally remain at eight to ten measurements close to the centre line.
To ensure representative snow depth estimates, we require at least five valid snow depth measurements for each total thickness measurement, which translates into at least 50 % coverage along the chord. Errors may occur locally, e.g., at cross-track 340 transition from a sea-ice floe to young ice in a newly refrozen lead as in the leftmost measurement points in Fig. 3, but we assume them to occur randomly and not cause systematic bias. However, uncertainties are reduced when averaged along-track over a length scale. Using the 800 m length scale, the resulting sea-ice bulk density uncertainty is generally less than 10 kg m −3 (Eq. (6)) but remain the highest for thin ice and low sea-ice freeboard where the relative uncertainties of the input parameters are the largest (see the size of the scatter points in Fig. 8). In turn, averaging over a length scale simplifies the natural variability 345 of sea ice. Figure 5 shows how a single 25 km satellite grid cell can contain already several sea-ice types: level, deformed, FYI, SYI, and MYI. However, assigning the sea-ice types is limited by the spatial (12.5 km) and temporal (weekly) resolution of the NSIDC sea-ice age product, which we try to compensate with the additional thickness-based conditions (see Sect. 2.5.2). In addition and similar to Alexandrov et al. (2010), the data is seasonally and regionally limited. Our measurements are confined to the western Arctic in early April and therefore, more measurements across the Arctic and the seasons are needed to evaluate 350 the spatial and temporal variability of sea-ice bulk density.

Impact on sea-ice thickness retrievals
Assuming all the other parameters for the conversion of freeboard to thickness remain the same, the average sea-ice bulk density estimates derived in this paper would result in 12.4 % and 16.7 % larger sea-ice thickness values for FYI and MYI, respectively, in comparison to A10. The effect is larger for thicker ice, for which snow depth plays a proportionally less important role.

355
Therefore, improving especially the MYI bulk density is important to derive accurate time series of sea-ice thickness and volume, as in the past, thicker MYI represented a larger fraction of the Arctic sea-ice cover. Kwok and Cunningham (2015) recognised the possibility of varying MYI density between the recent younger MYI and older MYI of the previous decades and discussed the impact of MYI density on sea-ice thickness and volume. Moreover, a potentially increasing degree of deformation may lead to a bias in the time series, as satellites underestimate sea-ice draft of deformed ice as shown by Belter et al. (2020) 360 and further discussed by Khvorostovsky et al. (2020). Deformed and unconsolidated sea-ice has an increased bulk density due to sea-water inclusions and thus, using a typical density of consolidated level sea-ice for deriving sea-ice thickness from satellite data will eventually lead to an underestimation. Given the thinner and younger sea-ice cover together with observed increase in sea-ice drift speed and deformation (Rampal et al., 2009;Spreen et al., 2011), there is a likely premise for systematic underestimation by current parametrisations of sea-ice density in regions where and at times when sea ice is deformed. A full 365 impact assessment of the sea-ice bulk density parametrisation on decadal sea-ice thickness data record is a logical next step but beyond the scope of this paper.

Outlook
To represent sea-ice bulk density range as a functional relationship to a parameter observable from space rather than fixed values based on sea-ice type classification, we parametrised sea-ice bulk density using sea-ice freeboard and obtained a significant 370 correlation. Opting for an exponential function instead of linear was beneficial for ensuring a better fit, capturing the high bulk density values at low freeboard values, and avoiding linear decrease to possibly unrepresentative and unphysical values at high freeboards. Parametrisation in Eq. (7) sets the limits of bulk density to 953.8 kg m −3 at zero sea-ice freeboard and approaching 881.8 kg m −3 at high freeboards. Figure 9 shows the parametrisation in Eq. (7) applied to the AWI Level-3 Collocated CryoSat-2 Sea Ice Product (Hendricks 375 and Ricker, 2020) from the winter 2018/2019 converting the monthly gridded sea-ice freeboard to sea-ice density. The resulting sea-ice density distribution had a smoother transition between ice types compared to the current ice-type dependent density classification of the retrieval algorithm. Difference between the density parametrised with Eq. (7) and A10 was positive overall during the winter except in the Central Arctic Ocean and locally in the Fram and Bering Straits in spring. The density difference was the largest on MYI in proximity to FYI but decreased toward spring.

380
Only less than 3 % of the airborne data set has a sea-ice freeboard value larger than 0.5 m with considerable spread in bulk density values and thus, introducing uncertainty to the parametrisation at high freeboard values. Constraining the parametrisation at high freeboards would require more data in deformed and multi-year sea-ice environments. However, that needs to coincide with a sufficient amount of open leads to ensure accurate conversion of surface elevations to freeboard from the ALS.
With the limitations of the current method, it is also not feasible to investigate cases of negative sea-ice freeboard due to the 385 possible presence of liquid water and altered dielectric properties affecting the retrieval of snow depth.
Our parametrisation improves upon the previous formulations of sea-ice density given the significantly larger number of data points, larger areal coverage, the variety of ice types including deformed sea ice, and the choice of predictor variable it is based on. Kovacs (1997)  space, or in the case of effective freeboard not even with a single in situ measurement, which makes them difficult to apply. Fig. 9, our freeboard-based parametrisation has potential for future applications in satellite altimetry. While we acknowledge that, due to the variability in snow mass and sea-ice thickness, our parametrisation may not be applicable 395 on sub-kilometre scales as reflected by the scatter in Fig. 8, we think that a sea-ice density parametrisation is a significant improvement upon a single value or fixed values based on ice type. In this paper, we decided to adopt a single-predictor parametrisation for the sake of simplicity. However, for future studies it could be worthwhile, e.g., to apply machine learning algorithms to the full parameter space to discover possible multi-variable relationships. Given the effect of sea-ice deformation on bulk density, including sea-ice surface roughness, a multi-variable approach could explain more of the variability.

Conclusions
The unique, collocated, multi-sensor measurements of the Arctic sea ice from the AWI IceBird campaigns allow us not only to observe sea-ice thickness, freeboard, and snow depth in high-resolution on regional scales, but also for the first time to estimate sea-ice bulk densities of different ice types from airborne measurements. Despite measuring the sea-ice bulk density indirectly by deriving it from other measurements, we are able to capture the effects of deformed ice on FYI bulk density. In the 405 current Arctic, the average late-winter FYI and MYI bulk densities are higher than and do not differ as much as earlier studies suggested partly due to including deformed ice in the analysis. Alexandrov et al. (2010) derived a difference of 34.7 kg m −3 whereas our measurements show only 26.1 kg m −3 providing yet one more indication and consequence that the Arctic seaice cover is getting younger. Satellite altimetry sea-ice thickness retrieval algorithms need to adapt to these changes in order to capture the sea-ice thickness and volume accurately, and to account for changes over the satellite radar altimetry record 410 spanning almost three decades. Taking advantage of the abundant measurements collected over different sea-ice types during two late-winter airborne campaigns, we are able to provide a parametrisation of sea-ice bulk density using sea-ice freeboard.
The single-variable exponential function presented here yields a smaller RMSE than the uncertainty of density values fixed by sea-ice type currently in use in large extent. With potential applications in sea-ice thickness retrieval from satellite radar altimetry, a density parametrisation alone does not completely solve the uncertainty problem in the freeboard-to-thickness 415 conversion. Together with improved knowledge of snow loading, they provide a path to decrease the uncertainty in observing sea-ice thickness and volume where the recent (CryoSat-2/ICESat-2 orbit resonance) and future (CRISTAL mission) advances in dual-altimetry will play a key role. In situ and airborne multi-sensor observations of various sea-ice parameters across the seasons will remain important to validate new approaches.
Data availability. Collocated sea-ice parameter data products, including measured total thickness, snow freeboard, snow depth, and surface temperature and derived parameters such as sea-ice thickness, sea-ice freeboard, and sea-ice bulk density, are available on PANGAEA Collocated CryoSat-2 Sea Ice Product is available here: ftp://ftp.awi.de/sea_ice/product/cryosat2/v2p3/.
Author contributions. AJ carried out the analysis with support from SH and RR and prepared the manuscript with input from SH, RR, LvA, TK, and CH. AJ, SH, and RR collected the data. SH applied the density parametrisation to the CryoSat-2 data and produced Fig. 9. LvA processed the EM-Bird data to identify level and deformed ice. TK conducted the backward-tracking of sea ice with ICETrack.