Estimation of Arctic Winter Snow Depth, Sea Ice Thickness and Bulk Density, and Ice Freeboard by Combining CryoSat-2, AVHRR, and AMSR Measurements

Information on snow depth on sea ice and bulk sea ice density is required to convert CryoSat-2 radar freeboard (<inline-formula> <tex-math notation="LaTeX">$F_{r}$ </tex-math></inline-formula>) into sea ice thickness (SIT). It is difficult to obtain their information on an Arctic basin scale; therefore, most CryoSat-2 SIT products largely rely on the distributions of snow depth and bulk sea ice density derived from parameterizations, which are based on sea ice type and climatological values. Several observational studies have found that the distributions of parameterized variables are inaccurate compared to the actual distributions. This study aims to develop a new type of retrieval algorithm for snow depth, SIT and bulk density, and ice freeboard in the Arctic winter by synergizing active CryoSat-2 with passive microwave and infrared measurements. Two parameterizations for the snow–ice thickness ratio and bulk sea ice density were combined with the hydrostatic balance and radar wave speed correction equations. Consequently, solutions for the four target variables were obtained and applied to different CryoSat-<inline-formula> <tex-math notation="LaTeX">$2~F_{r}$ </tex-math></inline-formula>, derived from empirical and waveform-fitting (WF) retracker algorithms. The retrieved thickness-related parameters based on <inline-formula> <tex-math notation="LaTeX">$F_{r}$ </tex-math></inline-formula> from the Lognormal WF retracker algorithm showed good agreement with the airborne snow depth, total freeboard, and mooring ice draft measurements. The retrieved multiyear sea ice bulk density was significantly higher than the value of 882 kg <inline-formula> <tex-math notation="LaTeX">$\cdot ~\text{m}^{-3}$ </tex-math></inline-formula>, which was used in the previous density parameterization, showing a higher agreement with values from in situ measurements. The spatial and interannual variabilities of SIT increased when the results from this study were compared with those based on previous parameterizations.

on F r from the Lognormal WF retracker algorithm showed good agreement with the airborne snow depth, total freeboard, and mooring ice draft measurements. The retrieved multiyear sea ice bulk density was significantly higher than the value of 882 kg · m −3 , which was used in the previous density parameterization, showing a higher agreement with values from in situ measurements. The spatial and interannual variabilities of SIT increased when the results from this study were compared with those based on previous parameterizations.

I. INTRODUCTION
T HE diminishing Arctic sea ice has been presented as clear evidence of rapid Arctic surface warming. Arctic sea ice extent has been continuously monitored using satellite passive microwave (PMW) measurements, showing significant decreasing trends, over the past four decades (1979-2019), regardless of the season; the trends for 1979-2019 were −2.7% and −12.9% per decade for March and September, respectively, relative to the corresponding climatological extent for the 1981-2010 period [1]. This shrinking sea ice is closely related to global climate change, extreme weather in the midlatitude region, and changes in the local turbulent flux between the ocean and atmosphere [2], [3], [4], [5].
Sea ice thickness (SIT) is as important as sea ice coverage because it is required to assess changes in the total amount of sea ice (i.e., sea ice volume) [6], [7]. The variability in sea ice volume affects freshwater content, which is essential for investigating the variability of the thermohaline circulation in the Arctic Ocean [8]. Moreover, when the largest area of the Arctic Ocean is covered by sea ice (especially in winter), the amount of energy transported from the warmer ocean to the colder atmosphere depends on the SIT (i.e., the insulation effects of sea ice). For example, a climate modeling study showed that over the period of 1982-2013, the Arctic amplification factor, which is the relative speed of the Arctic surface air temperature rise compared to the global mean temperature rise, increased by 37% under sea ice thinning conditions compared to that with fixed SIT [9]. As SIT is important in weather and climate studies, satellite missions have been conducted to construct Arctic basin-scale SIT data records [10], [11], [12], [13].
The Synthetic Aperture Radar Interferometer Radar Altimeter (SIRAL) onboard CryoSat-2, operated by the European Space Agency Paris, France, has provided continuous data on SIT since 2010 [12]. The SIRAL/CryoSat-2 provides ranging measurements to estimate ice freeboard, defined as the height from the sea surface to the snow-ice interface (or sea ice surface under snow-free conditions), which is converted into SIT through several additional steps in sophisticated algorithms [14], [15], [16]. The SIRAL/CryoSat-2 emits Ku-band radar waves to the snow-covered sea ice and sea surfaces. It then records the waveform, which is the intensity of the backscattered waves by the surface layer as a function of time. A retracker algorithm is used to analyze the measured waveform to obtain the range, which is defined as the distance between the satellite sensor and the target surface layer [14], [15], [16]. The radar freeboard, which is the ice surface elevation initially estimated from the retracker algorithm, is converted into ice freeboard, by correcting the reduced radar wave propagation speed in the snow layer on the sea ice. Finally, the ice freeboard is converted into SIT under the hydrostatic balance assumption with a priori information on snow depth on sea ice and bulk densities of materials, such as snow on sea ice, sea ice, and seawater, as inputs.
The sensitivity of CryoSat-2-estimated SIT to the input variables is investigated by several studies. Results showed that the snow depth and bulk density of sea ice are the primary contributors to uncertainty [13], [14], [17], [18], [66]. In order to obtain a priori information on snow depth and bulk sea ice density, most CryoSat-2 SIT are based on the following parameterizations.
A modified version of snow depth climatology constructed by Warren et al. [19] (hereafter referred to as W99) is the basis for most CryoSat-2 products. The W99 is based on snow depth measurements at Soviet drifting stations over the Arctic Ocean from 1954 to 1991. However, several studies found that the W99 snow depth was significantly thicker than the actual snow depth over the first-year sea ice (FYI) area in the 2010s [17], [18], [20]. Based on this finding, the original climatological value was halved over the FYI area, thus resulting in the modified W99 snow depth climatology (hereafter referred to as mW99). The mW99 has been commonly used to estimate SIT from the CryoSat-2 radar freeboard over the central Arctic basin [21], [22].
Recently, in addition to the FYI area, the snow depth during the freezing period on multiyear sea ice (MYI) has also reduced significantly compared to the W99 snow depth [23], [24], [25], [26]. This implies that the mW99 also causes an overestimation of the snow depth over the MYI area. Moreover, snow precipitation on Arctic sea ice has a large interannual variability that cannot be described using climatological data. In mW99, this snow depth variability exists only in regions where the sea ice type changes. Additionally, recent modeling and satellite observational studies have shown that snow depth trends differ according to region, with an increase and decrease over the western and eastern Arctic Ocean, respectively [26], [27]. Therefore, snow depth that varies spatiotemporally (also referred to as "dynamic snow depth") is preferable for the estimation of SIT [28], [29].
The widely used values for the bulk density of sea ice are 916.7 kg · m −3 for FYI and 882 kg · m −3 for MYI, as adopted from Alexandrov et al. [30] (hereafter referred to as A10). This density parameterization is based on the difference in physical characteristics between MYI and FYI. In general, the sea ice above the waterline of MYI has a lower density than that of FYI, due to more air pockets in MYI. Thus, the bulk density value of FYI is greater than that of MYI. However, a recent study by Jutila et al. [31] reported that the bulk density of sea ice varies horizontally and is not a constant value for certain ice types; additionally, it can be parameterized by sea ice freeboard. Moreover, it is unclear whether a bulk density of 882 kg · m −3 for the MYI is reasonable for the estimation of SIT [26], [32]. To calculate the representative MYI bulk density of 882 kg · m −3 , A10 assumed that the densities of the upper and lower layers (sea ice above and below the waterline, respectively) are 550 and 920 kg · m −3 , respectively, and the thicknesses of the upper and lower layers were assumed to be 0.3 and 2.6 m, respectively. However, 550 kg · m −3 for the upper layer implies that the porosity of the upper layer is approximately 40%, which would not be representative characteristics of MYI [31], [65].
These points show that there are problems with parameterizations for snow depth and bulk sea ice density that are commonly used to estimate SIT from the CryoSat-2 radar freeboard. Both parameters are dependent only on sea ice type; therefore, spatiotemporally varying snow depth and bulk sea ice density are not considered, indicating that these parameterized values are likely to have considerable uncertainty. Therefore, it would be beneficial to find an approach to estimate SIT from CryoSat-2 measurements, which includes more realistic considerations of the snow depth and bulk sea ice density.
It is preferable to use a dynamic snow depth that varies spatiotemporally. There are two types of dynamic snow depth products. The first is based on satellite PMW measurements, trained with Operation IceBridge (OIB) snow depth measurement data [24], which is available only for March and April. However, a rigorous validation was difficult to perform against OIB data, which provides snow depth measurements over a wide sea ice area and various sea ice types, because this PMW-based product is already dependent on OIB data. Model outputs can serve as an alternative to dynamic inputs [33], [34], [35]. However, their outputs are highly dependent on atmospheric reanalysis snowfall data, which possess large uncertainties and differences between models [36], [37].
Shi et al. [38] introduced a method that simultaneously estimates snow depth and SIT from combined CryoSat-2 and satellite radiometer measurements. This method uses a snow-ice thickness ratio (hereafter referred to as α) to convert ice freeboard into snow depth and SIT, without a priori information on snow depth. It has been shown that α is proportional (inversely proportional) to the temperature difference between the top and bottom of the snow (sea ice) layer during winter. Additionally, an empirical relationship that can estimate α from air-snow interface and snow-ice interface temperatures was obtained by analyzing monthly averaged buoy temperature profiles [38], [39]. The ability of α for simultaneous estimation was demonstrated by using measurements from an Advanced Microwave Scanning Radiometer (AMSR) and an Advanced Very High Resolution Radiometer (AVHRR) [26], [38]. Although the proposed method can provide plausible distributions of snow depth and ice thickness, it is still affected by the use of fixed bulk sea ice density values. Therefore, in addition to the use of α for SIT and snow depth estimation, a parameterization for bulk sea ice density that allows dynamic sea ice density should be considered for a more accurate estimation.
Therefore, the aim of this study is to develop a method that simultaneously estimates SIT, snow depth, ice freeboard, and bulk sea ice density, by synergizing active CryoSat-2 with passive AMSR and AVHRR measurements. In this study, it is shown that these four variables are achievable using synergized active and passive satellite measurements. The rest of this article is organized as follows. Section II describes the proposed algorithm and analysis methods. The descriptions of the data used are provided in Section III. The results and quality assessment of the SIT and bulk density, snow depth on sea ice, and ice freeboard are presented in Section IV. Sections V and VI present the discussion and conclusions of this study, respectively.

II. METHODS
This section describes the proposed approach for simultaneously estimating SIT, snow depth on sea ice, ice freeboard, and bulk sea ice density. The method assesses an analytical solution for the hydrostatic balance and radar wave speed correction equations with parameterizations for the snow-ice thickness ratio and bulk sea ice density. Each component is described in Sections II-A-II-G.

A. Hydrostatic Balance Equation
The hydrostatic balance equation describes the balance between the weight of the snow-ice column and buoyancy. It has been widely used for SIT estimation using satellite altimeter observations [13], [15], [22]. The hydrostatic balance equation is as follows: where H i , h s , and F i are the SIT, snow depth on the sea ice, and ice freeboard, respectively. The variables ρ w , ρ i , and ρ s are the bulk densities of the seawater, sea ice, and snow, respectively. In this study, the density of seawater (ρ w ) was assumed to be 1024 kg · m −3 . The details of ρ i and ρ s are provided in Section II-D.

B. Radar Wave Speed Correction for Snow Layer
The dependence of the radar wave speed on the medium should be considered when converting the radar freeboard F r into ice freeboard F i . The radar wave speed in the snow layer is slower than that in the atmosphere; therefore, the distance of the ice surface from the satellite sensor is estimated to be larger than the actual range, leading to an underestimation of the ice surface height. Thus, the corresponding correction should be applied to obtain the ice freeboard from the radar freeboard, using the following equation: where η s is the real part of the snow refractive index, which can be parameterized as a function of ρ s , according to Ulaby et al. [40] η s = (1 + 0.51ρ s ) 1.5 . (3)

C. Snow-Ice Thickness Ratio
In this study, the variable α, defined as the ratio of h s to H i (i.e., α ≡ h s /H i ), was used to constrain the hydrostatic balance equation, instead of the direct use of h s as input. To obtain α from satellite measurements, α is parameterized as a function of temperature at three different interfaces (such as air-snow, snow-ice, and ice-water interfaces) under the following two physical assumptions [38]: 1) the conductive heat flux is continuous at the snow-ice interface [41] and 2) the snow and sea ice layers both have a linear temperature profile on a monthly timescale during the winter [38]. Based on the above, the ratio α can be formulated as follows: where k s and k i are the thermal conductivities of the snow and ice layers, respectively; T snow is the temperature difference between the air-snow and snow-ice interfaces, and T ice is the temperature difference between the snow-ice and ice-water interfaces. The empirical relationship between α and the temperature difference ratio T snow / T ice was determined by analyzing 18 buoy thermistor string measurements with a 2-cm vertical resolution [39]. The suggested parameterization for α is As discussed, three interface temperatures are required to estimate α using (5). The two temperatures at the air-snow (T as ) and snow-ice interfaces (T si ) are available from satellite infrared and PMW measurements (see Section III). The ice bottom temperature was assumed to be −1.87 • C, which is a typical freezing temperature of Arctic seawater [39]. The monthly fields of α can then be estimated from the monthly averaged air-snow interface and snow-ice interface temperatures estimated using satellite data.

D. Bulk Densities of Sea Ice and Snow
In this study, instead of prescribing a constant bulk sea ice density, a parameterization for the sea ice density was introduced, allowing bulk density values to be changed for each target grid. It has been reported that the density of sea ice below the waterline is approximately uniform, regardless of ice type [42]. However, the density of sea ice above the waterline shows ice-type dependence [42]. The density of the upper layer of MYI is generally lower than that of FYI because air pockets are more prevalent in MYI than in FYI. Based on this, A10 suggested a parameterization that estimates bulk sea ice density by averaging the densities of the upper and lower layers of the sea ice, weighted by the thicknesses of the upper and lower layers, as shown in the following equation: where ρ u and ρ l are the densities of the upper and lower layers of the sea ice, respectively. Equation (6) implies that the bulk density of sea ice decreases as the ice freeboard fraction (i.e., F i /H i ) increases because ρ u is smaller than ρ l . In this study, the values of ρ u and ρ l were selected according to Timco and Frederking [42]. The lower layer density (ρ l ) was set to 920 kg · m −3 , regardless of the sea ice type. The ρ u values were set to 815 and 875 kg · m −3 for MYI and FYI, respectively. For snow density, a parameterization proposed by Mallett et al. [43] was introduced. Although this parameterization considers snow density to be spatially invariant over Arctic sea ice, it accounts for the seasonal evolution of snow density, as follows: where t is the number of months since October (e.g., t = 3 for January).

E. Combined System of Governing Equations and Parameterizations, and Its Solution
The hydrostatic balance equation contains four unknowns: ρ i , H i , h s , and F i . Three additional constraints can be considered: 1) the F i is a function of h s , as in (2) (i.e., wave speed correction). The refractive index of the snow layer, η s , is a function of ρ s ; 2) the estimated α provides information on the ratio of h s to H i ; and 3) the ρ i is a function of F i and H i , according to (6) (i.e., density parameterization). Consequently, there are four constraints (including the hydrostatic balance equation) and four unknowns; therefore, the equation system is now mathematically posed. The solution for H i can be written by combining (1)-(7), as follows: where the variables on the right-hand side are known (as previously mentioned). The derivation of (8) is provided in the Appendix. After the estimation of H i , the corresponding h s can be obtained by multiplying H i by α (i.e., h s = α H i ) The remaining variables are F i and ρ i , which can be obtained from (2) and (6), respectively, It is noted that H i , h s , and F i are dependent on both α and F r , whereas ρ i depends only on α, implying that the estimated ρ i is independent of the CryoSat-2 F r . Fig. 1 summarizes the data processing procedure described here. It shows that SIT, snow depth, ice freeboard, and bulk density of sea ice are estimated from CryoSat-2, AVHRR, and AMSR measurements.
The key point of the method proposed in this study is that the three unknown variables are not processed consecutively to obtain final estimates for SIT, but are simultaneously estimated by solving the combined equations. Therefore, the developed algorithm can provide spatiotemporally varying bulk sea ice density and snow depth. To validate the retrieval results against snow depth and total freeboard from OIB airborne measurements, and sea ice drafts from the Beaufort Gyre Exploration Project (BGEP) mooring measurements, the retrieved variables were converted into variables corresponding to the validation dataset. For example, the total freeboard (F t ) is calculated by adding F i and h s , and the sea ice draft (F d ) is obtained by subtracting F i from H i .

F. Retrieval Based on the Previous Parameterizations
As discussed in Section II-D, the proposed method is based on a different approach than the widely used algorithm for CryoSat-2 SIT. The SIT estimated in this study was compared with the SIT based on the mW99 snow depth and A10 sea ice density parameterizations. The latter was obtained from the identical CryoSat-2 radar freeboard (as used in this study), with mW99 snow depth, which was calculated by halving W99 over FYI, and the bulk sea ice density values from A10 for FYI and MYI, which were assigned according to the sea ice type classification (described in Section III-C). The remaining parameters were the same as in the proposed method.

G. Gaussian Error Propagation
The influences of the uncertainty in input parameters on the retrieval results were investigated using the Gaussian error propagation method, which can be written as follows: where Y is the retrieved variables such as h s , H i , F i , and ρ i , and σ is the uncertainty of the variable indicated as the subscripts (i.e., six input parameters: F r , T as , T si , ρ u , ρ l , and ρ s ). Each partial derivative in the above equation was numerically calculated and multiplied by typical uncertainties σ of inputs and assumed variables reported by the relevant references [38], [42], [44], [47], [64]. The spatially varying uncertainty of F r was obtained from the Threshold First Maximum Retracker Algorithm (TFMRA) with a 50% threshold value (TFMRA50) CryoSat-2 radar freeboard dataset (Section III-A). The uncertainty of T as was set to be 3.4 K, which is the standard deviation of the difference between satellite-derived T as and buoy-measured 2-m temperature [64]. The uncertainty of T si was set to be 1 K according to the validation result using the buoy-measured snow-ice interface temperature [47]. The uncertainty of ρ u was set to be 35 kg · m −3 for FYI and 95 kg · m −3 for MYI, which are half of the min-max range reported by Timco and Frederking [42]. Similarly, the uncertainty of ρ l was set to be 20 kg · m −3 . The uncertainty of ρ s was set to be 50 kg · m −3 [38]. Moreover, to investigate which input data or assumed parameters are the most responsible for the uncertainty in the retrieved parameters, the relative contribution of the six input parameters (C X ) was defined as follows: where the parameter X denotes the six input parameters. The results of the uncertainty analysis are provided in Section IV-D.

III. DATA A. CryoSat-2 Radar Freeboard
This study utilized the CryoSat-2 radar freeboard to simultaneously estimate SIT, snow depth, ice freeboard, and bulk ice density. Three different radar freeboards from different retracker algorithms (with their own biases) were considered to investigate which was the most compatible with the retrieval algorithm. The first is the TFMRA50 [44]. The TFMRA monthly radar freeboard data were obtained from the dataset titled "AWI CryoSat-2 Sea Ice Thickness (version 2.4)" available from the Alfred Wegener Institute, Bremerhaven, Germany, FTP site (ftp.awi.de/sea_ice/product/cryosat2/v2p4/). Another radar freeboard dataset used in this study is based on the waveform fitting (WF) retracker algorithm. This algorithm optimizes a simulated CryoSat-2 waveform to the observed waveform to obtain sea ice elevation and surface roughness [15]. As it assumes a Gaussian surface height distribution, this retracking method is hereafter referred to as a Gaussian WF. The dataset titled "CryoSat-2 Level-4 Sea Ice Elevation, Freeboard, and Thickness, Version 1 (RDEFT4)" was obtained from the National Snow and Ice Data Center (NSIDC), Boulder, CO, USA, website [21]. This dataset is provided on a 25km polar stereographic grid. Radar freeboard data were not included in this dataset. Therefore, the radar freeboard was reconstructed from the ice freeboard and snow depth data included in the dataset, according to the procedure described in the user's manual. The last dataset is also from the WF retracker algorithm, but it assumes a Lognormal surface height distribution (hereafter referred to as Lognormal WF) [16], [45]. The dataset is available on the British Antarctic Survey, Cambridge, UK, website [68]. The TFMRA50 and Lognormal WF radar freeboard data were regridded onto a 25-km polar stereographic grid. The January-February-March (JFM) datasets are available from 2011 to 2022, except for the Lognormal WF radar freeboard data, which is available until 2018. The differences in radar freeboard between the above three datasets used cannot be fully attributed to the retracking method itself; they are also related to differences in other parameters used in individual radar freeboard products, such as local sea level height.

B. AVHRR Air-Snow Interface Temperature
Arctic basin scale air-snow interface (i.e., snow surface) temperature data were obtained from satellite infrared radiometer measurements. Dybkjaer et al. [46] introduced a split-window algorithm to estimate the surface temperature over the Arctic Ocean, using AVHRR-measured brightness temperatures (TBs) at two window channels of 10.8 and 12.0 µm. Two daily surface temperature datasets were obtained from the Copernicus Marine Environment Monitoring Service (CMEMS; https://marine.copernicus.eu). For 2011-2019, the data were obtained from the "Arctic Ocean-Sea and Ice Surface Temperature REPROCESSED (SEAICE_ ARC_PHY_CLIMATE_L4_MY_011_016)" dataset. Data for 2020-2022 were from the "Arctic Ocean-Sea and Ice Surface Temperature (SEAICE_ARC_SEAICE_L4_ NRT_OBSERVATIONS_011_008)" dataset. These two datasets are available on the CMEMS website [69], [70]. The snow surface temperature data in these two datasets are based on an identical split-window algorithm, and the air-snow interface temperature data records for 2011-2022 are consistent. Daily air-snow interface temperature data in a 0.05 • grid format were regridded onto a 25-km polar stereographic grid.

C. AMSR-E/2 Snow-Ice Interface Temperature and Sea Ice Type
The snow-ice interface temperature data were obtained using AMSR-measured TBs at 6.925-GHz vertically and horizontally polarized channels, based on the method developed by Lee and Sohn [47]. In this study, ascending-and descendingaveraged L3 daily TB fields, measured from AMSR-E&2 for the JFM months of 2011-2022, were used. Due to an operational time gap between AMSR-E and AMSR2, data for JFM 2012 are not available. The AMSR TB fields in the 25-km polar stereographic grid format were obtained from the Japan Aerospace Exploration Agency (JAXA), Tokyo, Japan, through their FTP site (ftp.gportal.jaxa.jp).
To obtain consistent snow-ice interface temperature records for the study period, the following three preprocessing steps of TB fields were considered.
1) AMSR2 TBs were converted into AMSR-E-equivalent TBs, according to the intercalibration result reported by JAXA [48] in order to construct consistent AMSR-E and AMSR2 TB data records. 2) Atmospheric upwelling and surface-reflected downwelling radiance contributions within measured TBs were removed. Atmospheric TB contributions were calculated using the Satellite Data Simulator Unit (SDSU)-version 2.1 [49] with temperature and humidity profiles from the European Centre for Medium-Range Weather Forecasts (ECMWF), Reading, UK, ReAnalysis-5th Generation (ERA5) [50]. Detailed procedures for atmospheric correction can be found in [51]. 3) Data close to the coastline (within 100 km) were discarded to prevent land contamination effects due to the large footprint of the satellite PMW observations [52]. In this study, the sea ice type was identified using an AMSRderived emissivity difference between 10.65 and 18.7 GHz, based on a method introduced by Lee et al. [51]. The same preprocessing procedures described above were performed for L3 daily AMSR-E/2 TBs to calculate sea ice emissivity values. If the estimated vertically polarized emissivity at 10.65 GHz is greater than that at 18.7 GHz, the corresponding pixel is MYI; alternatively, it is FYI [51].

D. Sea Ice Concentration
To determine sea ice pixels, daily sea ice concentration (SIC) fields were obtained from the "NOAA/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration, Version 4" dataset [53]. The dataset consists of a combination of SICs from the National Aeronautics and Space Administration (NASA) Team [54] and the NASA Bootstrap algorithms [55] to synergize the strengths of each algorithm (detailed explanations of the dataset can be found in the user manual of [53]). The SIC data are available in a 25-km polar stereographic grid format. In this study, monthly fields of sea ice variables, including air-snow interface temperature, snow-ice interface temperature, and the emissivity difference between 10.65 and 18.7 GHz, were calculated by averaging daily fields under the condition that the SIC is greater than 98%.

E. Validation Datasets
To assess the accuracy of the retrieval results, the retrieved parameters were validated by utilizing measurements from the OIB mission and BGEP. The OIB mission is an aircraft mission that measured the total freeboard F t and snow depth h s over the Arctic Ocean, using a lidar altimeter (Atmospheric Topographic Mapper, ATM) and snow radar [56], [67]. All Locations of measurements used for the quality assessment of the retrieval results: (orange circles) OIB campaign, (green triangles) BGEP moorings, (red circles) Sever expedition over MYI, and (blue circles) Sever expedition over FYI.
available March OIB data from 2011 to 2018 were used in this study. Specifically, the OIB data for 2011-2013 were obtained from the "IceBridge L4 Sea Ice Freeboard, Snow Depth, and Thickness, Version 1" dataset [57], and data for 2014-2018 were obtained from the OIB Quick Look dataset [58]. The OIB measurements were discarded if a measured variable was missing, or if h s was greater than F t , because it is an incompatible condition for (1). Subsequently, daily gridded data were produced by collocating the OIB data on 25-km polar stereographic grids. This was done by assigning OIB measurements to the nearest 25-km grid and averaging the values of the collected measurements in the 25-km grid. The averaging process was performed only if a pixel had at least 500 OIB measurements, assuming that characteristics of sea ice and snow properties are spatially homogeneous (i.e., isotropic) in the 25-km satellite grid scale. In other words, it is assumed that 500 composites of 40-m-resolution OIB-measured snow and ice properties, which correspond to measurements along 20-km OIB track, represent general characteristics of snow and ice properties over the satellite grid scale. Fig. 2 shows the measurement locations of the OIB tracks (orange circles).
Additionally, sea ice draft data (defined as the thickness of sea ice below the waterline) from upward looking sonar (ULS) measurements from BGEP moorings were used to validate the obtained parameters. In this study, sea ice draft data from three BGEP mooring sites over the Beaufort Sea were used (see green triangles in Fig. 2). These mooring observations have provided year-round sea ice draft measurements from the ULS since 2003. The data from 2011 to 2021 were obtained from the BGEP website (https://www2.whoi.edu/site/beaufortgyre/). In order to obtain monthly data, the original ULS data, with a high temporal resolution of 2 min, were averaged for each month. The monthly ULS sea ice draft data were then compared with the nearest satellite retrievals.

F. Sever Expedition Dataset With Sea Ice Type Information
The plausibility of the obtained bulk sea ice density values was verified in Section IV-B. There are few available in situ sea ice density observations; therefore, the Sever expedition dataset analyzed in A10 was used. The Sever expedition provides in situ measurements of snow and sea ice parameters, such as snow depth, SIT, and ice freeboard for the period from the 1930 s to the 1980 s [59]. These three parameters can be converted to bulk sea ice density using (1), with the prescribed densities for snow and seawater. There is no information on sea ice type in the Sever dataset; therefore, A10 estimated the bulk density of FYI from the Sever dataset assuming that FYI is prevalent in the Eurasian-Russian Arctic Ocean, where the Sever expedition was performed. In the meantime, previous studies by Shalina and Sandven [60] and Jutila et al. [31] have shown that a 2-m SIT threshold can be applied to differentiate between MYI and FYI for Sever expedition measurements. For the years from 1980 to 1988, 644 measurements from Sever expedition provide necessary variables in this study, such as "prevailing snow depth," "prevailing sea ice thickness," "runway snow depth," "runway sea ice thickness," and "ice freeboard." By applying the 2-m threshold to prevailing SIT, which would represent the SIT of measurement sites, 160 and 484 cases were classified as MYI and FYI, respectively. For estimating the bulk sea ice density for each ice type, runway ice thickness, runway snow depth, and ice freeboard were analyzed. The runway data are the measurements collected near the runway of the aircraft, representing measurements for level ice; thus, data over the ice ridges and deformed ice are likely to be excluded. Fig. 2 shows the measurement locations of the Sever data and the corresponding sea ice types. The red and blue circles in Fig. 2 indicate MYI and FYI, respectively. A portion of 24.8% of the Sever data (160 out of 644 measurements) was from MYI, which is contradictory to the previous assumption that Sever provides FYI parameters made by A10.
Using the Sever measurements of SIT, snow depth, ice freeboard, and derived sea ice type, the bulk density values for MYI and FYI were calculated separately. This was performed by using (1), and assuming that the densities of seawater and snow were 1024 and 324 kg · m −3 , respectively, as suggested by Alexandrov et al. [30]. The density values obtained were compared with the satellite-estimated bulk sea ice density.

A. Retrieval Results of SIT, Snow Depth, Ice Freeboard, and Bulk Density of Sea Ice
The retrieval examples based on the Lognormal WF radar freeboard for January 2011 are shown in Fig. 3.
1) Sea Ice Thickness: The spatial distribution of retrieved SIT [ Fig. 3(a)] corresponds to a well-known characteristic of Arctic sea ice, which is that the sea ice is generally thicker for MYI than FYI. In January 2011, the mean values of the estimated SIT were 2.08 and 1.28 m for MYI and FYI, respectively. The standard deviations of SIT were 0.79 m for MYI and 0.42 m for FYI, indicating that the spatial variability of MYI thickness is larger than that of FYI. 2) Snow Depth: The spatial distribution of snow depth [ Fig. 3(b)] shows a high correlation with that of SIT [ Fig. 3(a)]. Over the MYI region, the mean snow depth was 0.19 m and the corresponding standard deviation was 0.10 m. In the FYI region, the mean snow depth was 0.09 m and the corresponding variability was 0.04 m, which are much smaller than those over MYI regions.
3) Ice Freeboard: Similar characteristics are also found in ice freeboard [ Fig. 3(c)]. Over the MYI region, the mean ice freeboard was 0.18 m with a standard deviation of 0.06 m. In the FYI region, the mean ice freeboard was 0.11 m with a standard deviation of 0.03 m, which are again much smaller than those over MYI regions.

4) Bulk Sea Ice Density:
The obtained pan-Arctic distribution of bulk sea ice density [ Fig. 3(d)] shows a spatially discrete pattern, whereas the three thickness-related variables discussed above show relatively continuous spatial distributions. It is because the ice-type-dependent upper layer density values were prescribed for the bulk sea ice density parameterization. The bulk density of MYI was lower than that of FYI, with means of 911.0 and 916.1 kg · m −3 , respectively. This corresponds to two recent findings that showed the following [31]: 1) the bulk density of MYI has a higher value than the value used in A10 (882 kg · m −3 ) and 2) the difference in bulk density between FYI and MYI is smaller than that based on A10. A higher variability of bulk density is found over MYI than over FYI. Specifically, the minimum and maximum values of bulk sea ice density are 909.5 and 912.6 kg · m −3 for MYI and 915.8 and 916.8 kg · m −3 for FYI, respectively.
A summary of these results and the corresponding statistics based on the TFMRA50 and Gaussian WF radar freeboards are provided in Table I

5) Dependence of the Retrievals to Radar Freeboard
Datasets: To investigate the dependence of the retrieval results on different types of CryoSat-2 radar freeboards, the SIT, snow depth, and ice freeboard estimated from the three radar freeboard datasets (discussed in Section III-A) were compared. The mean fields of each parameter during JFM of 2011-2018 were calculated, and the differences in the mean fields between the two different radar freeboards were obtained (Fig. 4). Compared to the retrieval results based on the Lognormal WF freeboard, the results based on the TFMRA50 and Gaussian WF radar freeboards exhibit different regional deviations. In general, the SIT, snow depth, and ice freeboard from TFMRA50 are thicker than those based on the Lognormal WF, over the regions adjacent to the north of the Canadian Archipelago and Greenland, where MYI is prevalent. This positive difference gradually increases toward land from the ocean. However, the retrieved parameters from the TFMRA50 radar freeboard were smaller over the rest of the area. The magnitude of the negative difference was smaller than that of the positive difference. In general, the retrieval results based on the Gaussian WF were smaller compared to those of the Lognormal WF, except for the relatively small areas near the coastal regions of Siberia and Alaska. These spatial differences in retrieval results between the three different radar freeboards are consistent with the study by Landy et al. [45]. This study investigated the difference between empirical and physical retracker algorithms, and the influence of ice surface roughness on the physical retracker algorithm.

B. Quality Assessment of the Estimated Variables
The comparison of snow depth and total freeboard between the retrievals and OIB measurements [ Fig. 5(a)  . The Gaussian WF radar results showed the largest negative bias compared to the OIB measurements (snow depth: −6.4 cm and total freeboard: −11.83 cm). It is interesting to note that the Lognormal WF radar results showed the smallest root-meansquared-difference (RMSD) against the OIB measurements (snow depth: 8.66 cm and total freeboard: 14.31 cm), although the results from the Lognormal WF radar freeboard have a larger negative bias than the results from the TFMRA50 radar freeboard. Although the overall statistics of TFMRA50 are superior to that of the other radar freeboards, it cannot be concluded that it is the best input for the developed algorithm because the OIB measurements can inherit their own bias, and cover only March.
To enhance the validity of the retrieval results, comparisons with ice draft measurements from BGEP moorings were also performed [ Fig. 5(g)-(i)]. The ice draft measurements are available for all months; therefore, the retrieved parameters can be compared with those for January and February, in addition to March. The comparison showed that the Lognormal WF showed the best statistics among the three radar freeboards (mean difference = −0.02 m, RMSD = 0.4 m, R = 0.47, and slope = 0.8). The mean difference in the TFMRA50 results was the poorest, which was the best in the OIB comparison. Considering the overall consistency with the OIB and BGEP measurements, it can be concluded that the Lognormal WF radar freeboard is the most compatible CryoSat-2 radar freeboard product for the developed algorithm (compared with the products based on TFMRA50 and Gaussian WF).
The bulk density values of sea ice estimated for 2011-2022 were compared with those derived from the Sever dataset (see Section III-F), and those from airborne measurements conducted by Jutila et al. [31] (hereafter referred to as J22). The monthly mean and standard deviation of the sea ice densities from these three data sources are shown in Fig. 6. The number of Sever data samples is small; therefore, outliers can considerably influence the calculation of the mean and standard deviation. Therefore, data beyond two times the standard deviation of monthly Sever data were discarded. The mean values of bulk sea ice density obtained in this study were 916.2 and 911 kg · m −3 for the FYI and MYI, respectively, showing small monthly variations. It should be noted that the retrieved MYI bulk density was smaller than the FYI bulk density by approximately 5 kg · m −3 . For March, Sever data are available for comparison with the retrieved densities from this study although there is a large time gap between the two datasets, with the Sever data (1984)(1985)(1986)(1987)(1988)(1989) and this study (2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021)(2022). The mean MYI bulk density for March obtained from the Sever dataset was ∼905 kg · m −3 , which is close to the value obtained in this study. The mean FYI bulk density for March obtained from Sever dataset was slightly larger than that obtained in this study. In general, the FYI bulk density based on Sever was larger than that obtained in this study. The FYI densities of A10 and this study were similar, whereas the MYI bulk density from A10 was significantly smaller than that from this study (∼30 kg · m −3 ). The MYI bulk density obtained in this study is far more consistent with the MYI bulk density from the Sever dataset than that from A10, which implies that the value of 882 kg · m −3 (widely used for MYI bulk density) is substantially smaller than the actual values. A similar conclusion was drawn when the Sever and J22 bulk sea ice density values (for April and May) were compared with those from A10.
The variations in bulk sea ice density based on the proposed method are small when compared with those based on Sever and J22. As shown in (6), the bulk sea ice density is a function of the ice freeboard fraction (i.e., F i /H i ) and the fixed sea ice density values for the upper and lower layers. The typical range of ice freeboard fraction is narrow, 7.91% ± 3.71% within the Sever dataset. Therefore, the influence of ρ u and ρ l on the bulk density calculation should be relatively large, so that the use of a fixed ρ u and ρ l is responsible for the small spatial variability in the bulk sea ice density. Consequently, a spatiotemporally varying parameterization for ρ u and ρ l should be further developed in order to reproduce the considerable variability observed in the Sever and J22 datasets. Additionally, the regional variations observed by satellites are likely to be smoothed, compared to in situ or airborne measurements. This is because there is a difference in the spatial scale between satellite observations and in situ or airborne measurements. Fig. 7(a) illustrates the spatial distributions of the mean SIT during JFM (2011-2018) based on the Lognormal WF radar freeboard. SIT ranges from ∼0.5 m over FYI to ∼4.1 m over MYI regions. In order to analyze the interannual variability of the retrieved parameters, the standard deviations of the retrieved parameters were calculated after removing the linear trend of the retrieved parameters during the study period. The spatial distribution of the detrended standard deviation of the estimated SIT over the study period is shown in Fig. 7(d).

C. Comparison With the Results Based on Previous Parameterizations
The interannual variability of SIT over MYI regions ranges from ∼0.4 to ∼1.2 m, greater than that over FYI regions which ranges from ∼0.1 to ∼0.3 m. The mean and detrended standard deviation fields of snow depth and bulk sea ice density are also shown in Fig. 7. The spatial distribution of snow depth variability is correlated strongly with SIT variability. The largest variabilities in bulk sea ice density are observed at the boundaries between the FYI and MYI where a sea ice type transition occurs. It is noted that the retrieved bulk sea ice density shows interannual variability over regions where the sea ice type did not change.
The spatial distribution of the three parameters, the same as in Fig. 7 except for the use of mW99 and A10 parameterizations, is shown in Fig. 8. The SIT in Fig. 8 exhibits a similar distribution to Fig. 7, but with a narrower SIT range [0.93-3.32 m, shown in Fig. 8(a)]. This difference is clearly depicted in the spatial distribution of the difference in mean SIT [ Fig. 9(a)]. Similarly, the magnitude of interannual variability of the obtained SIT from this study was smaller for the results based on the previous parameterizations [ Fig. 9(d)]. The detrended standard deviation of SIT with previous parameterizations ranges from ∼0.1 to ∼0.4 m [ Fig. 8(d)]. The mean snow depth based on mW99 showed a spatially different pattern compared with the retrieval results from this study [Figs. 7(b) and 8(b)]. Additionally, the mW99 distributions generally showed thicker snow than the snow depth distribution in this study. Over the Beaufort Sea, the difference between the two was up to ∼20 cm [ Fig. 9(b)]. However, over the Atlantic part of the central Arctic Ocean, the mW99 snow distributions were thinner than the snow distribution in this study [ Fig. 9(b)]. The distributions of mean bulk sea ice density from the previous and new parameterizations correlated well because they are both largely dependent on the sea ice type [Figs. 7(c) and 8(c)]. However, the bulk density values for MYI are considerably smaller in the A10 parameterization [ Fig. 9(c)]. The most important aspect is that there are differences in the interannual variability of snow depth and bulk sea ice density between the previous and the proposed parameterizations. No temporal variabilities were found in the mW99 and A10 parameterizations for the areas where a sea ice type change did not occur [ Fig. 8(e) and (f)], while variabilities exist over the corresponding areas in the retrieval results of this study.
In summary, the SIT estimated in this study showed larger spatial and interannual variabilities than those based on mW99 and A10 and the variability of estimated snow depth was generally larger than that of mW99 [ Fig. 9(d) and (e)]. This finding agrees well with the conclusions of other studies that the use of dynamic snow depth (instead of mW99) increases the spatiotemporal variability of the estimated SIT [28], [29], [61]. It is noted that the relatively larger interannual  variability of bulk sea ice density from the A10 parameterization over sea ice type transition regions has minimal contributions to difference in the interannual variability of SIT [ Fig. 9(d) and (f)].

D. Uncertainty Budget Analysis
The uncertainty in SIT (σ H i ) estimated from this study shows similar spatial patterns with the SIT itself (Fig. 10).
The magnitude of σ H i is generally greater over MYI than FYI. The highest σ H i is found over the north of the Canadian Archipelago, about 2.5 m. It has been found that the most responsible input parameter for σ H i is ρ l . It contributes approximately 65% of the total uncertainty over most Arctic regions. The uncertainty contribution of ρ l (i.e., Cρ l ) over the low-latitude regions is smaller than that over high-latitude regions because uncertainty contribution due to radar freeboard (C Fr ) increases over lower latitude regions due to increasing  random error in radar freeboard measurements. The greatest value of C Fr is found over the Beaufort Sea. The contributions of the T as , ρ u , and ρ s have a similar magnitude, and the T si contribution is negligible.
The spatial pattern of snow depth uncertainty (σ hs ) looks similar to σ H i ; however, major contributors to the total uncertainty are significantly different (Fig. 11). The T as is found to be the most significant uncertainty contributor to the retrieved snow depth. The reason why snow depth is highly dependent on T as is that the ratio α, which determines the snow depth portion, is sensitive to the T as . In a similar manner, T si has up to 15%-20% of contributions to σ hs . The second major contributor is again ρ l . These two parameters are responsible for the most uncertainty in retrieved snow depth. It is observed that the contribution of T as is greater over the FYI region and smaller over MYI regions, compared to that of ρ l .
In the case of ice freeboard uncertainty (σ Fi ), indeed, F r seems the most dominant contributor, especially over FYI  regions reaching up to 100% (Fig. 12). On the other hand, over the MYI region, the T as , ρ l , and ρ s equally contribute to σ Fi .  Bulk sea ice density shows higher uncertainty over the MYI region, about 22.5 kg · m −3 , than the FYI region, about 19.5 kg · m −3 (Fig. 13). The two major contributors are ρ l and ρ u . The contribution of ρ l is much greater than that of ρ u although ρ u has greater uncertainty compared to ρ l . It is because the thickness of the lower layer of sea ice is significantly greater than that of the upper layer of sea ice. However, it is noted that the contribution of ρ u relatively increases, especially for MYI regions. The ρ u contributes about 20% over the MYI regions. The uncertainty contribution of F r is zero because bulk sea ice density does not depend on the radar freeboard.
To generalize the analysis for January 2011, the uncertainty values of the estimated variables for each month and year, averaged according to sea ice type, are shown in Fig. 14 Similarly, regarding the relative contributions of input variables to the estimated uncertainty, the time series of the relative contributions are shown in Fig. 15 for FYI and MYI separately.
It was found that the analysis made for January 2011 is valid in general through the whole period, even though some temporal variabilities were found when looking into each month and year.

V. DISCUSSION
Although current bulk sea ice density retrieval results are largely dependent on sea ice type and the fixed density values of upper and lower layers, it is worthwhile to attempt estimating bulk sea ice density from satellite measurements to provide an Arctic basin-scale example for that widely used MYI bulk density value is significantly low, consistent with in situ and airborne observations which have limited coverage. While some studies and in situ measurements demonstrated that the actual bulk density of MYI is much greater than the widely used value of 882 kg · m −3 [31], [65], the authors hardly found studies or products related to CryoSat-2 SIT algorithm that reflect such recent findings, except for cases using identical density values of 915 kg · m −3 for both MYI and FYI assuming that MYI density is not as low as 882 kg · m −3 [21], [32]. Moreover, it is scientifically meaningful to reduce systematic bias originating from the low MYI density value, which can be canceled by overestimated snow depth on MYI. In short, by introducing the simple sea ice bulk density parameterization, this study showed that the mean value of MYI density estimated from satellite measurement is not as low as the widely assumed value in many CryoSat-2 SIT retrieval studies.
Further studies are required to parameterize the ρ u and ρ l to enhance the variability of the retrieved sea ice density. The use of scattering optical depth (SOD) at microwave frequencies to parameterize ρ u should be considered. This is because the SOD obtained from AMSR measurements reflects the amount of scattered radiant energy by air inclusion in the upper sea ice layer [62]. Since the density of sea ice decreases as the air amount increases, it is theoretically possible to relate ρ u to SOD. It is also recommended to introduce total freeboard measurements from ICESat-2 [63], as this may increase the degree of freedom in the proposed equation system by introducing additional total freeboard-related equation. Since the simultaneous estimation method using α is applicable to both CryoSat-2 and ICESat-2 [38], by tuning the fixed parameters pixel by pixel by constraining retrieval results from ICESat-2 and CryoSat-2 to be identical would provide the opportunity to obtain more realistic range of spatial variability of snow and sea ice density.
In addition, it was found that T si from satellite PMW measurements is not major uncertainty source although 6.9-GHz radiation may penetrate significantly into MYI. According to relative uncertainty contribution (C T si ), the uncertainty of the retrieved parameter due to potential error in T si about 1 K is not a major uncertainty contributor compared to other inputs except for snow depth on sea ice. In the case of snow depth, the relative uncertainty contribution of T si can reach about 15%-20%. It is because snow depth is related to the snowice thickness ratio, which depends on T si .
Lastly, the retrieval results are notably different in accordance with the radar freeboards used for retrieval, implying that the characteristics within the radar freeboard are as important as snow depth or bulk sea ice density. Therefore, it is important to determine a suitable combination of radar freeboard and other parameterizations (e.g., snow and sea ice properties) to neutralize their bias and error. For example, the widely used combinations found in two CryoSat-2 SIT products are: 1) the TFMRA50 radar freeboard, with mW99 snow depth and A10 sea ice parameterization [44] and 2) the Gaussian WF radar freeboard with mW99 snow depth, and a fixed sea ice density value of 915 kg · m −3 [21]. For both combinations, improving only one input parameter for the hydrostatic balance equation would not improve (or could even reduce) the accuracy of SIT estimations. Therefore, the quality of all input variables, including snow depth, bulk sea ice density, and radar freeboard, should be collectively enhanced to improve the estimation accuracy of SIT.

VI. CONCLUSION
This study developed a method to simultaneously estimate snow depth, SIT, ice freeboard, and bulk sea ice density using CryoSat-2, AVHRR, and AMSR measurements. Two parameterizations have been introduced in this study. These are the snow-ice thickness ratio α and the ice freeboard fraction-dependent bulk sea ice density, replacing the mW99 snow depth and the A10 sea ice density parameterizations, which are solely dependent on sea ice type. The ratio α was parameterized in terms of the interface temperatures based on buoy measurements. To obtain α from satellite observations, the monthly averaged fields of air-snow interface temperature from AVHRR, and snow-ice interface temperature from AMSR were used in this study. The bulk sea ice density was parameterized as a function of the ice freeboard fraction, and the updated representative density values for the upper and lower layers of the sea ice were used. By combining these two parameterizations with hydrostatic balance and radar wave speed correction equations, the analytical solutions for snow depth, SIT, ice freeboard, and bulk sea ice density were obtained.
These four parameters were obtained using the three types of CryoSat-2 radar freeboard data (i.e., TFMRA50, Gaussian WF, and Lognormal WF) and were highly dependent on the radar freeboards. The radar freeboard data from the Lognormal WF retracker algorithm are the most compatible with the proposed method, when the retrieval parameters are compared with the OIB-measured total freeboard, snow depth, and the BGEP-observed ice draft. The estimated bulk sea ice density values were in agreement with those derived from the Sever dataset. The mean values of the retrieved sea ice densities obtained in this study are more consistent with the Sever and J22 values than the values used in the A10 parameterization, particularly for MYI.
In this study, the SIT based on previous parameterizations was compared to that based on the new parameterizations. The estimated SIT based on the new parameterizations showed a larger spatiotemporal variability than that based on the previous parameterizations. Similar characteristics were observed for the estimated snow depth. The mW99 snow depth and A10 sea ice density are dependent on only the sea ice type; therefore, their temporal variabilities exist only over areas where an ice type transition occurs. However, the retrieved snow depth and bulk sea ice densities showed a relatively larger variability over wider regions compared to those variables based on mW99 and A10 parameterizations. Therefore, it can be concluded that the limited variabilities within the mW99 snow depth and A10 sea ice densities cause the smaller variability in the corresponding estimated SIT, under the condition of identical radar freeboard.