Winter Arctic sea ice thickness from ICESat-2: upgrades to freeboard and snow loading estimates and an assessment of the first three winters of data collection

Reliable basin-scale estimates of sea ice thickness are urgently needed to improve our understanding of recent changes and future projections of polar climate. Data collected by NASA’s ICESat-2 mission have provided new, high-resolution, estimates of sea ice freeboard across both hemispheres since data collection started in October 2018. These data have been 15 used in recent work to produce estimates of winter Arctic sea ice thickness using snow loading estimates from the NASA Eulerian Snow On Sea Ice Model (NESOSIM). Here we provide an impact assessment of upgrades to both the ICESat-2 freeboard data (ATL10) and NESOSIM snow loading on estimates of winter Arctic sea ice thickness. Misclassified leads were removed from the freeboard algorithm in the third release (rel003) of ICESat-2 freeboard data, which increased freeboards in January and April 2019, and increased the fraction of low freeboards in November 2018, compared to rel002. 20 These changes improved comparisons of sea ice thickness (lower mean biases and standard deviations, higher correlations) with monthly gridded thickness estimates produced from ESA’s CryoSat-2 (using the same input snow and ice density assumptions). Later releases (rel004 and rel005) of ICESat-2 ATL10 freeboards result in less significant changes in the freeboard distributions and thus thickness. The latest version of NESOSIM (version 1.1), forced by CloudSat-scaled ERA5 snowfall, has been re-calibrated using snow depth estimates obtained by NASA’s Operation IceBridge airborne mission. The upgrade from NESOSIM v1.0 to v1.1 results in only small changes in snow depth which have a less significant impact on thickness compared to the rel002 to rel003 freeboard changes. Finally, we present our updated monthly gridded winter Arctic sea ice thickness dataset and highlight key changes over the past three winter seasons of data collection (November 2018 - April 2021). Strong differences in total winter Arctic thickness across the three winters are observed, linked to clear differences in the multiyear ice thickness at the start of each winter. Interannual changes in snow depth provide significant 30 impacts on our thickness results on regional and seasonal scales. Our analysis of recent winter Arctic sea ice thickness variability is provided online in a novel Jupyter Book format to increase transparency and user engagement with our derived gridded monthly thickness dataset.


Introduction 35
NASA's Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) is providing significant advances in our ability to monitor Earth's fast-changing sea ice cover. The Advanced Topographic Laser Altimeter System (ATLAS) onboard ICESat-2 measures surface elevation at high resolution (individual laser footprints of ~11 m, Magruder et al., 2020) and high precision (< 2 cm over sea ice flat surfaces, Kwok et al., 2019a), with dense along-track sampling (70 cm along-track from the 10 kHz pulse repetition rate, Neumann et al., 2019). ATLAS was designed in part to obtain accurate and routine estimates of sea ice 40 freeboard, the vertical extension of sea ice above local sea level, across the polar oceans (Markus et al., 2017). Sea ice freeboard can typically range from millimeters to tens of centimeters depending on the region or season profiled. ICESat-2 benefits from extensive polar coverage (profiling up to 88 degrees N/S, monthly sub-cycle) and has collected year-round data with minimal downtime since production started in October 2018. ICESat-2 sea ice height and freeboard data are provided in the official ATL07 (Kwok et al., 2021a) and ATL10 (Kwok et al., 2021b) products respectively. The first winter 45 season of ICESat-2 Arctic Ocean sea ice freeboards (ATL10) was presented in Kwok et al., (2019b), highlighting the regional and seasonal freeboard distributions obtained by ICESat-2.
Validation of the ATL07 and ATL10 products is on-going. ATL07 sea ice heights showed very strong agreement (0 cm mean differences, correlation coefficients of 0.97 to 0.98) with coincident airborne data collected by NASA's Operation IceBridge north of Greenland and the Canadian Archipelago in spring 2019 (Kwok et al., 2019a). The freeboard agreement 50 was more modest (mean differences of 0 to 4 cm), although the comparisons were hindered by the lack of available leads to reliably determine a local sea surface in either product. Additional analysis of the ATL07/10 surface classification scheme using imagery collected by the Copernicus Sentinel-2 mission, provided evidence of high skill in lead classification (Petty et al., 2021), a key part of the freeboard determination procedure. However, both the spring 2019 OIB Arctic campaign comparisons (Kwok et al., 2021d) and Sentinel-2 imagery assessments (Petty et al., 2021) highlighted errors in the 'dark 55 lead' classification in ATL07/10. Briefly, it was hypothesized that low/optically thin clouds in these regions attenuate the photon rate around these segments due to increased atmospheric scattering, tricking the empirical threshold-based classification algorithm into characterizing height segments over sea ice as dark leads. High photon rate specular leads are now the only lead types used to derive sea surface and thus freeboard in the Release 003 and subsequent sea ice products (Release 005 at the time of writing) while a possible filter for the dark-lead segments is developed and tested. The impact of 60 this change was an increase in freeboard in ATL10 of 0 to 3 cm depending on the season/region analyzed, as well as a decrease in coverage due to the reduction in sea surface tie-points (Kwok et al., 2021d).
Measurements of sea ice freeboard are typically collected as a means of estimating sea ice thickness (see schematic in Figure 1). This is conventionally achieved by combining freeboard measurements with ancillary estimates of snow loading (snow depth and density), sea ice density and an assumption of hydrostatic equilibrium (e.g. Giles et al., 2007;Kwok 65 and Cunningham, 2008;Laxon et al., 2013;Kwok, 2018). Sea ice thickness was estimated from Release 002 ATL10 freeboards using external snow loading estimates from the NASA Eulerian Snow on Sea Ice Model (NESOSIM) v1.0 and https://doi.org/10.5194/tc-2022-39 Preprint. Discussion started: 16 February 2022 c Author(s) 2022. CC BY 4.0 License. modified versions of the Warren climatology (Petty et al., 2020. The Feb/March 2019 ICESat-2 thicknesses were ~10 cm thinner than Feb/Mar 2008 ICESat thickness estimates, alluding to a possible decline in end-of-winter Arctic sea ice thickness over this 11-year period. However, the P2020 thickness estimates were also significantly thinner than those 70 produced using radar freeboard measurements from ESA's CryoSat-2 using the same input assumptions (tens of cm biases depending on the month and product analysed). Significant biases still exist in satellite-derived estimates of sea ice thickness, even those based on the same satellite sensor (e.g. radar altimetry data from ESA's CryoSat-2 mission, Sallila et al., 2019;Petty et al., 2020) which have limited their utility to-date, e.g. for constraining or calibrating polar climate projections (e.g. SIMIP Community, 2020). 75 The thickness results presented in P2020 used NESOSIM v1.0 snow loading forced by ECMWF ERA-Interim (ERA-I) snowfall (Dee et al., 2011). However, ERA-I production ended in August 2019 and was superseded by ERA5 (Hersbach et al., 2020). While ERA5 total precipitation is similar to ERA-I over the Arctic Ocean (Wang et al., 2019;Barrett et al., 2020), ERA5 produces relatively more snowfall (and thus less rainfall) compared to ERA-I, especially in the Atlantic sector of the Arctic (Wang et al., 2019). Additional developments and calibration of NESOSIM have been carried out to 80 upgrade NESOSIM (v1.0 to v1.1) and extend the derived ice thickness product beyond the first winter season (2018/2019) presented in P2020 which we present here. The significant changes in ATL10 freeboards and the availability of updated NESOSIM snow loading warrants an updated winter Arctic sea ice thickness assessment. ATL10 and NESOSIM v1.1 output are now also available from fall 2018 through to spring 2021, providing three winter seasons of data to assess. The main objectives of this paper are to: (i) 85 highlight upgrades to the ICESat-2 ATL10 freeboard product and NESOSIM v1.1 snow loading and assess their impact on winter Arctic sea ice thickness; (ii) carry out updated comparisons against CryoSat-2 derived thickness estimates; and (iii) assess monthly gridded thickness data from the past three winter seasons across the entire Arctic Ocean. The monthly gridded thickness analysis is also available online in a Jupyter Book format (https://nicolekeeney.com/icesat2-book) to increase transparency and user engagement in our analysis of these data. 90 2 Input data and upgrades

ICESat-2 ATL10 freeboards
We use the ICESat-2 ATL10 sea ice freeboard product (currently at Release 005, rel005), which is disseminated through the National Snow and Ice Data Center (NSIDC) (Kwok et al., 2021a). ATL10 is the end result of a series of algorithms that convert the primary geolocated photon product (ATL03, Neumann et al., 2019), to sea ice height and type (ATL07, Kwok et 95 al., 2021a), and then sea surface height and freeboard (ATL10, Kwok et al., 2021b). Briefly, the ATL07 algorithm subtracts a mean sea surface and time-varying ocean tide and inverted barometer corrections from ATL03, then aggregates and windows 150 photons around this corrected surface along each beam independently. ATL07 then extracts a best-guess Gaussian height distribution (convolved with the expected system response) to the photon height histogram to determine a single 'segment' height and various metrics summarizing the goodness of fit and radiometry (e.g., photon rate) of each 100 segment. This photon aggregation results in data with variable segment lengths of, on average, ~15 m for the strong beams and ~60 m for the weak beams (Kwok et al., 2019b). The spatial resolution of the individual segments can be estimated by adding the individual laser footprint size of ~11 m (Magruder et al., 2020) to the segment length, i.e., a mean of ~25 m for the strong beams and 70 m for the weak beams. An empirically based decision-tree algorithm is used to discriminate the height segments as either sea ice or sea surface/lead (Kwok et al., 2016). More details of the surface classification scheme 105 are available in Kwok et al., (2021d) and Petty et al., (2021), while the complete processing methodology is available in the Algorithm Theoretical Basis Document (ATBD) for sea ice products (Kwok et al., 2021c).
ATL10 converts adjacent sea surface segments into lead groups (to reduce noise in the lead height estimate) and then averages these into 10 km along-track reference sea surface height estimates along each beam. Sea ice freeboard is calculated as the difference between the individual ice height segments and the local sea surface height, independently for 110 each beam. The ICESat-2 beams are arranged in 'strong' and 'weak' beam pairs with each beam pair separated by ~3.3 kilometers in the across-track direction and the strong/weak beams separated by ~90 m across-track and ~2.5 km alongtrack. The weak beams are around 4 times lower energy (lower photon rate) than the strong beams. In this study we utilize only the strong beams to ensure the highest possible data quality. 115

ATL10 upgrades
The ICESat-2 sea ice products are continuously being updated as new assessments on the data are undertaken. All ICESat-2 products currently follow the same nominal release schedule (~6-12 months), so release updates are not necessarily based on the significance of the changes or improvements made to the given product. All new release data are processed and released from the start of the mission (October 14 th , 2018) onwards, until the production of a new release begins. The sea ice 120 thickness results presented in P2020 utilized rel002 ATL10 data, and differences with thickness estimates produced using rel001 ATL10 were noted to be negligible. As discussed earlier, in rel003 ATL10 (and subsequent releases), dark leads have been removed as possible sea surface height segments, since false positive classifications were found in the presence of clouds, resulting in 0-3 cm freeboard changes at basin-scales and some loss of coverage, especially within the more consolidated central Arctic ice pack (Kwok et al., 2021d). This is arguably the biggest change in the ICESat-2 sea ice 125 products to-date. The rel003 ATL10 data also included a relaxing of the height/freeboard quality flag (from 3 to 4), which means height segments with a poorer fit (generally segments from ridges with a more variable and complex height profile) are now included to increase retrieval counts over ridged ice regimes.
In rel004, most of the updates involved changes related to the treatment of the solid earth tides (a transition of ATL07 into a tide-free system to be consistent with ATL03). This caused a significant change in the magnitude of the 130 heights reported in ATL07 and ATL10, but as freeboard is a relative measurement, this was not expected to impact the reported freeboards. In rel005 ATL10, the only changes relevant to freeboard determination include improved calculation of the 10 km reference surface location to the centre of each section (effectively a bug fix). The rel005 data now also includes https://doi.org/10.5194/tc-2022-39 Preprint. Discussion started: 16 February 2022 c Author(s) 2022. CC BY 4.0 License. data from previously held granules where known satellite calibration scans were occurring somewhere along the granule.
New automated pointing angle and calibration scan filters were introduced in rel005 to ensure only data within each granule 135 experiencing degraded performance are filtered out, instead of withholding entire data granules. Most other developments in rel003 to rel005 ATL10 can be categorized as minor bug fixes and are listed in the ATBD change log (available since rel004). It should also be noted that every new release of ATL07/10 has coincided with a new ATL03 release, meaning ATL07/10 release upgrades also reflect upgrades to the underlying ATL03 photon heights (e.g., improvements in geolocation). 140 In the Supplementary Information ( Figure S1) we provide an assessment of the coverage change from rel002 to rel005 by counting the number of 10 km reference sea surface points available across the 4 releases from all data collected by the strong beams between November 2018 and April 2019. The analysis provides further evidence of the decline in coverage between rel002 and rel003. The rel003 to rel005 coverage differences are sporadic and linked mainly to the inclusion of calibration scan data granules. Calibration scans occur mainly over lower latitudes but can occasionally extend 145 over the Arctic ice pack -data during these scans are generally considered degraded i.e., heights with sub-nominal geolocation quality. Automated calibration scan filtering was introduced in rel005 to exclude these data more reliably and ensure only the highest quality height returns are utilized. In Figure S2 we provide a beam coverage assessment over the same time period using rel005 data only, highlighting the consistently higher coverage provided by the strong beams compared to the weak beams (in this first winter of data collection), and especially the middle strong beam. 150

NESOSIM
We use snow depth and density estimates from the NASA Eulerian Snow On Sea Ice Model (NESOSIM) (Petty et al., 2018a which is publicly available on GitHub (https://github.com/akpetty/NESOSIM). NESOSIM was developed primarily in preparation for the launch of ICESat-2, to enable timely production of snow depth and density estimates for seaice thickness retrievals using a simple snow accumulation model framework. NESOSIM includes two vertical snow layers 155 and several simple parameterizations (accumulation, wind packing, advection-divergence, blowing snow loss) to represent the expected primary sources and sinks of snow on Arctic sea ice during the accumulation season. Summer melt processes are currently neglected, so the model is typically run between September and the end of April. NESOSIM v1.0 was first presented in P2018 and the output using this v1.0 framework was used in P2020 to produce snow loading needed to convert ATL10 freeboards (rel002) to sea ice thickness from October 2018 to April 2019. The NESOSIM v1.0 output used in P2020 160 was forced with snowfall, winds and near-surface air temperature (to scale the initial snow conditions) from ERA-I (Dee et al., 2011), sea ice concentrations from the NASA Climate Data Record (CDR) version 3 (Meier et al., 2017), and ice drifts from the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) Ocean and Sea Ice Satellite Application Facility (OSI SAF) (Lavergne et al., 2010) which were all regridded to a 100 km x 100 km Arctic Ocean domain. 165

NESOSIM upgrades
Here we describe recent upgrades made to NESOSIM which has been tagged as a new version 1.1 (v1.1) code release (https://github.com/akpetty/NESOSIM/releases/tag/v1.1, archived at https://doi.org/10.5281/zenodo.4448356). Key updates in NESOSIM v1.1 include: CloudSat snowfall scaling (described more below, Cabaj et al., 2020), a new wind-atmosphere blowing snow loss term, an extended Arctic domain to cover the full extent of the Arctic peripheral seas, an improved 170 smoothing filter to reduce noise in the dynamic snow budget terms, an upgrade to Python 3, and various minor bug fixes.
Much of the NESOSIM v1.1 development was motivated by the need to recalibrate NESOSIM using ERA5 forcings (Hersbach et al., 2020), now that ERA5 has succeeded ERA-I (ERA-I data production ended in August 2019) and given reports of increased ERA5 snowfall compared to ERA-Interim (Wang et al., 2019;Cabaj et al., 2020). ERA5 is thought to offer improvements over ERA-I related to improved cloud representation, an updated assimilation scheme and higher spatial 175 resolution (Hersbach et al., 2020). Regardless, whether ERA5 exhibits a high snowfall bias over the Arctic or ERA-I a low bias is still uncertain and likely regionally dependant. Cabaj et al., (2020) used snowfall estimates from CloudSat to calibrate several reanalyses snowfall estimates, including ERA5, within the NESOSIM framework -reducing the spread in snowfall from the chosen reanalyses, although not significantly changing the magnitude of the ERA5 snowfall in the North Atlantic region, where winter snowfall rates are highest overall. On average, ERA5 reports more snowfall over the Arctic basin than 180 what is observed by CloudSat measurements, so the scaling tends to slightly decrease the overall magnitude of the snowfall and the resulting snow depth in NESOSIM (Cabaj et al. 2020). The CloudSat-reanalysis scaling coefficients are now included in the NESOSIM v1.1 code repository.
The other significant code development was the introduction of a new blowing snow loss term. The parameterization of blowing snow lost to leads in NESOSIM v1.0 (Eq. 10 in P2018) has been challenged due to uncertainties 185 around how much snow might be lost to leads under windy conditions, rather than sublimated (lost to the atmosphere) or transported either within or to adjacent grid-cells (Liston et al., 2020). Motivated by this, we introduced an additional blowing snow loss to the atmosphere term, which is a similar function of wind speed and snow in the top 'new' snow layer to the loss-to-leads term, but is not also a function of sea ice concentration: where ℎ ! is the snow depth in the top 'new' snow layer, U is the wind speed, % is the number of seconds in the daily timestep, is the wind action threshold, is the blowing snow loss coefficient and is a new (unitless) atmosphere snow loss coefficient. This parameterization, which provides a simple mechanism for increasing snow loss under given atmospheric 195 conditions independent of sea ice conditions, requires calibration of an additional free parameter, .
Additionally, NESOSIM v1.1 was forced with daily sea ice concentrations from the NASA Climate Data Record (CDR) version 3 (Meier et al., 2017), daily ice drifts from both the NSIDC Polar Pathfinder version 4 dataset (Tschudi et al., https://doi.org/10.5194/tc-2022-39 Preprint. (Lavergne et al., 2010) from September 2019 to April 2021 due to contrasting data availability. As noted in P2018, the impact on snow depth from ice drift forcing is generally second order to snowfall, although this can have first-order impacts at more regional scales. All forcings were regridded to our updated 100 km x 100 km North Polar Stereographic (EPSG: 3413, https://epsg.io/3413) Arctic Ocean model domain.
We recalibrated NESOSIM v1.1 considering the new forcings and model changes described above, by targeting 205 estimates of spring Arctic snow depths derived from Snow Radar data collected during NASA's Operation IceBridge as used in P2018: the snow radar layer detection (SRLD) product (Koenig et al., 2016), the NASA Goddard Space Flight Center (GSFC) empirical threshold based product  and the Jet Propulsion Laboratory (JPL) product (Kwok et al., 2017). Our approach differs from the calibration approach used in P2018, which calibrated NESOSIM v1.0 against Soviet Station drifting station data collected in the 1980s (Warren et al., 1999) then assessed these results against OIB-derived snow 210 depths. Here we choose instead to recalibrate NESOSIM v.1.1 against the spring OIB snow depth data from 2010 to 2015 to provide a more reliable snow depth representation focussed on our contemporary period of interest. We retain, however, the density values for the new 'top' and old 'bottom' layer snow (Table 1) which were derived from the Soviet Station calibration effort.
As noted in P2018 and presented in (Kwok et al., 2017), there is a large spread between the available OIB snow 215 depth products due to various challenges in interpreting Snow Radar data. To account for this large inter-product uncertainty, we develop a 'consensus' gridded OIB spring snow depth product. We take all raw (~7 m along-track resolution) snow depth measurements from the three snow depth retrieval algorithms for a given day, bin them to the 100 km x 100 km NESOSIM v1.1 Arctic Ocean model domain using a simple binning procedure (average of all snow depths in the given grid-cell in each day), then take the median snow depth value at each daily grid-cell across the three OIB products. Quick-Look (QL) snow 220 depths are available for the more recent years (2012 to 2019), using the GSFC waveform fitting approach (https://nsidc.org/data/NSIDC-0708/versions/1). However, it was noted in Kwok et al., (2017) that these estimates tend to exhibit a low bias compared to the other OIB products. A low bias in the GSFC QL product was also shown based on in-situ measurements collected in March 2014 (King et al., 2015). These biases were confirmed in our own analysis comparing our consensus OIB snow depths with the GSFC QL product (2013)(2014)(2015), showing mean biases of ~6 cm (QL thinner than our 225 consensus product, see Supplementary Figure S4), motivating us to exclude these from our model calibration efforts here.
We heuristically calibrated NESOSIM v1.1 using the daily OIB consensus gridded snow depths with the aim of removing the mean bias relative to OIB when using the default NESOSIM v1.0 parameter settings ( Figure 2a). Current work is exploring more automated calibration approaches (Cabaj et al., 2021), but here we were able to find a solution that reduced the mean bias to 0 cm by halving the blowing snow coefficient, tuning the new atmosphere snow loss coefficient, , 230 and extending the model initialization date to September 1 (instead of August 15) as shown in Figure 2b. In the absence of contemporary ground-truth data, we view the initial conditions (either their distribution or the representative start date) as another tuning parameter, within reason. For example, the SnowModel-LG framework presented in Liston et al., (2020) instead initializes with zero snow on August 1st of each year, assuming the snowpack melts out completely each summer based on interpretation of the Warren snow climatology (Warren et al., 1999, W99). As NESOSIM includes no snow melt 235 terms, we prefer instead to initialize later in the year (Sep 1st) and prescribe an expected end of August snow depth. NESOSIM v1.1 was run from 1980 to 2021 and is expected to be updated in future years to enable continued thickness processing from ICESat-2. The output from this v1.1 model framework from 1980-2021 has been archived on Zenodo (https://doi.org/10.5281/zenodo.5164314). The NESOSIM v1.0 output from P2018 was originally released from 2000 to 2015 only but was extended for the 2018/2019 winter to produce snow depths used in the initial P2020 ICESat-2 sea ice 240 thickness processing. To place the ICESat-2 period results in broader context, Figure 3 shows the monthly mean NESOSIM v1.1 snow depth distributions as violin plots, with the recent ICESat-2 years overlaid. During the initial months of winter 250 (September/October), recent years show similar or deeper than average snow, while the middle/end-of-winter months (November to April) show thinner than average snow. The 2019-2020 and 2020-2021 snow depths especially are at or near record low values across most of these months. A more detailed snow depth trend analysis is beyond the scope of this study.
Finally, planning is on-going to release a Quick Look (~2-3-day latency) ATL10 sea ice freeboard product (the current public release latency of ATL10 is ~2-3 months). A Quick Look ATL07 product was recently made available on the 255 NSIDC (https://nsidc.org/data/ATL07QL). Producing near-real time NESOSIM output to derive near-real time thickness estimates is challenging due to the reliance on various forcing datasets from different international groups released at different latencies, so we also produce and describe here the production of a contemporary NESOSIM-derived snow depth/density climatology towards the goal of near real-time thickness assessments with ICESat-2. We take the mean NESOSIM snow depth and density value for each day of the year in each grid-cell averaged across the period September 1, 260 2010, to April 30, 2020. These climatology results (averaged across the Inner Arctic Ocean domain, Figure S3) are also highlighted in Figure 3 and have also been archived on Zenodo (https://doi.org/10.5281/zenodo.5164314).

ICESat-2 sea ice thickness data upgrades
We use the same approach as in P2020 to generate estimates of winter Arctic sea ice thickness and an associated uncertainty estimate. Briefly, thickness is calculated assuming hydrostatic equilibrium and input estimates of sea ice density, snow depth and snow density. The coarse resolution (~100 km) snow depth input estimate, primarily from NESOSIM, is redistributed to the high-resolution (~30 m) ATL10 freeboards using a piecewise functional fit obtained from snow depth and freeboard data collected by NASA's Operation IceBridge mission. Uncertainties are calculated by propagating errors through the hydrostatic equilibrium equation with contributions from random errors (estimates based on previous studies) and systematic errors (estimates based on the spread in applied input assumptions). Small differences in our thickness processing to that 270 presented in P2020 include a bilinear interpolation scheme instead of nearest neighbour to assign NESOSIM data to the ATL10 freeboard segments. Nearest neighbour interpolation was originally used to reduce processing time but introduces unphysical step changes. We also fixed some minor bugs in the freeboard uncertainty calculation and have incorporated the new NSIDC regional mask of the Arctic Ocean (courtesy W. Meier and S. Stewart, NSIDC, see Supplementary Information Figure S3). As in P2020 we use daily estimates of ice type from the European Organization for the Exploitation of 275 Meteorological Satellites (EUMETSAT) Ocean and Sea Ice Satellite Application Facility (OSI SAF, www.osi-saf.org) (Breivik et al., 2012) to classify each segment as either first-year ice (FYI) or multiyear ice (MYI). Ice type information is needed in-part to derive the modified Warren snow depth estimates (see Section 2.2.2. in P2020), so our approach is to assume all ice is MYI unless the OSI SAF product explicitly characterizes the segment as FYI. Thus, in September when OSI SAF does not provide any ice type estimate due to added uncertainties in the end-of-summer retrievals, we assume all 280 our ATL10, and derived thickness data are MYI. The along-track data, both raw and 10 km means, are being made available through the NSIDC (IS2SITDAT4, link to be provided upon completion of peer review) and will be updated as new ATL10 data are released.
In producing the monthly gridded dataset, we use all three strong beams to increase coverage and lower expected uncertainties, compared to the single strong beam used in P2020. The use of all three strong beams was also motivated by 285 the reduction in data coverage in rel003 (and onwards) ATL10 data processing (described in Section 2.1.1 and noted in Figure S1). Our gridding approach is slightly different from the official gridded ICESat-2 freeboard product (ATL20, https://nsidc.org/data/ATL20) as we bin all data within a given month for each grid-cell, as opposed to producing daily gridded composites then monthly gridded composites from the daily gridded data. Our monthly gridded data includes ancillary data variables representative of the mean day of the month for each grid-cell calculated as the mean date of the 290 input ATL10 data, and the number of ATL10 freeboard segments used in the monthly grid-cells to enable sampling bias assessments. The monthly gridded data also includes monthly NOAA/NSIDC Version 4 Climate Data Record (CDR) sea ice concentrations (Meier et al., 2021), the new NSIDC regional mask of the Arctic Ocean (courtesy W. Meier and S. Stewart, NSIDC) and the OSI SAF ice type mask (sub-sampled by ICESat-2 then gridded monthly). The data are projected on to the NSIDC North Polar stereographic grid (EPSG: 3411, https://epsg.io/3411) and binned onto a 25 km x 25 km grid. 295 To further mitigate some of the data gaps and spatial sampling biases, we now also generate smoothed and interpolated variables of freeboard, snow depth and thickness. Our approach is as follows: i) set monthly mean grid-cells to zero where the monthly CDR concentration is valid and less than 15%, ii) apply linear interpolation using Delaunay triangulation (https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.LinearNDInterpolator.html) on all grid- cells to interpolate between valid grid-cells, iii) smooth using a Gaussian filter with a kernel width of 0.5 standard deviations 300 in x and y directions; iv) mask all grid-cells greater than 50 km away from grid-cells containing data in the original monthly gridded dataset using a KDTree algorithm, v) mask interpolated/smoothed data where the monthly NOAA/NSIDC CDR concentration is <50%. We expect that future work will explore more sophisticated interpolation procedures and blending with other thickness datasets (discussed more in the summary).
The initial version of this monthly gridded dataset described in P2020 (no interpolated variables, no region mask or 305 sea ice concentration data, based on rel004 ATL10) was made available through the NSIDC (IS2SITMOGR4 Version 1, https://nsidc.org/data/IS2SITMOGR4). An updated dataset using the updated rel005 ATL10 data and including the interpolated variables, updated NSIDC Arctic region mask and CDR sea ice concentrations shown in this study are being made available as a new version 2 (v2) release of the IS2SITMOGR4 dataset which we expect to keep updated along with the along-track product (IS2SITDAT4) as new ATL10 data are released. 310

CryoSat-2 sea ice thickness estimates
Following P2020 we compare our winter Arctic sea ice thickness estimates with those generated from the European Space Agency (ESA) CryoSat-2 mission from four different groups: NASA's Goddard Space Flight Center (GSFC, Kurtz and Harbeck, 2017), the Jet Propulsion Laboratory (JPL, Kwok and Cunningham, 2015), the Center for Polar Observation and Modelling (CPOM, Laxon et al., 2013;Tilling et al., 2018) and the Alfred Wegener Institute (AWI, Hendricks and Ricker, 315 2016). As the different products make different assumptions regarding snow loading and sea ice, for these comparisons we use the same snow loading assumptions in our ICESat-2 thickness processing to generate direct thickness comparisons (instead of the NESOSIM based thickness estimates). As in P2020 we re-grid the monthly gridded CS-2 estimates to the NSIDC 25 km x 25 km North Polar Stereographic grid using a simple nearest neighbor interpolation scheme and compare these with our gridded ICESat-2 sea ice thickness estimates that have been produced using the same snow loading and ice 320 density assumptions as the given CS-2 product, as summarized in Table 2 in P2020 (e.g., modified versions of the Warren climatology snow loading). Note that in the CryoSat-2 comparisons we produce monthly gridded data using strong beam #1 only to be consistent with the results shown in P2020.
Modified versions of the Warren snow depth climatology (mW99, Warren et al., 1999) were used by all four of these CryoSat-2 thickness products (it is worth noting that most of these groups are actively working on incorporating more 325 sophisticated snow loading models). P2020 noted strong differences between NESOSIM v1.0 and mW99 (mW99, snow depths halved over first-year ice, see Supplementary Figures S2 to S4 in P2020). Generally, snow depths are similar over the thicker multiyear ice, but mW99 is thinner later in the year, due primarily to the reduced snow over first-year ice. Figure 3 shows the mean Inner Arctic Ocean snow depth from mW99 (snow depths halved using observed OSI SAF ice types from 2010 to 2019), showing similar values to NESOSIM in October but thinner mW99 snow in April. Again, we use the same 330 snow loading when producing comparison ICESat-2 thickness estimates, but it is worth noting that these thickness estimates will be different to our NESOSIM-derived thickness product due to these differences in snow.

PIOMAS sea ice thickness estimates
In this study we additionally compare our gridded winter Arctic sea ice thickness estimates with those generated from the Pan-Arctic Ice-Ocean Modeling and Assimilation System (PIOMAS, v2.1;Zhang & Rothrock, 2003). PIOMAS is an ice-335 ocean model that generates estimates of sea ice thickness, constrained predominantly by the assimilation of sea ice concentration and sea surface temperature. PIOMAS data is commonly used in the sea ice community for longer-term assessments of Arctic sea ice thickness change. PIOMAS ice thicknesses estimates have been shown to exhibit differences on the order of tens of centimeters compared to satellite-derived estimates, although this depends strongly on the season and region analyzed (Schweiger et al., 2011;Zygmuntowska et al., 2014;Petty et al., 2018b). 340

ATL10 freeboards, NESOSIM snow loading and sea ice thickness distributions
In Figure 4 we show probability distributions of winter Arctic sea ice freeboard calculated using rel002 (as used in P2020) through to rel005 ATL10 data. We  cm) and a 3.4 cm increase in April 2019 (35.9 cm to 39.3 cm). In contrast, the rel003 to rel005 freeboard distribution differences across these three months are small or negligible (<0.4 cm). Rel005 generally shows the highest freeboards from the four releases. As discussed earlier, this was largely expected due to the significant algorithm change in rel003 (Kwok et al., 2021d) and the lack of significant algorithm changes related to freeboard derivation in rel004 and rel005. The 1.2 cm 355 mean freeboard reduction in November 2018 is due to a more significant primary freeboard mode and a less significant secondary peak in the rel003 to rel005 freeboard distributions, while January 2019 and April 2019 distributions exhibit a clear increase in the unimodal freeboard in rel003 onwards. Kwok et al., (2021d) analysed gridded freeboard distributions in January 2019, June 2019, October 2019 and found 3 cm, 1 cm and 2 cm increases respectively between rel002 and rel003, in-line with the differences observed here. As discussed in Section 2.1.1 and demonstrated in Figure S1, the different 360 releases also include changes in coverage, especially between rel002 and rel003, which may influence these differences along with changes in the freeboard determination algorithm.  Figure 4 (right column) also shows an analysis of the inter-beam differences across the three strong beams for the same time periods and Inner Arctic Ocean region for rel005 data only. The inter-beam differences are small (<1 cm), similar to the rel003-rel005 differences. Each strong beam is separated by ~3 km across track, so more significant differences are 365 expected at local scales, however an in-depth analysis of spatial length-scales is beyond the scope of this study. We instead note that at basin/monthly scales, the beams provide similar freeboard distributions despite the small differences in coverage (as discussed in Section 2.1.1 and highlighted in Figure S2), increasing our confidence in using all three beams to extend coverage across the Arctic.  The difference in mean sea ice thickness across the three months and NESOSIM frameworks is greatest in November 2018: Nv1.1 thicknesses are 10 cm thinner than Nv1.0 while the Nv1.1clim thicknesses show a higher primary peak but longer tail.
The thickness differences across the three NESOSIM frameworks are 4 cm in January 2019 and 2 cm in April 2019. In 385 general, the impact on thickness from the choice of NESOSIM framework (including the NESOSIM climatology) is less significant than the impact from rel002 to rel003 freeboard changes.

Comparisons with CryoSat-2
In Figure 6 we show comparison statistics (correlation coefficient, mean bias and standard deviation) of monthly gridded sea ice thickness estimates using both rel002 and rel003 ATL10 freeboards with those produced from ESA's CryoSat-2 using 390 the same input assumptions, i.e., modified versions of the Warren climatology for snow loading (See Table 2 in P2020). As in P2020 we only use strong beam #1 and mask all data below 0.25 m in both datasets and outside of an Inner Arctic Ocean domain before producing these comparisons.
In general, the agreement between ICESat-2 and CryoSat-2 using rel003 ATL10 freeboards is improved compared to those based on rel002 ATL10 freeboards in terms of the correlation coefficient, mean bias and standard deviation across  Figure 6. More work is needed to better reconcile these datasets and assess sources of bias (as discussed more in the summary) but these results represent an encouraging initial development in terms of the ICESat-2 sea ice freeboard and thickness product development. 405

Gridded Arctic sea ice thickness and comparisons over the last three winters.
In Figure 7 we show an example output of our updated monthly gridded ICESat-2 winter Arctic sea ice thickness product (IS2SITMOGR4 version 2, v2) for April 2021 using the latest default thickness processing configuration (rel005 ATL10 and NESOSIM v1.1 snow loading). Spatial coverage is high across all months (not shown) despite the concerns expressed in Kwok et al., (2021d) related to reduced lead/sea surface height segments and thus freeboard determination. This was partly 410 mitigated by our use of 3 (strong) beams (coverage changes discussed in Section 2.1.1). Despite our use of three strong beams, there are still large regions of missing data in our monthly gridded dataset, e.g., the missing data in the Laptev/East Siberian Sea shown in Figure 7 despite the monthly CDR ice concentrations showing concentrations greater than 50% in that same region. Data drop-out is often caused by the presence of clouds and the resultant atmospheric scattering impacts on ATLAS retrievals. Our interpolated/smoothed variables of freeboard, snow depth and thickness data (Figure 7i-j) do not 415 substantially increase coverage in these regions, which was in-part by design to avoid over-extrapolation of our thickness estimates. In general, the monthly gridded data gaps are limited but should be considered when using these data to assess regional and basin-scale thickness variability.
Our following analysis of the monthly gridded ICESat-2 winter Arctic sea ice thickness data is a subset of the analysis presented in an online interactive Jupyter Book: https://nicolekeeney.com/icesat2-book. The Jupyter Book consists 420 of a series of Jupyter Notebooks that provide all code and analysis output (written in the Open-Source Python programming language) for demonstrating and sharing our thickness analysis workflow. The development of the Jupyter Book was motivated by the desire for transparency and the broader goals of facilitating more open science, but also the desire to provide a simple mechanism for interested users to explore regions and time periods not shown here. For example, the Jupyter Book allows users to adapt the code interactively (either locally or using Binder) to select months and regions of 425 interest to explore characteristics of this dataset beyond the core figures we show here.
In Figure 8 we show the seasonal evolution of winter freeboard, snow depth/density and sea ice thickness from our IS2SITMOGR4 v2 monthly gridded dataset. The results shown in Figure 8 are again restricted to our Inner Arctic Ocean domain using the included NSIDC Arctic region mask ( Figure S3) to simplify interpretation and avoid regions of higher 430 uncertainty within the more peripheral seas of the Arctic. The data within this domain in September and October is generally lower concentration (40-60% on average, Figure 8e), than the proceeding months (~90% in November and near 100% in December through April), so changes in early winter are still strongly influenced by the changing coverage of sea ice as the ice pack refreezes. We thus mainly focus on analysing November to April changes, the period in which we also have full monthly data available across all three winters. results in a similar mean thickness. It is also worth noting that the fraction of multi-year ice is inversely related to the thickness rankings -i.e., the 2018-2019 winter shows the lowest mean fraction of multi-year ice in this three-year period, but also shows the highest freeboard, snow depths and thickness. Three years is not a long enough record to establish true 460 relationships, but the results highlight the potential pitfalls of inferring thickness from ancillary quantities such as freeboard or multi-year ice fraction if one is interested in tracking interannual changes.
To assess the spatial distribution of the winter changes discussed above, Figure 9 shows maps of winter mean (November to April) freeboard, snow depth and thickness from IS2SITMOGR4 v2, while Figure 10  The strong negative Barents Sea snow depth anomalies mitigate much of the negative freeboard anomalies observed by ICESat-2, meaning the impact on the 2020-2021 thickness is strong negative anomalies throughout much of the Central Arctic region and into the Chukchi Sea. The low freeboards within the Barents Sea region and the increased potential for surface flooding in this region (Granskog et al., 2017), a process which is not currently simulated by NESOSIM, means those results should also be treated with caution (this region is mostly excluded from our Inner Arctic Ocean domain, Figure  485 S3). Monthly maps and anomalies of all IS2SITMOGR4 v2 variables have been generated and provided in the relevant Jupyter Book page (https://nicolekeeney.com/icesat2-book/sea_ice_characteristics).
In Figure 11 we show a comparison of these three winters of Inner Arctic Ocean sea ice thickness (and associated uncertainties) with thickness estimates derived from PIOMAS. As opposed to the CryoSat-2 comparisons, we focus here on basin-scale comparisons to assess the broad level of agreement when using both datasets for assessing the winter seasonal 490 cycle and interannual thickness variability -arguably the primary use of PIOMAS data to-date. In general, both seasonal timeseries show good agreement, however PIOMAS is generally thicker than ICESat-2, especially by the end of winter in all three years of our analysis (20 to 40 cm thicker). These differences are within the spread of our ICESat-2 systematic uncertainty estimates, however. Spatial difference plots are generated and provided in the Jupyter Book (https://nicolekeeney.com/icesat2-book/pio_vs_is2) which highlight the stronger disagreement at regional scales, e.g., 495 PIOMAS not simulating the thicker ice north of the Greenland and CAA coasts shown in our ICESat-2 data.

Key drivers of winter thickness differences
The regional anomaly maps allude to a strong ice type dependency as some of the more significant winter anomalies are observed within the thicker/older ice of the Central Arctic. To explore this further, Figures 12 and 13 show the seasonal timeseries of IS2SITMOGR4 v2 but delineated by ice type (data still masked outside the Inner Arctic Ocean 500 domain). Figure 12 shows the mean seasonal time series of regions identified as first-year ice (FYI) only. The differences in FYI freeboard and snow depth in November/December are small (<1 cm) across the three winters, with the higher 2018-2019 FYI freeboard and snow depth only appearing later in the season compared to the 'all ice' analysis. The resultant FYI sea ice thickness winter timeseries comparison is notable for its consistency across the three winters (interannual thickness differences < 15 cm across all months). Figure 13 shows the mean seasonal time series of regions identified as multiyear ice 505 (MYI) only. The interannual MYI differences across most variables are higher than the FYI differences. The 2018 November MYI freeboards are 5 cm higher than the 2019 and 2020 Novembers. This difference largely persists through the winter until February onwards when the 2019-2020 freeboard increases and reduces the interannual spread, driven by the coincident strong increase in 2019-2020 snow depths. The result is MYI thickness that exhibits similar seasonal cycles across the three winters but with differences of 10 to 50 cm that largely persist across the three winters. 510 These ice type differences align with our general understanding of winter sea ice growth -thinner ice is more responsive to atmospheric forcing and can thicken rapidly due to its reduced insulation (the negative feedback of ice growth) so small differences in the thickness of thin FYI at the start of winter are not expected to be good predictors of end-of-winter thickness (Petty et al., 2018b). Conversely, MYI is significantly thicker at the start of winter, meaning thickness anomalies are more likely to persist through winter as the ice is more insulated and less sensitive to atmospheric forcing. The 515 differences in MYI thickness at the start of our three winters appears to provide a strong control on the total (combined MYI and FYI) winter thickness anomalies across all months, albeit in this limited record. Previous studies based on CryoSat-2 derived Arctic sea ice thickness estimates have highlighted the significant role of variable summer conditions in determining start of winter ice thickness anomalies and thus total winter thickness (and volume) anomalies (Tilling et al., 2015;Kwok, 2015). More specifically, a sharp increase in the start-of-winter 2013 Arctic thickness/volume was related to reductions in 520 the duration of the summer melt season (Tilling et al., 2015) and also to dynamically driven convergence of ice within the Central Arctic (Kwok, 2015). The observed positive autumn 2013 thickness anomaly persisted through winter months, as in our 2018/2019 results. We do not seek to provide a similar level of analysis in this study as our primary goal was to highlight and describe this new thickness dataset, but the agreement with our prior physical understanding is encouraging.
To better understand some of the more regional differences, the spatial thickness maps in Figure 10 have been 525 overlaid with winter mean ice drifts from the monthly OSI SAF global low resolution ice drift product (Lavergne et al., 2010). In general, the mean circulation across these three winters are similar -featuring anti-clockwise Beaufort Gyre circulations and Transpolar drifts, but with some key differences. For example, ice drifts through the southern Beaufort Sea in 2018-2019 winter were stronger than the 2019-2020 and 2020-2021 drifts, which is likely associated with the positive freeboard/snow depth/thickness anomalies observed in the Chukchi Sea and negative anomalies in the Beaufort Sea in 2018-530 2019. Disentangling cause from effect is challenging as ice drift is strongly influenced by the sea ice conditions (Petty et al., 2016), however the strong association between drift patterns and thickness anomalies is again encouraging. Stronger ice drift anomalies are apparent when assessing monthly (not seasonal) differences, which can be explored more in the relevant Jupyter Book page (https://nicolekeeney.com/icesat2-book/sea_ice_characteristics).
Finally, to briefly explore any possible winter atmospheric drivers of these differences, Figure 14 shows

Summary
In this study we provided an impact assessment of upgrades to the input data used to produce ICESat-2-derived winter Arctic sea ice thickness estimates shown in Petty et al., (2020), and an extended analysis of the upgraded monthly gridded winter Arctic thickness dataset across the three winters profiled since the launch of ICESat-2 in September 2018.
Input data upgrades include the ICESat-2 ATL10 freeboards (Release 002 to 005, rel002 to rel005) and NASA 550 Eulerian Snow On Sea Ice Model (NESOSIM, version 1.0 to version 1.1) snow loading. A key change in ATL10 data was the removal of misclassified leads from the determination of sea surface and thus freeboard in rel003. This was thought to be the primary cause of the increase in freeboard observed in January 2019 and April 2019 in rel003 data compared to rel002, together with the more significant primary mode of lower freeboards in November 2018 rel003 data. The rel003-derived monthly gridded winter Arctic ice thickness data show improved comparisons with thickness estimates produced from 555 ESA's CryoSat-2 using the same input assumptions across all 2018-2019 winter months (lower mean biases and standard deviations, higher correlations) compared to rel002-derived estimates. Later releases of ATL10 (rel004 and rel005) involved only minor changes to the freeboard algorithms and thus exhibit less significant changes in the observed freeboard distributions. The different releases also show slight differences in ATL10 data coverage, due primarily to the changes https://doi.org/10.5194/tc-2022-39 Preprint. Discussion started: 16 February 2022 c Author(s) 2022. CC BY 4.0 License. associated with dark lead usage, but also the inclusion/filtering of satellite calibration scan data. Our updated (version 2) 560 monthly gridded winter Arctic sea ice thickness dataset now utilizes all three strong beams and includes new interpolated/smoothed data variables to mitigate these coverage issues.
Our upgraded version of NESOSIM (version 1.1) presented in this study includes a new wind-driven atmosphere snow loss term, CloudSat-scaled ERA5 snowfall forcing (Cabaj et al., 2020) and some more minor bug fixes. NESOSIM v1.1 was also re-calibrated (heuristically) using spring Arctic snow depth estimates obtained by NASA's Operation 565 IceBridge airborne mission (a gridded consensus product derived in this study). NESOSIM v1.1 generally shows similar snow depths to NESOSIM v1.0, resulting in a less significant impact on Arctic winter sea ice thickness compared to the rel002 to rel003 freeboard changes. A new NESOSIM v1.1 climatology (from 2010-2020 daily means) was introduced towards the production of a more near-real time thickness retrievals from the forthcoming Quick Look ATL10 freeboard dataset. 570 Finally, we presented estimates of winter Arctic sea ice thickness from this updated monthly gridded winter Arctic sea ice thickness dataset (IS2SITMOGR4 v2) over the past three winter seasons of data collection (November 2018 -April 2021, September 2019 -April 2020 and September 2020 -April 2021). Our results showed clear differences in mean winter Arctic sea ice thickness within our Inner Arctic Ocean domain across the three winters profiled, due primarily to differences in the multiyear ice thickness across the three winters (multiyear ice thinning of 10 to 50 cm each year across the three 575 winters analysed). Interannual changes in snow depth provide significant regional/monthly impacts on our thickness resultsmitigating some (or in some cases all) of the impact from interannual differences in Arctic winter freeboards observed by ICESat-2. Our results provide further evidence of the importance of accurate snow representation when assessing interannual variability in winter Arctic sea ice thickness from satellite altimetry (Bunzel et al., 2018;Mallett et al., 2021). Specific regional thickness anomalies, e.g. in the Southern Beaufort and Chukchi seas, were associated with interannual ice drift 580 anomalies. Our mean Inner Arctic Ocean thickness estimates showed good agreement with those generated from the PIOMAS v2.1 reanalysis in terms of the seasonal cycle and interannual differences, although more significant differences were noted at regional scales.

Future work
ICESat-2/ATL10: Work is still on-going to re-introduce dark leads to the sea surface and freeboard algorithm in 585 ATL10, which requires a new filter to skilfully discriminate dark lead segments (low photon rate) from segments with photon attenuation driven by the presence of clouds. The variable properties of clouds and their impact on photon attenuation, together with the limited availability of coincident imagery for validation (as used in Petty et al., 2021) makes this development challenging. An additional near-term goal related to ATL10 is the plan to utilize all six beams (or at least the three strong beams) concurrently to produce two-dimensional interpolated fields of sea surface height, as opposed to the 590 independent beam processing currently utilized. However, residual height biases of several centimetres are still observed between the beams as of Release 005 (updated from the analysis shown in Bagnardi et al., 2021, not shown), hindering this development. More sophisticated sea surface interpolation methods should also be explored (Landy et al., 2021). Additional algorithm development efforts are on-going towards the next data release (Release 006), expected sometime later in 2022.

NESOSIM:
Work is on-going to utilize a Markov Chain Monte Carlo to further calibrate NESOSIM v1.1 -595 providing better constraints on the free parameters and also more robust uncertainty estimates of snow depth and density from this model framework (Cabaj et al., 2021). This approach benefits from the low computational cost of NESOSIM, allowing thousands of model simulations to be generated across plausible model parameter space. Additional physical upgrades are still desired, e.g., the introduction of a snow melt parameterization to extend NESOSIM through summer.
However, additional reliable ground-truth data at regional/basin-scales are needed to calibrate and validate such development 600 activities.
ICESat-2-derived sea ice thickness: Our primary focus of the three-winter thickness assessment was the monthly gridded winter Arctic thickness dataset (IS2SITMOGR4). However, raw (and 10 km smoothed) along-track data at the segment resolution of ATL10 (~20 m) are also available (IS2SITDAT4, data shown in Figure 5), which provide higher fidelity information regarding the sea ice state than the monthly gridded estimates. Work is currently on-going to assess the 605 winter Arctic sea ice thickness distribution from these data including comparisons with model-based estimates. Continued refinement and/or redevelopment of the snow redistribution scheme is expected. We also hope to combine these data with new ICESat-2-derived floe size estimates (Petty et al., 2021) towards a joint floe size-thickness distribution. This dataset is more computationally demanding but increasing access to high performance computing environments (e.g., cloud compute platforms) should help increase its usability. The extension of NESOSIM through summer months will help enable summer 610 preliminary production of summer Arctic thickness estimates.

Sea ice thickness reconciliation:
The improved correspondence between our ICESat-2 derived estimates of winter Arctic sea ice thickness and those generated from ESA's CryoSat-2 are encouraging. There are clear advantages (and disadvantages) from estimating sea ice thickness from either radar or laser altimetry, which need to be better considered and utilized for constraining total Arctic (and eventually Antarctic) sea ice volume. Radar altimeters, e.g., CryoSat-2, are highly 615 sensitive to leads and are unaffected by clouds, providing benefits to both the quality and coverage of data collected. In contrast, laser altimeters (e.g., ICESat/ICESat-2) generally provide higher resolution data and obtain more precise estimates of the snow-covered ice surface height (and thus total freeboard) compared to the arguably less distinct/certain ice-snow interface height (and thus ice freeboard) obtained by typical radar altimeters. The effective radar penetration depth at Ku/Kaband is generally considered to come from the ice-snow interface although recent studies continue to challenge this (Nandan 620 et al., 2017, King et al., 2018. In both cases, uncertainties in the derived freeboard estimates are combined with uncertainties in the various input assumptions (snow loading, sea ice density) to provide total thickness uncertainty estimates.
Constraining the various input uncertainties and residual biases remains challenging, which points to the need for improved exploitation of existing ground-truth data and further field and airborne campaigns considering the fast-changing Arctic.
Regardless, improvements to the underlying freeboard algorithms and input assumptions are urgently needed as we seek to 625 reconcile these datasets and hopefully move towards multi-sensor thickness assessments (increasing coverage and data quality). Planning is underway for a coordinated intercomparison exercise around new semi-synchronous along-track measurements available since the CRYO2ICE orbit alignment (https://earth.esa.int/eogateway/missions/cryosat/cryo2ice).

Code availability
Our analysis of the monthly gridded winter Arctic thickness data described above (Figures 8-14) have been summarized and 630 made available through an online Jupyter Book (https://nicolekeeney.com/icesat2-book/home.html). Interested users are able to view the various figures (and additional plots not shown here), view the code used to generate them, and also interactively run the analysis locally or online (using the associated Binder links, https://mybinder.org), e.g., changing the domain and/or months of interest. It is our expectation that this Jupyter Book will be updated as new IS2SITMOGR4 data are created and made public to enable continued assessments of winter Arctic thickness change. A version tagged version of this Jupyter 635 Book will be obtained and archived on Zenodo on completion of peer review.
The original sea ice thickness processing code presented in Petty et al., (2020) is available on GitHub (https://github.com/akpetty/ICESat-2-sea-ice-thickness, all in the open-source language Python). We plan to update this 640 using the small upgrades made to our processing chain on completion of this peer-review.

Data availability
The monthly gridded winter Arctic sea ice thickness data derived in this study (IS2SITMOGR4, version 2) is being made available through the National Snow and Ice Data Center (NSIDC) (https://nsidc.org/data/IS2SITMOGR4, currently showing version 1 from Petty et al., 2020). The along-track (raw and 10 km mean) Arctic sea ice thickness estimates are 645 also in the process of being ingested and made publicly available through the NSIDC (IS2SITDAT4, we expect the data portal will be made available during peer review).

Author contributions 665
AP led the study, produced the updated NESOSIM v1.1 framework/output and ICESat-2 thickness estimates (along-track and gridded). NK produced the Jupyter Book with assistance from AP and helped produce the new interpolated thickness variables. AC and PK helped AP produce the updated NESOSIM framework and output. MB generated the release and beam coverage assessments and provided key input on ATL10 upgrades. AP wrote the manuscript with input/edits provided from all authors. Wind action threshold (m s -1 ), ⍵ 5 5 Blowing snow coefficient (s -1 ), 2.9 x 10 -7 1.45 x 10 -7 Atmosphere loss coefficient, N/A 0.15 Wind packing coefficient, ⍵ (s -1 ) 5.8 x 10 -7 5.8 x 10 -7

Initial conditions
Start date August 15th September 1st