Impacts of snow data and processing methods on the interpretation of long-term changes in Baffin Bay sea ice thickness

In the Arctic, multi-year sea ice is being rapidly replaced by seasonal sea ice. Baffin Bay, situated between Greenland and Canada, is part of the seasonal ice zone. In this study, we present a long-term multi-mission assessment (2003-2020) of spring sea ice thickness in Baffin Bay from satellite altimetry and sea ice charts. Sea ice thickness within Baffin Bay is calculated from Envisat, ICESat, CryoSat-2 and ICESat-2 freeboard estimates, alongside a proxy from the ice chart stage of development that closely matches the altimetry data. We study the sensitivity of sea ice thickness results estimated from an array 5 of different snow depth and snow density products and methods for redistributing low resolution snow data onto along-track altimetry freeboards. The snow depth products that are applied include a reference estimated from the Warren climatology, a passive microwave snow depth product, and the dynamic snow scheme SnowModel-LG. We find that applying snow depth redistribution to represent small-scale snow variability has a considerable impact on ice thickness calculations from laser freeboards but was unnecessary for radar freeboards. Decisions on which snow loading product to use and whether to apply 10 snow redistribution can lead to different conclusions on trends and physical mechanisms. For instance, we find an uncertainty envelope around the March mean sea ice thickness of 13% for different snow depth/density products and redistribution methods. Consequently, trends in March sea ice thickness from 2003-2020 range from -23 cm/dec to 17 cm/dec, depending on which snow depth/density product and redistribution method is applied. Over a longer timescale, since 1996, the proxy ice chart thickness product demonstrates statistically significant thinning within Baffin Bay of 7 cm/dec. Our study provides further 15 evidence for long-term asymmetrical trends in Baffin Bay sea ice thickness (with -17.6 cm/dec thinning in the west and 10.8 cm/dec thickening in the east of the bay) since 2003. This asymmetrical thinning is consistent for all combinations of snow product and processing method, but it is unclear what may have driven these changes. 1 https://doi.org/10.5194/tc-2021-135 Preprint. Discussion started: 25 May 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
Arctic sea ice concentration and thickness has reduced significantly in recent decades (Kwok, 2018;Stroeve and Notz, 2018). 20 With a >50% decrease in multi-year ice (MYI) cover since the turn of the century, the Arctic is increasingly becoming dominated by seasonal ice (Kwok, 2018). The Arctic winter area coverage of seasonal ice has surpassed that of MYI, making the understanding of processes over first-year ice (FYI) to be as important as those over MYI (Jeffries et al., 2013).
Regions in the Arctic where the perennial ice (sea ice that remains at the end of summer) has declined rapidly include the Chukchi and East Siberian Seas, where there is nowadays virtually no perennial ice left (Stroeve and Notz, 2018). The Kara 25 Sea, Barents Sea and Baffin Bay are generally composed of seasonal ice. It is important to understand the interannual variation in the seasonal ice zone as these are generally regions of current and future shipping activity where thickness estimations are critical for safety (Christensen et al., 2018) and FYI regions have a strong influence on summer ice extent forecasting accuracy (Day et al., 2014;Serreze and Stroeve, 2015).
Most of the recent efforts to produce sea ice thickness estimates from satellite altimetry have focused on the central Arctic 30 sea ice pack which is dominated by MYI (Kwok and Cunningham, 2015;Laxon et al., 2013;Petty et al., 2020;Ricker et al., 2017;Tilling et al., 2018). This is largely due to increased density of in-situ snow measurements in the more central Arctic (Warren et al., 1999) and issues with waves/freeboard determination around the sea ice edge. There are thus generally fewer published records focused on providing and estimating ice thickness and trends in the peripheral seas (Mallett et al., 2020).
It is becoming particularly important to create a reliable sea ice thickness product for the seasonal ice zone with the ongoing 35 rapid replacement of MYI by seasonal ice.
Baffin Bay, situated between Greenland and the Canadian Arctic ( Fig. 1), is part of the seasonal ice zone. There is an import of multi-year ice from the Arctic Ocean through Nares Strait, but most of Baffin Bay is covered with first-year ice and the entire bay is ice free in summer. Baffin Bay plays a key role in modulating the freshwater flux from the Arctic basin to the Labrador Sea, which is a key location of North Atlantic Deep-Water formation, the deep-water component of the Atlantic Meridional 40 Overturning Circulation (Cuny et al., 2005;Curry et al., 2011;Fenty and Heimbach, 2013;Holland et al., 2001). Moreover, Baffin Bay is a busy shipping region (Christensen et al., 2018) and an important area for polar bear migration (Obbard et al., 2010).
Previous research has suggested an apparent asymmetry in Baffin Bay spring sea ice thickness, with a thicker ice pack in the west of the bay than in the east (Landy et al., 2017). Landy et al. (2017) also demonstrated possible long-term asymmetrical All current snow depth products typically have a much lower horizontal resolution (∼25-100 km) than the along-track freeboard measurements they are applied to (∼10-100 m for laser altimetry, ∼300-2000 m for radar altimetry). Previous studies have employed redistribution functions to represent small-scale variability not captured in the large-scale snow depth products such as new ice formation in leads and wind redistribution (Kurtz et al., 2009;Kwok and Cunningham, 2008;Petty et al., 2020). Kwok and Cunningham (2008) utilized a sigmoidal function of the large-scale snow depth and high-resolution 70 ICESat freeboard measurements. When the total freeboard is close to, or less than, the large-scale snow depth, the effective local snow depth is taken to be a fraction of the total freeboard as defined by a sigmoidal function. More recently, Petty et al. (2020) applied an updated version of Kurtz et al. (2009) on ICESat-2 measurements using an empirical approach of fitting a piecewise function that was determined based on airborne measurements from the NASA Operation IceBridge campaigns. Snow depth redistribution schemes have been applied to laser altimetry data, i.e. ICESat and ICESat-2, but have not yet been 75 tested for radar altimetry data.
In this study, we aim to reconcile the spring sea ice thickness derived from multiple satellite altimetry sensors and sea ice charts in Baffin Bay and produce a robust long-term record (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). We retrieve sea ice thickness from Envisat (2003-2012), ICESat (2003), CryoSat-2 (2011 and ICESat-2 (2019ICESat-2 ( -2020 to look at the regional long-term sea ice thickness change. We use this long-term record to analyse possible asymmetrical trends in sea ice thinning in Baffin Bay. 80 We investigate the impact of different snow depth and density products and redistribution methods on retrieved sea ice thickness from satellite altimetry along-track sea ice freeboard data. Moreover, we determine whether snow redistribution schemes need to be applied to radar freeboard data, as well as laser freeboard data, for accurate determination of the sea ice thickness. 2 Methods 2.1 Freeboard observations from satellite altimetry 85 We use freeboard data from the following satellite altimetry missions to determine March sea ice thickness: ICESat (2003( ), Envisat (2003( -2012, CryoSat-2 (2011CryoSat-2 ( -2020 and ICESat-2 (2019-2020) (Fig. 2). We limit our study to March only to focus on end-of-winter/early spring sea ice thickness when Baffin Bay is fully ice covered. Also, the ICESat mission only produced data for specific quasi-monthly campaign periods, including February-March. The freeboard observations from the various satellite sensors used in this study are described below.

ICESat
ICESat was launched in January 2003 by the National Aeronautics and Space Administration (NASA). The satellite orbits up to a latitudinal limit of 86°with a 91-day repeat period. The Geoscience and Laser Altimeter System (GLAS) on board ICESat provided surface elevation estimates within specific fall and spring campaign periods. The laser altimeter (wavelength 1064 nm had a laser footprint of ∼70 m in diameter spaced at ∼170 m intervals (Kwok et al., 2006). The laser stopped working in 2009.  Surface elevation relative to the TOPEX/Poseidon ellipsoid is provided by the National Snow and Ice Data Center (NSIDC) as the GLA13 product (Version 34, available at http://nsidc.org/data/icesat/ (Zwally et al., 2014)).
For the present study, elevation returns are filtered and corrected for geodetic and oceanographic biases, including geoid undulations, tides, dynamic topography of the ocean and the inverted barometer effects following Landy et al. (2017). We used elevation data where the sea ice concentration in the daily 25 km OSI-SAF global ice concentration reprocessing dataset 100 (Andersen et al., 2012) was above 20%. For days with more than 100 elevation samples within the Baffin Bay/Hudson Bay area, returns outside 4 times the standard deviation in relative elevation are considered outliers and are removed. Sea ice (floe) and sea surface (lead) elevation samples were separated based on differences in the reflective properties and relative elevation of these surface types using the approach in Landy et al. (2017), applying an adapted version of the algorithms from Kwok et al. (2007) and Kwok and Cunningham (2008). Leads are identified where the reflectivity is more than 0.2 units below the 105 background level in 25 km segments. A snow depth correction is applied to the elevation of the lead measurements (Kwok and Cunningham, 2008) and outlying leads are removed where the elevation is 3 times the standard deviation above the median or 2 times the standard deviation below the median. The sea level is interpolated between lead elevation measurements using cubic interpolation and smoothed using a 100 km window.
Total (snow plus sea ice) freeboard was calculated as the height difference between the ice surface elevation and the inter-110 polated sea level measurements. Negative total freeboard measurements were removed. Freeboard uncertainty depends on the distance to the closest lead elevation sample. Freeboard uncertainty for sea ice measurements within 25 km from the closest lead sample were determined from uncertainty in sea surface height obtained from a 25 km running standard deviation of lead measurements (Landy et al., 2017). For measurements located further than 100 km away from the closest lead, the freeboard uncertainty is constrained only by the height difference between the ICESat observation and the sea surface height. For mea-115 surements between 25-100 km from the closest lead sample the freeboard uncertainty increases linearly with distance between these values.

Envisat
Envisat was launched in March 2002 by the European Space Agency (ESA). The satellite orbits up to a latitudinal limit of 81.5°with an orbit period of 100 minutes with a repeat cycle of 35 days. The intertrack spacing at the Equator is 80 km. The 120 RA-2 altimeter on board Envisat includes a Ku-band pulse-limited radar altimeter (320 MHz), which theoretically penetrates snow if it is present on sea ice and measures sea ice elevation. The RA-2 altimeter has a nominal footprint of ∼2km (Connor et al., 2009). Communication with Envisat was lost in April 2012.
Envisat radar freeboard data on the satellite swath for March 2003-2012 were collected from ESA Climate Change Initiative (http://cci.esa.int/seaice) (Hendricks et al., 2018b). The surface-type classification is based on backscatter, leading-edge width, 125 and pulse peakiness with the thresholds being determined with a combination of unsupervised clustering and supervised classification of waveforms from the central Arctic (Hendricks et al., 2018b;Paul et al., 2018). Retracking of the sea ice waveforms was done by applying a threshold first-maximum retracker algorithm where the threshold was calibrated to CryoSat-2 radar freeboard results in the Central Arctic. The product includes a radar freeboard uncertainty estimate.

CryoSat-2 130
CryoSat-2 (CS2) was launched in April 2010 by ESA. The satellite orbits up to a latitudinal limit of 88°on a 369-day repeat period with a 30-day subcycle. The intertrack spacing at the Equator is 7.5 km (Wingham et al., 2006). The SIRAL altimeter (Ku-band) on board CryoSat-2 combines a pulse-limited radar altimeter, with synthetic aperture and interferometric signal processing (320 MHz). The footprint of CS2 is pulse-Doppler-limited ∼300 m along track and pulse-limited ∼1500 m across the track of the beam, with measurements spaced at ∼300 m intervals (Wingham et al., 2006). CryoSat-2 is still active.

135
CryoSat-2 sea ice freeboard data for March 2010-2020 were obtained from the Lognormal Altimeter Retracker Model (LARM) that accounts for variable roughness as described in (Landy et al., 2020). The sea ice freeboards obtained from this algorithm have been demonstrated to perform well over thin sea ice and are therefore appropriate for a predominantly seasonal ice regime such as Baffin Bay (Landy et al., 2020).
CryoSat-2 radar freeboard data are also available from the ESA CCI record and we have obtained radar freeboard data on 140 the satellite swath for March 2011-2017 from https://climate.esa.int/en/odp/#/dashboard (Hendricks et al., 2018a). This dataset shows some differences to the LARM processing, including significant variability between the two products along the coast.
However, it is included to examine the impact of different CryoSat-2 freeboard products on derived ice thickness and to provide a direct comparison with the earlier Envisat CCI data. The CryoSat-2 CCI product includes a freeboard uncertainty estimate.

ICESat-2 145
ICESat-2 (IS2) was launched in September 2018 by NASA. The satellite orbits up to a latitudinal limit of 88°on a 91-day repeat period with subcycles of 29 and 33 days. ICESat-2 carries the Advanced Topographic Laser Altimeter System (ATLAS) laser altimeter that operates in a split-beam configuration of three beam pairs which each include a strong and a weak beam (all at a wavelength of 532 nm). Each beam pair is separated by ∼3 km across-track, with a pair spacing of 90 m (Markus et al., 2017). The multiple beam pairs and high along-track sampling rate (10 kHz shots every 70 cm) provide improved spatial coverage over existing satellite altimetry missions over sea ice. The small footprint diameter of each of the beams (∼11 m) is useful for sea surface height measurements in the often narrow leads needed for sea ice thickness retrievals (Kwok et al., 2019).
The along-track freeboard height is computed by differencing sea ice heights from individual segments (produced through an aggregation of 150 photons along each of the beams) and local reference sea surface height determined from the lead height estimates within each beam (Kwok et al., 2020).

Snow data and redistribution
2.2.1 Snow depth and density 160 We apply three snow depth products on the satellite altimetry data: the Warren climatology (W99), a passive microwave (PMW) snow depth product and modelled snow depth forced by reanalysis snowfall. The Warren climatology (Warren et al., 1999, W99) ( Fig. 3a) is based on field measurements collected in the central Arctic Ocean. Webster et al. (2014) showed that the climatology is not representative of the Arctic-wide snow cover characteristics of recent years when compared against contemporary measurements from NASA's Operation IceBridge (OIB) airborne surveys. Analysis of the OIB snow depth 165 estimates suggested that snow depth on first year ice is approximately 50% the snow depth shown by W99 in those same regions (Kurtz and Farrell, 2011). As the climatology does not cover Baffin Bay, the regionally adapted W99 snow depth estimate for Baffin Bay is taken as the mean of snow depth on Pan-Arctic FYI in the W99 climatology and applied as a constant value over space and time. We multiply the mean FYI W99 snow depth climatology by 0.5 to account for the thinning of the FYI snow pack following Laxon et al. (2013). A constant snow depth value is unrealistic, but it acts as a reference for 170 other snow products and to compare with previous studies in which W99 has been used as the default snow depth dataset.
The PMW snow depth product was created from DMSP SSM/I-SSMIS brightness temperature measurements. DMSP brightness temperatures were chosen over AMSR brightness temperate data, because DMSP provides a continuous record over the study period and the PMW snow depth from DMSP better reconciles sea ice thickness from CryoSat-2 and ICESat-2 over the coinciding period of these two sensors (2019-2020). Moreover, it has been shown that DMSP provides a similar 175 snow depth distribution and there are no anomalous biases compared to AMSR-E in Baffin Bay (Landy et al., 2017). To create the PMW snow depth product over seasonal ice regions (Fig. 3b), the DMSP SSM/I-SSMIS Pathfinder Daily EASE-Grid (25km x 25km) Brightness Temperatures V2 at 37 GHz and 19 GHz vertical polarization were retrieved from NSIDC (https://nsidc.org/data/NSIDC-0032/versions/2). Daily PMW snow depth (2003-2020) was retrieved from the brightness temperatures according to Markus and Cavalieri (1998). This method relies on a strong correlation between decreasing brightness 180 temperature with increasing snow depth from increased volume scattering in the snow pack (Rango et al., 1979). Scattering decreases with increasing wavelength; therefore the ratio between brightness temperatures at 37 GHz and 19 GHz is used to determine snow depth (Markus and Cavalieri, 1998).
Modelled snow depths used in this study come from the recently developed Lagrangian snow-evolution model SnowModel- LG (Liston et al., 2020;Stroeve et al., 2020) (Fig. 3c). The model output used here were forced by MERRA2 atmospheric 185 reanalyses and NSIDC sea ice parcel concentration and trajectory datasets to produce daily snow distributions on a 25-km x 25-km grid (Liston et al., 2020). Processes including snow blowing, snow density evolution and ice dynamics are accounted for (Liston et al., 2020). We also utilized snow depth estimates from the NASA Eulerian Snow On Sea Ice Model (NESOSIM) (Petty et al., 2018), which uses similar large-scale reanalysis and satellite data input to SnowModel-LG, but more simple parameterizations of snow accumulation and loss. NESOSIM snow depth estimates were applied as a second check on the 190 effect of dynamic snow model results on sea ice thickness results in Baffin Bay.
A fifth snow depth product was created from CryoSat-2 and ICESat-2 freeboard observations, following the method of (Kwok et al., 2020) (Fig. 3d). Snow depth is expressed as the difference between the total freeboard (measured by ICESat -2) and sea ice freeboard (the radar freeboard measured by CryoSat-2 plus a correction for the propagation speed through the snowpack). This method can only be applied for the overlapping period between the two satellites (2019-2020) and is thus not 195 used for sea ice thickness processing but used as one of multiple references for evaluating the other snow and derived sea ice thickness products.
To calculate the sea ice thickness from the freeboard data, also the snow density is needed. Two snow density products were used. Firstly, a spatially constant snow density for each given day of the year was retrieved from the density observations from Warren et al. (1999). Secondly, the daily snow density outputs from SnowModel-LG were applied.

Snow redistribution
As the horizontal resolution of the snow depth products (∼25-100 km) is much lower than the resolution of the along-track freeboard measurements (∼15m -2km), snow redistribution provides a simple, albeit rudimentary, attempt to represent smallscale variability not captured in the large-scale snow depth products.
In this study we apply and test the effect of two snow redistribution methods together with a non-redistribution of snow.

205
The first redistribution function is the sigmoidal function by Kwok and Cunningham (2008). This function determines the effective snow depth as a fraction of the total freeboard. First, one determines the total freeboard-snow depth ratio. If this ratio is above 1.3, the effective snow depth is the same as the large-scale snow depth. Below 1.3, the effective snow depth-snow depth ratio is determined from the sigmoidal function based on the freeboard-snow depth ratio (Kwok and Cunningham, 2008). The functional fit was based on heuristic ideas regarding snow accumulation/loss (thinner ice is generally younger so has had less 210 time for snow to accumulate on) and is not based on empirical analysis (high-resolution snow/freeboard data was generally lacking when it was first developed). This function generally removes snow and does not attempt to conserve snow depth (i.e., mass) so can be considered as an aggregate snow loss-term on the large-scale snow depth.
The second redistribution we applied is the piecewise function that was determined based on airborne measurements of total freeboard from the Airborne Topographic Mapper (ATM) and snow depth from Snow Radar from NASA's Operation IceBridge airborne campaigns (Petty et al., 2020). This function was first adopted by Kurtz et al. (2009) to represent the full snow depth distribution on the regional scale, based on empirical evidence of a linear relationship between freeboard and snow depth up to a certain (threshold) freeboard. This redistribution function is applied iteratively to conserve snow depth within the given large-scale grid cell. This piecewise function was shown to improve the representation of sub grid-scale variability over no redistribution or other functional fits in the OIB data (Petty et al., 2020).

220
The sigmoidal function (Kwok and Cunningham, 2008) and the piecewise function (Petty et al., 2020) were designed to redistribute snow depth for high-resolution laser altimetry data. However, snow depth redistribution schemes have not yet been tested for radar altimetry data. We assessed whether there is a requirement and any benefit in using a radar freeboard snow redistribution scheme. We follow the methods outlined in Kurtz et al. (2009) to evaluate this, which is further described in Supplementary Materials S1.

225
The radar freeboard estimates did not show a relation with the snow depth on the length-scale of the pulse-limited footprint of the Envisat or the SAR-focused footprint of CryoSat-2 radar altimeters (Fig. S1). This shows that snow redistribution on radar altimetry freeboard data would not improve the conversion from ice freeboard to thickness, so we did not apply a redistribution function to the snow depth. For Envisat and CryoSat-2 data we simply used the snow depths at their native resolution (25 km) to convert freeboards to ice thicknesses.

Sea ice thickness
Along-track freeboard samples from all satellite altimetry products were discarded where the distance to land was less than 10 km or where the distance to the closest lead was more than 200 km to avoid land contamination and uncertain sea surface height from interpolation between leads.

Laser altimetry 235
When the snow is not redistributed on laser altimetry data, the large-scale snow depth can be larger than the small-scale total freeboard at individual along-track footprints, resulting in negative sea ice freeboards. Sea ice freeboard can be negative when the snowpack provides such a load that it depresses the ice surface below sea level. However, when the snow depth h s exceeds the ratio ρw (ρw−ρs) h f total -with ρ w density of sea water (1024 kg m −3 ), ρ s the density of snow, and h f total the total sea ice plus snow freeboard -this results in a physically impossible negative sea ice thickness. We assume that the snow depth measurement 240 for these samples is invalid and decrease the snow depth to be equal to the total freeboard measurement. This is true for 17-25% of the ICESat measurements and 8-9% of the ICESat-2 measurements when no redistribution was applied.
We calculated sea ice thickness (Fig. 4) from the along-track laser freeboard data with the W99 snow depth, the PMW snow depth and the SnowModel-LG snow depth, with the W99 and SnowModel-LG snow density products, and without redistribution of snow, with the sigmoidal function of effective snow depth (Kwok and Cunningham, 2008) and with the 245 piecewise snow depth redistribution (Petty et al., 2020).
Sea ice thickness h i was estimated from ICESat and ICESat-2 along-track freeboard observations using: where h f total is the total sea ice plus snow freeboard, h s is the depth of the snow layer, and ρ w , ρ i , and ρ s are the densities of sea water, sea ice, and snow respectively. Sea water density was taken as 1024 kg m −3 and sea ice density was obtained from 250 an ice thickness-dependent parameterization: ρ i (h i ) = 936 − 18h 0.5 i kg m −3 (Kovacs, 1996).

Radar altimetry
Sea ice thickness was estimated from Envisat and CryoSat-2 along-track radar freeboard observations with the W99 snow depth, the PMW snow depth product and the SnowModel-LG snow depth, with the W99 and SnowModel-LG snow density products, and without redistribution of snow using: where h fi is the ice-only freeboard, because the radar is assumed to penetrate the snow pack. To account for the lower propagation speed of the radar wave through the snowpack, a correction is applied to the radar freeboard: where c s is the speed of light through snow, parameterized by c s = c(1 + 0.51ρ s ) −1.5 (Ulaby et al., 1982). Sea water density 260 (ρ w ) was taken as 1024 kg m −3 and sea ice density was obtained from an ice thickness-dependent parameterization: Kovacs, 1996).

Uncertainty and gridding
The along-track sea ice thickness values were gridded on a 12x12 km (CryoSat-2 and ICESat-2) and 25x25 km (ICESat and Envisat, which provide a lower number of observations) EASE-grid using a mean filter inverse-linearly weighted by the sample 265 uncertainty and distance (Geiger et al., 2015) up to a distance of twice the resolution of the grid.
The sea ice thickness sample uncertainty was estimated by accounting for individual uncertainties in freeboard, snow depth and snow, sea ice and ocean water density (Landy et al., 2017). The individual uncertainty estimates for snow depth, snow density, sea ice density, and ocean water density from Landy et al. (2017) were applied. The random individual uncertainties in sea ice density and ocean water density and speckle noise (for the radar data) were assumed to be uncorrelated and could 270 be scaled by the number of observations within the grid cell. The individual uncertainties in the snow depth and snow density products were assumed to be correlated within the grid cell, as they were all interpolated from gridded products. The freeboard uncertainty for each of the sensors was assumed to be correlated within a single track, but the separate tracks are uncorrelated, so the uncertainty was scaled by the number of tracks. The individual uncertainties were combined, using Gaussian propagation of uncertainty, to generate a single uncertainty estimate for each sea ice thickness grid cell. The uncertainty in the mean 275 climatological sea ice thickness was found by taking the spatial average of the uncertainty estimates.

Canadian Ice Service Charts
The Canadian Ice Service (CIS) Charts provide a continuous estimate of sea ice concentration, stage of development (e.g. new ice, first-year ice (FYI), multi-year ice (MYI)) and forms of ice in the Canadian Arctic from 1968 to present. The ice charts are available weekly in summer (June-November) since the start of data collection) and monthly , bi-weekly (2006- ship and airborne observations. Over time the CIS integrated several data sources in the ice charts including satellite remote sensing data. The largest improvement in accuracy was achieved in 1996 with the inclusion of near-real-time satellite data from the Canadian remote sensing Earth observation satellite program RADARSAT. We make an estimate of sea ice thickness from the stage of development estimates from the CIS charts. The CIS charts were 285 retrieved from https://iceweb1.cis.ec.gc.ca/Archive/page1.xhtml?lang=en. We convert CIS .e00 vector files into 5 km resolution gridded 'raster' datasets. The charts provide the five most present stages of development and their partial concentration within each polygon. The stages of development are given with a range of thicknesses and the averages of these ranges are taken (Table 1). To estimate sea ice thickness from the ice charts, the mean of the thickness for each stage is taken weighted by the partial concentration of the stages of development. The uncertainty is estimated by the range of thicknesses given (Table 1).

290
The resulting product does not provide the actual sea ice thickness and the absolute values may not be reliable, but the product can be used to investigate relative spatial and temporal patterns in estimated sea ice thickness. Here we only use ice chart data following the launch of RADARSAT in 1996, so that results are consistent across the analysed time period.

Mean sea ice thickness and trends 295
As there are very little in-situ observations of snow depth or sea ice thickness available within Baffin Bay, it is not possible to determine unambiguously which of the snow depth products and redistribution methods give the most accurate estimation of sea ice thickness from satellite altimetry. We therefore look at the mean of the processing methods to illustrate the long-term inter-mission sea ice thickness in Baffin Bay (Fig. 5).
The March mean sea ice thickness of all the satellite records ( Fig. 5a-d) confirms a strong east-west asymmetry in sea 300 ice thickness across the bay (Landy et al., 2017). However, the magnitude of sea ice thickness is quite different between the There is no significant trend in the mean sea ice thickness (mean of all processing methods) in Baffin Bay from ICESat and 310 CryoSat-2 LARM (2003-2020).

Spread in results from processing methods
March mean ICESat sea ice thickness products from the eighteen different processing techniques show a similar spatial pattern of thick sea ice in the west of Baffin Bay and thinner sea ice in the east (Fig. 5a). However, the magnitude of ice thickness depends on the chosen snow depth product and redistribution method (Table 2). Generally, the PMW product results in thinner 315 sea ice than the constant W99 snow depth value applied for FYI (16% with no redistribution, 9% with the sigmoidal function, and 18% with the piecewise function). SnowModel-LG product results in slightly thinner sea ice thickness than the W99 product (9% with no redistribution, 5% with the sigmoidal function, and 10% with the piecewise function). Redistributing the snow depth results in a thinner sea ice pack with the piecewise function and a thicker sea ice pack with the sigmoidal function for all snow depth products. ICESat March mean sea ice thickness using the PMW snow depth product with the sigmoidal 320 and the piecewise redistribution gives respectively 12% thicker and 2% thinner sea ice than without the redistribution. The SnowModel-LG snow density product results in on average 4.6% thinner sea ice than the W99 snow density product.
The chosen snow depth product has a significant effect on the spatial pattern of the CryoSat-2 March sea ice thickness (see The SnowModel-LG snow density product results in on average 10% thinner sea ice than the W99 snow density product. The choice of CryoSat-2 radar freeboard product has a significant effect on the magnitude of sea ice thickness. Sea ice thickness is ∼30% thicker when generated with the CCI freeboard product than with the LARM freeboard product (see sup-  Moreover, the spatial pattern of sea ice thickness trends different for different combinations of snow depth and density products and snow redistribution methods (Fig. 7). The W99 snow depth and density with sigmoidal redistribution SIT (Fig.   7a) shows a strong asymmetrical thinning pattern, with strong thinning in the west of Baffin Bay and no thinning in the east.
The PMW snow depth with piecewise redistribution and W99 snow density (Fig. 7b) shows thickening in almost all of Baffin Bay, with some thinning in the North Water Polynya region. The SnowModel-LG snow depth and density with piecewise 345 redistribution (Fig. 7c) shows thinning in the northwest of the bay and thickening in the east and south.

Asymmetry and asymmetry trends
The sea ice thickness asymmetry across Baffin Bay is determined from the difference in mean sea ice thickness in western Baffin Bay and eastern Baffin Bay (Fig. 1). Although there is apparent east-west asymmetry in sea ice thickness across Baffin Bay in all products, the level of east-west sea ice thickness asymmetry across Baffin Bay depends on snow depth product   Table 4. Trend in west-east asymmetry for ICESat and CryoSat-2 LARM (2003-2020) for different snow depth products and snow redistribution methods. The first figure is with the W99 snow density value and the italic figure is with the SnowModel-LG snow density product.
The trends in asymmetry range from -41.6 cm/decade to -9.4 cm/decade (results in bold are significant within p<0.05 and redistribution technique (Fig. 8). There is also yearly variability in the strength of asymmetry. Trends in asymmetry, determined from the combined ICESat and CryoSat-2 LARM sea ice thickness data record using different snow depth products and redistribution methods, are summarized in Table 4. The mean of all the processing methods gives a significant decrease in west-east asymmetry over time of -28.3 cm/dec, caused by a thinning in the west of -17.6 cm/dec and a thickening in the east of 10.8 cm/dec. The trend in asymmetry depends on which snow data has been used. When the snow depth is not redistributed, 355 both the PMW and the SnowModel-LG snow depth products with W99 snow density result in a significant decrease in sea ice thickness asymmetry across the bay over the period 2003-2020. The SnowModel-LG product also results in a significant decrease in sea ice thickness asymmetry over time with redistribution of snow depth. However, the PMW product with the sigmoidal and with the piecewise snow depth redistribution function do not result in a significant reduction in asymmetry. This shows that any trends in asymmetry are dependent on the chosen snow depth and density product and redistribution method.

360
The asymmetry trends when combining ICESat and CryoSat-2 CCI sea ice thickness results are slightly stronger and range from -9 cm/dec (PMW, Piecewise) to -42 cm/dec (SnowModel-LG, no redistribution).

Canadian Ice Service Chart record
The estimated thickness and spatial patterns from the CIS ice charts are consistent with the ICESat and CryoSat-2 sea ice thickness results (Fig. 9a). Trends in FYI and MYI can be determined by the change in partial concentration of these ice types from the CIS charts ( Fig.   10). There has been a negative trend in MYI in the western part of the bay from ∼37% to ∼9% between 1996 and 2020, with this older ice being replaced by increasing FYI in most of Baffin Bay. Near the west coast, there has been some increase in

Mean state and trends
The March mean sea ice thickness of all the satellite records showed a strong east-west asymmetry in sea ice thickness across the bay. The characteristic east-west asymmetry in sea ice thickness that was found within Baffin Bay -both in the altimetry 380 thickness fields and the sea ice charts -can be explained by a combination of stage of development of the ice and freezeup date.  In spring the main stage of development within the bay is first-year ice, with a band of multi-year ice generally present along the west coast (Fig. 9b). As the bay is ice free in summer, this ice must be imported through Nares Strait and the Canadian Arctic Archipelago channels into Baffin Bay and drift with the Baffin Island Current along the west coast of the Bay (Bi et al., 2019;Kwok, 2005). Secondly, the freezeup date in the western part of the bay is earlier in the season than in the eastern part 385 of the bay (Fig. 9c), resulting in a longer period over which the ice can grow and thicken and more snow to accumulate. The sea ice generally drifts southwards along the western side of the bay (Tang et al., 2004), following the Baffin Island Current, and does not cross the bay to the east.
The CIS charts show a replacement of MYI with FYI within the bay in the first decades of the 21st century, with a decrease of MYI in the northwest of ∼37% to ∼9% between 1996 and 2020. This might be associated with a thinning of spring sea 390 ice within Baffin Bay, but no significant change in sea ice thickness was found from the altimetry products within the error margins for all inter-mission processing methods for the limited period multi-mission altimetry data are available. It is possible the MYI entering and transiting Baffin Bay has always included decayed and thinner ice floes that are not appreciably thicker than FYI floes measured with altimetry.
The freezeup onset displays a trend with earlier freezeup in the north-east and later freezeup in the south-east of Baffin Bay 395 over the past two decades. The freezeup date in the west of Baffin Bay does not display any obvious trend.
The stage of development from the CIS ice charts show strong agreement with the satellite altimetry sea ice patterns and offers the potential to be used as a proxy for the mean and regional distribution in sea ice thickness and spatial pattern (see CryoSat-2 and ICESat-2 (mean difference in sea ice thickness for overlapping years of 0.25 m, 24%) than Envisat and ICESat (mean difference in sea ice thickness for overlapping years of 0.77 m, 69%).
The CryoSat-2 sea ice thickness product shows a thick region in Melville Bay to the northeast of Baffin Bay in March for 7 of the 10 covered years (e.g. Fig. 8c), which is not as present in the sea ice thickness products from the other missions. A discussion of this feature is presented in the Supplementary Materials (S4).

Different processing methods
The long-term record of sea ice thickness from multiple satellite altimetry missions processed with different snow depth and density products and snow redistribution methods can lead to significantly different results, and therefore influence conclusions on total and regional sea ice thickness trends. The selected snow depth product influences both the found mean sea ice thickness, spatial patterns and trends. Processing decisions can introduce regional biases that, at least in thin ice areas, obscure decadal 425 trends and patterns. Care must be taken to estimate true product uncertainty envelopes. For instance, across the wide range of processing options tested here (for ICESat and CryoSat-2) we find a mean uncertainty envelope around the March mean sea ice thickness of 13% and a range of possible multi-mission trends of -23 cm/decade to +17 cm/decade.
When comparing the snow depth products to a determination of snow depth from the difference between CryoSat-2 radar freeboard and ICESat-2 laser freeboard in March 2019 and 2020 (Kwok et al., 2020) (Fig. 3d), the PMW snow depth product 430 shows the most similar pattern and magnitude. SnowModel-LG shows thicker snow depths along the coasts and the sea ice margin, and much thinner snow depths in the centre of the bay. SnowModel-LG also does not capture the west-east asymmetry that the PMW snow depth product shows and would be expected because of the longer period for snow to accumulate on the older ice in the west. Another snow model (NESOSIM, Petty et al. (2018)) was compared and shows a similar pattern in snow depth to SnowModel-LG but with thicker snow depth in the north of Baffin Bay (see S2). As there are no direct observations 435 of snow depth within Baffin Bay, a selection of the best snow depth product cannot be made. None of the snow depth products show a good reconciliation in sea ice thickness between the temporally overlapping CryoSat-2 and ICESat-2 products.
The trends in west-east sea ice thickness asymmetry strongly depend on the chosen processing methods and snow depth product, for this thin seasonal sea ice region. Most of the processing methods exhibit a significant trend in west-east sea ice thickness asymmetry across the bay, agreeing with the finding of Landy et al. (2017). The only method that does not show a 440 significant trend is the PMW snow depth with piecewise snow redistribution sea ice thickness.

Conclusions
This study produced a long-term multi-mission record of spring Baffin Bay sea ice thickness (2003-2020) with multiple snow depth and density products. The record shows asymmetrical sea ice thickness for all satellite altimetry products, with thicker sea ice in the west and thinner sea ice in the east of the bay. There is no significant trend in mean sea ice thickness. However,

445
this long-term record shows a significant decrease in the west-east asymmetry of -28.3 cm/dec. This paper showed that there are inter-mission biases in sea ice thickness in this region. The Envisat sea ice thickness results (1.86 ± 0.44 m) appear to overestimate the sea ice thickness compared to all other products (1.12 ± 0.47 m, 1.13 ± 0.47 m, and 1.30 ± 0.42 m) in Baffin Bay. This suggests that the CCI Envisat freeboard may not be as effective over the thinner FYI of Baffin Bay as it is in the Central Arctic and suggests further processing work on historic radar altimetry data is needed to create 450 reliable sea ice thickness products in the seasonal ice zone or in zones that have transitioned between multi-and first-year sea ice in recent decades.
Baffin Bay is part of the seasonal ice zone and it is becoming increasingly important to understand the trends and variability of sea ice in this region due to the rapid replacement of MYI by FYI that is both seen in Baffin Bay (a decrease of the concentration of MYI from 37% to 9% between 1996 and 2020) in this study from the CIS ice charts, and wider throughout 455 the rest of the Arctic (Kwok, 2018;Mallett et al., 2020).
We have compared the effects of applying different snow depth products and snow redistribution methods on sea ice thickness calculations from satellite altimetry. We show that different data processing techniques in satellite altimetry can lead to significantly different results in March mean sea ice thickness (ranging by ∼13%), the spatial pattern of sea ice thickness, and trends in sea ice thickness (ranging from -23 cm/dec to +17 cm/dec). Decisions on which snow depth product to use or whether 460 to use a redistribution function can influence the results and conclusions on physical mechanisms driving changes in the ice.
Having identified more consistent sea ice thickness distributions and magnitudes for the two years of CryoSat-2 and ICESat-2 overlap, it is clear that mission overlaps are vital for ensuring long-term SIT trends are robust. None of the used snow depth products provide a good reconciliation between CryoSat-2 and ICESat-2 in mean sea ice thickness nor spatial variability. This shows that more observations on snow depth and density in the seasonal ice zone are necessary to create and validate a suitable 465 snow depth product for this region.
Author contributions. IAG and JCL conceptualised the study. IAG carried out the main analysis and wrote the paper. JCL, AAP and NTK contributed to the interpretation of the results. AAP and NTK helped develop the methodology. All authors contributed to revising the manuscript. JCS provided SnowModel-LG data and AAP provided NESOSIM data.