Advances in altimetric snow depth estimates using bi-frequency SARAL and CryoSat-2 Ka–Ku measurements

Although snow depth on sea ice is a key parameter for sea ice thickness (SIT) retrieval, there currently does not exist reliable estimations. In the Arctic, nearly all SIT products use a snow depth climatology (the modified Warren-99 climatology, W99m) constructed from in situ data obtained prior to the first significant impacts of climate change. In the Antarctic, the lack of information on snow depth remains a major obstacle in the development of reliable SIT products.

In this study, we present the latest version of the altimetric snow depth (ASD) product computed over both hemispheres from the difference of the radar penetration into the snow pack between the Ka-band frequency SARAL/Altika and the Ku-band frequency CryoSat-2. The ASD solution is compared against a wide range of snow depth products including model data (Pan-Arctic Ice-Ocean Modelling and Assimilation System (PIOMAS) or its equivalent in the Antarctic the Global Ice-Ocean Modeling and Assimilation System (GIOMAS), the MERCATOR model, and NASA's Eulerian Snow On Sea Ice Model (NESOSIM, only in the Arctic)), the Advanced Microwave Scanning Radiometer-2 (AMSR2) passive radiometer data, and the Dual-altimeter Snow Thickness (DuST) Ka–Ku product (only in the Arctic). The ASD product is further validated in the Arctic against the ice mass balance (IMB) buoys, the CryoSat Validation Experiment (CryoVEx) and Operation Ice Bridge's (OIB) airborne measurements. These comparisons demonstrate that ASD is a relevant snow depth solution, with spatiotemporal patterns consistent with those of the alternative Ka–Ku DuST product but with a mean bias of about 6.5 cm. We also demonstrate that ASD is consistent with the validation data: comparisons with OIB's airborne snow radar in the Arctic during the period of 2014–2018 show a correlation of 0.66 and a RMSE of about 6 cm. Furthermore, a first-guess monthly climatology has been constructed in the Arctic from the ASD product, which shows a good agreement with OIB during 2009–2012. This climatology is shown to provide a better solution than the W99m climatology when compared with validation data. Finally, we have characterised the SIT uncertainty due to the snow depth from an ensemble of SIT solutions computed for the Arctic by using the different snow depth products previously used in the comparison with the ASD product. During the period of 2013–2019, we found a spatially averaged SIT mean standard deviation of 20 cm. Deviations between SIT estimations due to snow depths can reach up to 77 cm. Using the ASD data instead of W99m to estimate SIT over this time period leads to a reduction in the average SIT of about 30 cm.

1 Introduction

Since the launch of CryoSat-2 (CS-2) in 2010 (Wingham et al., 2006; Parrinello et al., 2018), sea ice thickness (SIT) observations are routinely derived from altimetric measurements. The principle is to measure the fraction of the sea ice above the sea level, called the sea ice freeboard, from differences between the heights in leads (cracks in the ice referring to the local sea level) and the heights of the ice floes (Laxon et al., 2003). By integrating such sea ice freeboard estimations in the hydrostatic equilibrium equation, several SIT products have been computed (e.g. Laxon et al., 2013; Kwok and Cunningham, 2015; Guerreiro et al., 2017; Paul et al., 2018; Landy et al., 2019; Laforge et al., 2020).

Among the parameters involved in the SIT calculation, snow depth (sd) over sea ice is one of the most significant contributors adding to the overall SIT uncertainty (e.g. Giles et al., 2007; Zygmuntowska et al., 2014; Guerreiro et al., 2016). For example, it is necessary to account for the snow loading (Laxon et al., 2013) and for the decrease in altimetric radar speed as it penetrates into the snow pack (Kwok and Cunningham, 2015; Mallett et al., 2020). In fact, variabilities in snow properties affect radar signals and then the snow depth measurements. More generally, snow cover has strong impacts on the sea ice (e.g. Massom et al., 2001; Powell et al., 2005; Bin et al., 2008; Sturm and Massom, 2009; Ricker et al., 2014) that affects the entire climate system (e.g. Ingram et al., 1989; Ledley, 1991; Eicken et al., 1995; Singarayer et al., 2006). Because of its high albedo and a low thermal conductivity, the snow regulates the transfer of solar heat energy penetration across the ice–ocean interface (e.g. Grenfell and Maykut, 1977; Sturm et al., 1997). It acts as an insulator, slowing down sea ice melt in summer and slowing down sea ice growth in winter (e.g. Perovich et al., 2003; Sturm and Massom, 2016). Such processes of sea ice formation and melting govern sea ice physical and chemical properties that impact the biological processes in sea ice (Van Leeuwe et al., 2018). The vertical distribution of light under sea ice that controls biological processes and biogeochemical fluxes is also strongly linked with the snow depth (e.g. Perovich, 2007; Arndt and Nicolaus, 2014; Arndt et al., 2017). In addition snow accumulated over sea ice from precipitation represents a tank of freshwater likely to be carried into the ocean. Recent increases in seasonal ice have promoted the amount of snow water discharge in the ocean, which impacts the freshwater budget (Andersen et al., 2019; Overland et al., 2019). Snow cover also modifies surface roughness that impacts the air–ice drag coefficient and transfer coefficients of latent and sensible heat fluxes (Andreas et al., 2005).

Despite evidence of its major importance, snow cover over sea ice is still not sufficiently understood. Until now, nearly all SIT products in the Arctic are computed using the snow depth of the modified Warren-99 climatology (W99m; Warren et al., 1999). W99m is mainly constructed using snow depth observations obtained from the Soviet North Pole drifting stations based in the central Arctic, collected during the last century (1950–1980s). Considering the ongoing rapid modifications of the Arctic due to climate change, several studies have shown these data to be outdated, even with the modification of snow depth reduction (by 50 %) applied over the seasonal ice zone (Kurtz and Farrell, 2011; Kern et al., 2015; Kwok et al., 2011). The climatology of Forsström et al. (2011) provides snow depths in areas outside the central Arctic where the W99m climatology is not working properly. More recently, Shalina and Sandven (2018) have created an improved snow depth climatology taking into account snow depth observations obtained during campaigns over seasonal ice. However, it is still based on data collected mainly in the 1960s, 1970s and 1980s. Apart from the data used to construct these climatologies, several campaigns have provided snow depth measurements in the Arctic. Between 1997 and 1998, the SHEBA (Surface HEat Budget of the Arctic Ocean) project has highlighted the complex temporal and spatial snow depth variations (Perovich et al., 1999; Sturm et al., 2002). The ice mass balance (IMB) buoys, originally deployed by the US Army CRREL-Dartmouth Mass Balance Buoy Program during the SHEBA project have measured snow depths since the 2000s. It allows for monitoring changes in the sea ice volume in key areas of the Arctic (Richter-Menge et al., 2006). Since 2003, the CryoSat Validation Experiment (CryoVEx) campaigns (e.g. Haas et al., 2006; Helm et al., 2006; Skourup et al., 2013; Hvidegaard et al., 2006) have provided data with the main goal of investigating radar penetrations into ice and snow cover. Their measurements include bi-frequency altimetry snow depth estimations as of 2017 (see Sect. 3.3). In 2003, the AMSR-Ice03 campaign (Sturm et al., 2006) was carried out in the Beaufort Sea region to validate the Advanced Microwave Scanning Radiometer for the Earth Observing System (AMSR-E, from 2002 to 2011; Kelly, 2009) passive microwave snow depth data. The same year, an equivalent mission (ARISE) took place in East Antarctica (Massom et al., 2006). Since 2009, and until their final campaign in 2020, the Operation Ice Bridge (OIB; see Sect. 3.3) campaigns have provided airborne snow depth measurements each year between March and April in the Arctic (Kurtz et al., 2013; Koenig et al., 2010) using a frequency modulated continuous wave (FMCW) snow radar (Kurtz and Farrell, 2011). Note also that Webster et al. (2014) assess spring snow depth distribution on the Arctic sea ice from 2009–2013 airborne OIB radar observations. In parallel, in situ snow thickness data were collected in 2009 at the Danish GreenArc 2009 ice camp in the north of Greenland (Farrell et al., 2011). Between February and March 2015, snow depth was measured in the Nansen Basin during the N-ICE expedition (Granskog et al., 2016), demonstrating among others the impact of heavy precipitation on thin ice growth (Merkouriadi et al., 2017). From September 2019 to October 2020, the international MOSAIC expedition monitored the central Arctic (Shupe et al., 2020) providing, among others, data on snow precipitation, snow water equivalent and Magnaprobe snow depth measurements (e.g. Munoz-Martin et al., 2020). Hence, several campaigns have been conducted during the last two decades in the Arctic, all of which have provided crucial information and proved invaluable as reference data.

Compared to the Arctic, snow depth measurements, and more generally sea ice in situ data, are very limited in the Antarctic. Actually, several CryoVEx and OIB missions have been conducted in the Antarctic, but the data are still only available in raw L1b level and are difficult to interpret. To the best of our knowledge, only the meereis data portal (https://seaiceportal.de/en/, last access: March 2021) provides snow buoy data of snow accumulation (Grosfeld et al., 2016) (for both hemispheres), but the comparison with snow depth is not direct. It is, for instance, limited by the flooding of ice floes due to the heavy snow loading occurring in the Antarctic. Currently, the main expedition providing snow depth data in the Southern Hemisphere is the ASPeCt programme (Worby et al., 2008a), established in 1996 to model the role of the Antarctic sea ice in the coupled atmosphere–ice–ocean system. Recent ASPeCt standardised data (Kern, 2020) cover the period 2002–2019 with an amount of information that is too sparse for a reliable assessment. Future validations in the Antarctic will also benefit from the recent AWI IceBird campaigns carrying a snow radar since 2019 (https://www.awi.de/en/science/climate-sciences/sea-ice-physics/projects/ice-bird.html, last access: March 2021).

At a basin scale, several sea ice models produce snow accumulation estimations from atmospheric reanalyses. These are converted into snow depth over sea ice by considering sea ice drifting and thermodynamic transfers of the coupled sea ice and snow system (Blanchard-Wrigglesworth et al., 2015; Lecomte et al., 2011). Likely due to biases in precipitation inputs and omission of snow processes, these global sea ice models hardly reproduce consistent on-ice snow depths (e.g. Holland et al., 1993; Serreze et al., 2000; Déry and Tremblay, 2004; Leonard and Maksym, 2011; Blazey et al., 2013). Notz (2012) states that an improved representation of the snow on sea ice in models is the most urgent task, and recently, Kaminski et al. (2018) argued that assimilating snow depth products could considerably increase sea ice model performances.

Therefore, a solution allowing us to provide large-scale information on the spatial and temporal variation in snow depth and to improve the performance of seasonal sea ice forecasts is to use satellite data. For now, only a few datasets exist. The first snow depth estimates come from passive radiometric data. Among the various sensors that have been deployed (e.g. the SMMR, Chang et al., 1987, or the SSM/I, Grody, 1991), the most relevant estimations of snow depth over sea ice have been calculated from spectral vertical gradients of brightness temperatures (Comiso et al., 2003; Worby et al., 2008b) of AMSR-E and its successor AMSR2 (Lee et al., 2015). In spite of the relevant long-term availability of data, the commonly used snow depth retrieval algorithm is not adapted for multi-year ice (MYI) zones (Brucker and Markus, 2013). Consequently, AMSR's official data in the Arctic are only calculated in first year ice (FYI) regions. As far as we know, the AMSR2 product version of the University of Bremen (AMSR2B; Rostosky et al., 2018) and very recently the SPICES project snow depth retrieval method (Mäkynen et al., 2020) are the only attempts to compute snow depths in the Arctic MYI zones. The dataset developed from the work of Rostosky et al. (2018) is available but has the inconvenience of being re-calibrated from OIB data, which limits their availability to the months of March and April. With the same disadvantage of depending on OIB, Maaß et al. (2013) and Zhou et al. (2018) have proposed algorithms to retrieve snow depth from SMOS (Soil Moisture Ocean Satellite), with mitigated results. Recently, Braakmann-Folgmann and Donlon (2019) proposed to compute snow depth by combining AMSR2 and SMOS data with a neural network approach.

Studies have also estimated snow depth on sea ice in the Antarctic (e.g. Kacimi and Kwok, 2020; Maksym and Markus, 2008). AMSR-E and AMSR2 also provide snow depths in the Antarctic. Compare to the Arctic it covers the whole domain since all sea ice is considered as seasonal ice. The recent work of Kacimi and Kwok (2020) estimated the Antarctic snow depth from the difference between the total freeboard obtained from ICESat-2 lidar measurements and sea ice freeboards of CS-2. They show a very low variability in CS-2 freeboards compared to ICESat-2 and that ICESat-2 freeboards explained >90 % of the variance in snow depth. They also highlight that the validation sea ice parameters in the Antarctic remains a challenge since no seasonally or regionally diverse datasets from field records can be used to assess the large-scale satellite retrievals. Based on this, they urge that more sustained and extensive field measurements must be conducted. Note that such lidar minus Ku-band snow depths can also be computed in the Arctic. Finally, we note that a recent study of Fons et al. (2021), an update from the former approach experimented by Fons and Kurtz (2019), directly retrieves snow freeboard from CS-2 using a two-layer scattering model and shows promise in estimations of SIT in the Antarctic.

An alternative to passive microwave radiometry is to derive snow depth estimations from bi-frequency altimetric measurements. Recently, Armitage and Ridout (2015) have demonstrated this possibility by considering the difference of penetration between the Ku-band frequency CS-2 (13.5 GHz), assuming it is reflected near the snow–ice interface, and the Ka-band frequency SARAL/Altika (35.7 GHz), assuming it is reflected near the top of the snow pack, i.e the air–snow interface. Thereafter, a preliminary altimetric snow depth (ASD) version covering the 2013–2016 winter period has been developed at the LEGOS during the ESA CryoSeaNICE project (Guerreiro et al., 2016). Meanwhile, Lawrence et al. (2018) also developed a bi-frequency altimeter snow depth product (DuST) available for the Arctic. Similar to AMSR2B, the main difference between ASD and DuST is that DuST relies on a re-calibration of the Ka–Ku freeboards using OIB data, in this case to account for the difference of operating mode between CS-2 and SARAL (synthetic aperture radar (SAR) versus low-resolution mode (LRM) modes). These two products are the only existing publicly available snow depth products from altimetry.

A recent study of Zhou et al. (2021) presented an inter-comparison of available snow depth products from re-analyses, passive radiometry and altimetry (DuST). Similar, this paper reviews the state of the art by comparing current main snow depth estimations. Yet, the main objective is to present and assess the upgraded version of the ASD product (see Sect. 2), covering the 2013–2019 period in both hemispheres. The article fits in with the upcoming HPCM CRISTAL mission (Kern et al., 2020) and aims to demonstrate the potential of such snow depth data to further specific studies on, for instance, improved sea ice volume estimations, freshwater budgets, snow properties or data assimilation. Except from an analysis of the impact of the snow depth uncertainty on SIT retrieval, it does not explicitly address these questions.

The paper is organised as follows:

• First, we detail the methodology to process the ASD product and present all the datasets used in this study.

• Section 4 compares ASD with the other main existing snow depth satellite and model data in both hemispheres.

• The datasets are then assessed against OIB, CryoVEx and IMB validation data in the Arctic in Sect. 5.

• To circumvent the temporal limitation induced by SARAL, we propose in Sect. 6 a preliminary snow depth climatology based on ASD.

• Section 7 aims to quantify the SIT level of uncertainty due to the snow depth from an ensemble of SIT estimations calculated from the satellite and model snow depth datasets presented in the previous sections.

• We finally discuss and conclude, emphasising current needs for snow depth data in sea ice studies.

2 Data processing of ASD

The Altimetric Snow Depth (ASD) product presented here is an upgraded version of the data presented in Guerreiro et al. (2016). It has been developed at the LEGOS laboratory in the framework of the ESA CSAO+ (http://cryosat.mssl.ucl.ac.uk/csao/index.html, last access: March 2021) and Polar+ Snow on Ice projects. It has been freely available from the AVISO+/ODATIS national data centres since mid-2020 (http://ctoh.legos.obs-mip.fr/data/sea-ice-products, last access: March 2021). The ASD product includes data from March 2013 to October 2019 provided on a monthly basis for the 6 winter months (from November to April in the Arctic and from May to October in the Antarctic). It is projected onto a 500×500 EASE2 grid with a 12.5 km pixel size resolution. Snow depth calculation is based on the difference of penetration between Altika, the Ka-band range altimeter of SARAL (which is assumed to be reflected at the air–snow interface), and SIRAL, the Ku-band range altimeter of CS-2 (which is assumed to be reflected at the snow–ice interface). The main assumption is that the difference between these two surface elevations is only due to the penetration of the Ku radar in the snow pack and that the Ku radar penetrates fully to the snow–ice interface. Note that the validity of this hypothesis is strongly linked with the space and time variability in snow properties such as snow density, grain sizes, the liquid water content or the surface roughness. We use the CS-2 Ku-band waveforms in pseudo LRM (PLRM) extracted from geophysical ocean product (GOP) (https://earth.esa.int/documents/10174/125272/CryoSat-Baseline-C-Ocean-Product-Handbook, last access: March 2021) generated by the ESA CryoSat Ocean Processor (Bouffard et al., 2018a). The Ka-band SARAL/Altika waveforms are extracted from the Centre national d'études spatiales (CNES) SGDR-T (Sensor Geophysical Data Record) official products. Note that SARAL's orbit is limited to 81.5^∘ in latitude, which reduces the coverage of snow depth data in the Arctic. The use of a degraded version of the CS-2 SAR waveforms (the PLRM mode) allows us to have a footprint in accordance with the LRM mode of SARAL/Altika waveforms, which avoids having different impact of the surface roughness. Since the latest Baseline-C PLRM GOP was only available from 2017 at the time we computed the ASD data, we have used the Baseline-B for the period 2013–2016. It does not impact ASD since we use only the L1b product levels which have identical waveforms on both baselines. However, the Baseline-B does not include the SARIN (interferometric synthetic aperture radar) data. The next version of the ASD product will be produced with only the Baseline-C PLRM GOP to include all SARIN mode zones.

The first step to compute ASD data consists in extracting heights, H, from both Ku-band and Ka-band satellite waveforms following Eq. (1):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}H& =\mathrm{alt}-(\mathrm{range}-({\mathrm{tropo}}_{\mathrm{dry}}+{\mathrm{tropo}}_{\mathrm{wet}}+\mathrm{iono}\\ & +\mathrm{MSS}+{\mathrm{tide}}_{\mathrm{ocean}}+{\mathrm{tide}}_{\mathrm{bar}}\left)\right),\end{array}\end{array}$$

where alt is the satellite altitude, range is the altimetric range retracked from the waveforms, tropo_dry is the dry tropospheric correction, tropo_wet is the wet tropospheric correction, iono is the ionospheric correction, MSS is the DTU15 mean sea surface correction (Andersen and Knudsen, 2015), tide_ocean is the oceanic tide corrections, and tide_bar is the barometer tide correction.

The surface classification between heights corresponding to ocean surfaces (leads) and heights of sea ice surfaces (ice floes) is performed using the waveform pulse-peakiness (PP) criteria calculated from the waveforms (Eq. 2):

$$\begin{array}{}\text{(2)}& \mathrm{PP}={\displaystyle \frac{\max \left(\mathrm{WF}\right)}{{\sum }_{i=\mathrm{1}}^{N_{\mathrm{WF}}}{\mathrm{WF}}_i}},\end{array}$$

where WF is the discrete waveform echoes, and N_WF is the number of along-track measurements. Following Guerreiro et al. (2017), surfaces with a PP larger than 0.3 are considered as leads, and surfaces with a PP smaller than 0.1 are considered as sea ice floes. Observations with a PP between 0.1 and 0.3 are considered as mixed echoes and are discarded. This criteria is the same for SARAL and CS-2.

Altimetric ranges (range parameter in Eq. 1) of leads and floes are thereafter calculated using the TFMRA (Helm et al., 2014) with a 50 % threshold. A 25 km median filter is then applied to the 20 Hz along-track surface height of sea ice floes and sea ice leads. For each observation, the sea level anomalies (SLAs) under floes are calculated as the median value of all leads in a 25 km radius around the considered along-track point. Freeboard heights are computed by $\mathrm{fb}=H_{\mathrm{floe}}-H_{\mathrm{lead}}$, where H_floe is the heights of the ice floes and H_leads the heights of the leads. Finally, a 25 km radius median smoothing is applied to the retrieved freeboards. The reason is that we assume that the ASD product should not be able to provide relevant information at smaller scales. However, additional analyses and comparisons with validation data would be necessary to properly characterise this point. Note that the ability to consistently observe small scales would certainly be significantly improved in the future dual-frequency snow depth products from the CRISTAL mission. In the first ASD data version (Guerreiro et al., 2016), snow depths were only calculated at CS-2 and SARAL crossing track points. For this version, snow depths are calculated from the difference between monthly EASE2 gridded freeboards as ${\mathrm{sd}}_{\mathrm{r}}={\mathrm{fb}}_{\mathrm{Ka}}-{\mathrm{fb}}_{\mathrm{Ku}}$, where sd_r is a radar snow depth. Furthermore, to compute snow depth (sd), we need to take into account the decrease in the Ku radar echo velocity as it penetrates into the snow pack (Ulaby et al., 1986; Kwok and Cunningham, 2015):

$$\begin{array}{}\text{(3)}& \mathrm{sd}={\mathrm{sd}}_{\mathrm{r}}\times (\mathrm{1}+\mathrm{0.51}\times {\mathit{\rho }}_{\mathrm{s}}{)}^{(-\mathrm{1.5})},\end{array}$$

where ρ_s is the snow density, set to 300 kg/m³. Note this formulation agrees with the conventional interpretation to correct for the slower wave propagation speed through snow, as recently shown in Mallett et al. (2020).

To estimate snow depth uncertainties, we assume that errors are unbiased, uncorrelated and follow the Gaussian propagation law. Snow depth uncertainty is then given by

$$\begin{array}{}\text{(4)}& {\mathit{\delta }}_{\mathrm{sd}}=\sqrt{({\mathit{\delta }}_{{\mathrm{sd}}_{\mathrm{r}}}\times C{)}^{\mathrm{2}}+({\mathrm{sd}}_{\mathrm{r}}\times B\times {\mathit{\delta }}_{{\mathit{\rho }}_{\mathrm{s}}}{)}^{\mathrm{2}}},\end{array}$$

with

$$\begin{array}{rl}& C=(\mathrm{1}+\mathrm{0.51}\times {\mathit{\rho }}_{\mathrm{s}}{)}^{-\mathrm{1.5}},\\ & B=-\mathrm{1.5}\times \mathrm{0.51}\times (\mathrm{1}+\mathrm{0.51}\times {\mathit{\rho }}_{\mathrm{s}}{)}^{-\mathrm{2.5}},\\ & {\mathit{\delta }}_{{\mathrm{sd}}_{\mathrm{r}}}=\sqrt{{\mathit{\delta }}_{{\mathrm{fb}}_{\mathrm{Ka}}}^{\mathrm{2}}+{\mathit{\delta }}_{{\mathrm{fb}}_{\mathrm{Kur}}^{\mathrm{2}}}}\mathrm{and}{\mathit{\delta }}_{{\mathit{\rho }}_{\mathrm{s}}}=\mathrm{3.2}.\end{array}$$

The freeboard uncertainties ${\mathit{\delta }}_{\mathrm{fb}}=\sqrt{{\mathit{\delta }}_{H_{\mathrm{leads}}}^{\mathrm{2}}+{\mathit{\delta }}_{H_{\mathrm{floes}}}^{\mathrm{2}}}$ are calculated along-track from the statistical dispersion of heights in a 12.5 km radius. Since sea ice topography can significantly vary within this area, we assume that the statistical dispersion of floes is the same as that of the leads.

$$\begin{array}{}\text{(5)}& {\mathit{\delta }}_{H_{\mathrm{floes}}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{H_{\mathrm{leads}}}}{N_{\mathrm{floes}}}}\mathrm{and}{\mathit{\delta }}_{H_{\mathrm{leads}}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{H_{\mathrm{leads}}}}{N_{\mathrm{leads}}}}\end{array}$$

Using this methodology, Fig. 1 presents an example of ASD snow depth estimations and its associated uncertainties in the two hemispheres. In the Antarctic, thicker snow is located in the Weddell Sea and, to a lesser extent, in the coastal zones of the Bellingshausen and Amundsen seas. It is also relevant to identify the very low snow values associated with the Ross Ice Shelf. In the Arctic, the snow distribution follows well the dynamics of sea ice, the most characteristic of which is the export of MYI in the Beaufort Gyre. This figure also shows that the ASD data are very different from the W99m climatology, in which the W99m climatology tends to exhibit thicker snow layers over sea ice. The ASD mean level of uncertainties (of about 4 cm on average with standard deviation of 1 cm) are smaller than the deviations between these two products (on average about 6 cm with standard deviation of 7 cm).

.

Figure 1 Maps of the altimetric snow depth (ASD) annual mean snow depth (a) and its uncertainties (d) in the Antarctic and the Arctic (b and e) for 2015. The third column presents the annual mean snow depth map of the modified Warren-99 climatology (c) and its difference with the ASD product (ASD-W99m) in the Arctic (f) in 2015.

3 External datasets

This section presents the various snow depth datasets used in this study. For both hemispheres, the time period of model and satellite products ranges from 2014–2019 for the Arctic and 2013–2019 for the Antarctic (with limitations for some data products for both hemispheres), as explained in Fig. 2. Table A4 specifies the time period of the different satellite and model data.

3.1 Satellite data

3.1.1 DuST

The DuST (Lawrence et al., 2018) data are provided in the Arctic on a 1.5^∘ longitude × 0.5^∘ latitude grid by the CPOM/UCL (http://www.cpom.ucl.ac.uk/DuST/, last access: March 2021). They are available for the 6 winter month (November to April) until 2018. Similar to ASD data, DuST relies on the difference of penetration between the Ka-band SARAL/Altika and the Ku-band CS-2 radar altimeters. However, unlike the ASD product (which uses the CS-2 PLRM data; see Sect. 2), the DuST product uses SAR mode freeboard estimates of CS-2. The CS-2 and SARAL (LRM data) data are then calibrated against data from National Aeronautics and Space Administration's (NASA's) OIB airborne measurements to align the freeboard observations from SARAL/Altika with the snow–air interface and the freeboard observations from CS-2 with the snow–ice interface. Although this method provides a generic approach adaptable to different footprints (e.g. ICESat-1 and 2), it has the disadvantage of being dependent on the OIB data which are mainly located near the Canadian Archipelago and the Beaufort Sea and are only obtained during March and April.

3.1.2 AMSR2

The Advanced Microwave Scanning Radiometer-2 (AMSR2) is a passive radiometer on board the JAXA/GCOM‐W (Japan Aerospace Exploration Agency's Global Change Observation Mission – Water) satellite launched in 2012. Snow depth on sea ice retrieval is based on the gradient ratio of vertically polarised brightness temperature at 19 and 37 GHz (Markus and Cavalieri, 1998), and two snow depth products are available. The first product (AMSR2-NSIDC), described in Meier et al. (2018), is accessible on the NSIDC (https://nsidc.org/data/AU_SI12/versions/1, last access: March 2021). It provides daily L3 snow depth data in 12.5 km × 12.5 km stereopolar grids constructed as 5 d running averages. Data are only calculated over FYI (with MYI concentration of less than 20 %). This is due to the fact that MYI has a spectral signature similar to snow cover over FYI. Because of this limitation, we mainly focus on the Bremen AMSR2B product (Rostosky et al., 2018) developed at the Institute of Environmental Physics of the University of Bremen. This product is available on a daily basis from November 2012 to April 2018 on a polar stereographic grid with a 25 km × 25 km resolution. Unlike in AMSR2-NSIDC, March and April snow depth data are also calculated over MYI for AMSR2B. Other months are only available over FYI. Note that the algorithm used to calculate snow depth over MYI is also calibrated with observations from OIB.

3.2 Model data

3.2.1 PIOMAS

The pan-Arctic Ice-Ocean Modelling and Assimilation System (PIOMAS; Global Ice-Ocean Modeling and Assimilation System, GIOMAS, in the Antarctic) is a pan-Arctic coupled ocean and sea ice model developed for climate applications (Zhang and Rothrock, 2003). Snow depth data are provided on a daily basis in a 360 × 120 generalised curvilinear coordinate system from the PIOMAS model version 2.1 (http://psc.apl.uw.edu/research/projects/arctic-sea-ice-volume-anomaly/data/model_grid, last access: March 2021). The atmospheric surface forcing fields are issued from NCEP–NCAR reanalysis. Ocean and Sea Ice Satellite Application Facility (OSI SAF) sea ice concentration from EUMETSAT is assimilated in near real time. A complete description can be found in Schweiger et al. (2011).

3.2.2 MERCATOR

In the framework of the Copernicus Marine Environment Monitoring Service (CMEMS), Mercator Ocean implements a real-time global analysis and forecasting system at $\mathrm{1}/\mathrm{12}$^∘ resolution. This system is based on the NEMO ocean platform (Madec et al., 2015) and is driven by atmospheric conditions from the ECMWF Integrated Forecast System (IFS) analysis and forecasting system. The sea ice model component is the Louvain-la-Neuve sea ice model (LIM-2) (Fichefet and Maqueda, 1997; Vancoppenolle et al., 2012). The unique sea ice quantity assimilated within the system is the near-real-time sea ice concentration from OSI SAF provided by CMEMS. The accumulation of the snow depth onto sea ice originates from ECMWF snowfall forcing fields. A complete description can be found in Lellouche et al. (2018). Snow depth data used in this study have been provided on a monthly basis on the ORCA tripolar grid. This LIM-2 model configuration will be referred to as the MERCATOR model.

3.2.3 NESOSIM

The NASA Eulerian Snow On Sea Ice Model (NESOSIM) is a 3D snow budget model configured to produce snow depth and density over sea ice in the Arctic Ocean. Data used in this paper are snow depth monthly mean maps provided on a 100 km × 100 km stereographic polar grid issued from the NESOSIM 1.0 configuration (Petty et al., 2018). It is freely available at https://earth.gsfc.nasa.gov/cryo/data/nasa-eulerian-snow-sea-ice-model-nesosim (last access: March 2021). Snowfall data forcing fields are the ECMWF ERA-Interim global reanalysis and a median of three reanalyses: the 55-year reanalysis of the JRA, NASA's MERRA and ASR, version 1.

Figure 2 Times series of annual mean snow depth of the different products in both hemispheres. Annual means are calculated as the average of the spatial means of all monthly grid maps. Only the 6 winter months are considered in the calculation (November–April in the Arctic and May–October in the Antarctic). Note that the AMSR2 data are not represented in the Arctic since the AMSR2-NSIDC product is only available over FYI, and the AMSR2B product is only available with its full spatial coverage (FYI and MYI) in March and April.

3.3 Validation data

The validation data presented in this section are only compared with ASD in the Arctic. As mentioned in the introduction CryoVEx and OIB data in the Antarctic are still only L1b. IMB buoys are only in the Arctic. The geographical locations of the validation data are shown in Sect. 4.3.

3.3.1 CryoSat Validation EXperiment (CryoVEx)

These airborne data are a joint effort of the DTU National Space Institute and ESA (in cooperation with AWI). The main goal of the project is to quantify and validate CS-2's ability to measure SIT. Since 2003, the CryoVEx missions have provided both laser (ALS, airborne laser scanner) and Ku-band radar (ASIRAS, ESA airborne Ku‐band interferometric radar) altimetric measurements. In addition, in 2017, the data include for the first time along-track total freeboards from the Ka-band KAREN altimeter (Haas et al., 2017), allowing us to compute Ka–Ku snow depth estimations. These data have been developed by AWI in the context of the ESA CryoSeaNICE project. For further information, please see the technical report at https://earth.esa.int/eogateway/documents/20142/1526226/CryoVEx2017-final-report.pdf (last access: March 2021).

3.3.2 Operation Ice Bridge (OIB)

OIB is one of the largest airborne missions in polar regions aiming to determine sea ice properties. Among others, it carries a snow radar measuring the snow–air and snow–ice interfaces of the signal scattered by the area illuminated beneath the aircraft. All information concerning the different campaigns and instruments can be found in the literature (e.g. Newman et al., 2014; Kurtz et al., 2013; Armitage and Ridout, 2015). OIB snow depth data presented in this paper are the NSIDC OIB Quicklook version (Kurtz et al., 2012; King et al., 2015) validated from in situ data of the BROMEX campaign (Nghiem et al., 2013). These data are available at https://daacdata.apps.nsidc.org/pub/DATASETS/ICEBRIDGE/Evaluation_Products/ (last access: March 2021). For comparison with the ASD dataset, we use all spring campaigns during the period of 2014–2018. For comparison with the ASD climatology, we use all OIB spring campaigns during 2009–2012.

3.3.3 Ice mass balance (IMB) buoys

IMB data (Perovich et al., 2021) are issued from the US Army CRREL-Dartmouth Mass Balance Buoy Program. They provide point measurements of the ice mass balance characterised by a combination of snow and ice conditions (Richter-Menge et al., 2006; Perovich and Richter-Menge, 2015). Snow depth data are collected at 4 h intervals from acoustic sounders. In this study we analyse the data of three winters (2013–2014, 2014–2015 and 2015–2016). Note that we do not use the data of 2017 and 2018 because they are still provisional and subject to revision.

4 Comparison between snow depth data

4.1 Methodology

Model and satellite snow depth estimations are projected onto a 12.5 km pixel size EASE2 grid (similar to ASD) using a linear 2D multivariate interpolation. For each product we compute monthly mean maps for the 6 winter months (November–April in the Arctic and May–October in the Antarctic) for both hemispheres from March 2013 to October 2019. These maps are provided for 2015 in the Supplement (Figs. S1 to S13). Note that to improve the consistency of the analyses we always use the larger common time period for all the compared products. Thus, the covered time periods may vary depending on the availability of the snow depth products. Three statistical diagnoses are used to characterise the differences between snow depth datasets:

1. The climatic annual mean is simply the average of all monthly snow depth maps from all years.

2. The mean annual variability (MAV) is the average, for the years y, of the annual snow depth standard deviation (SD_y(sd)) maps computed from the six monthly mean maps. The MAV is calculated at each grid point following Eq. (6):

   $$\begin{array}{}\text{(6)}& \begin{array}{rl}& \mathrm{MAV}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{y}}}}\sum _{\mathrm{y}}\left({\mathrm{SD}}_{\mathrm{y}}\right(\mathrm{sd}\left)\right)\text{with}\\ & {\mathrm{SD}}_{\mathrm{y}}\left(\mathrm{sd}\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{m}}}}\sum _{\mathrm{m}}({\mathrm{sd}}_{\mathrm{y},\mathrm{m}}-\stackrel{\mathrm{‾}}{{\mathrm{sd}}_{\mathrm{y}}}{)}^{\mathrm{2}}},\end{array}\end{array}$$

   where sd_y, m is the monthly snow depth map of the month m ($\mathrm{m}=\mathrm{11},\mathrm{12},\mathrm{01},\mathrm{02},\mathrm{03},\mathrm{04}$ in the Arctic and $\mathrm{m}=\mathrm{05},\mathrm{06},\mathrm{07},\mathrm{08},\mathrm{09},\mathrm{10}$ in the Antarctic) for the year y, and $\stackrel{\mathrm{‾}}{{\mathrm{sd}}_{\mathrm{y}}}$ is the snow depth annual mean map of the year y. N_m and N_y are the number of months (six in our case) and the number of year in the considered time period, respectively.

3. The mean interannual variability (MIV), which is the average of the interannual standard deviation SD_m(sd) maps for the 6 winter months m over the 2013–2019 time period. The MIV is calculated at each grid point following Eq. (7):

   $$\begin{array}{}\text{(7)}& \begin{array}{rl}& \mathrm{MIV}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{m}}}}\sum _{\mathrm{y}}\left({\mathrm{SD}}_{\mathrm{m}}\right(\mathrm{sd}\left)\right)\text{with}\\ & {\mathrm{SD}}_{\mathrm{m}}\left(\mathrm{sd}\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{y}}}}\sum _{\mathrm{y}}({\mathrm{sd}}_{\mathrm{y},\mathrm{m}}-\stackrel{\mathrm{‾}}{{\mathrm{sd}}_{\mathrm{m}}}{)}^{\mathrm{2}}},\end{array}\end{array}$$

   where $\stackrel{\mathrm{‾}}{{\mathrm{sd}}_{\mathrm{m}}}$ is average snow depth map of the month m.

For now, the lack of available relevant snow depth data in the Antarctic limits these comparisons to the Northern Hemisphere. Comparisons with all in situ and validation data are performed with a comparable methodology: snow depth model and satellite gridded maps are projected along the aircraft trajectories using a bilinear interpolation. Airborne and in situ are generally daily data, and the comparisons are performed with the mean maps of the month to which that day belongs. In order to achieve similar spatial scales in the comparisons, we have applied a 25 km window rolling mean to smooth the external data. Due to the non-Gaussianities within 25 km sections, we have noticed that the smoothing has a mean effect which slightly tends to reduce the mean value of data. It is important to consider that this allows us to correct potential bias that would be induced by a difference in resolution. Note that the space and time consistency in resolution between all model and satellite data is also ensured by projection onto similar grids. Since snow depth products and validation data are available over different time periods (see table A4), the comparisons are also performed for different time periods. Our approach is to take, for each comparison, the time period common to all the compared data. The aim is to provide reliable statistics. For each comparison, a label specifies the considered time period.

4.2 Comparison between satellite and model estimations

4.2.1 Results in the Arctic

Here, we compare the ASD, DuST and AMSR2B snow depth data in the Arctic by considering the months of March and April over the 2014–2018 period (Fig. 3). These climatic “seasonal” means and interannual variability presented in Fig. 3 are the mean and the standard deviation of all monthly maps of March and April over this time period. Statistical diagnoses are summarised in Table A1.

Figure 3 Comparisons of satellite snow depth estimations in the Arctic. (a, b, c) Bi-monthly March and April “climatological” snow depth mean maps calculated over the period 2014–2018 for the ASD data (a, d), the DuST data (b, e) and the AMSR2B data (c, f). (d, e, f) The standard deviation between all the data used to calculate the mean snow depth maps (equivalent to the mean interannual variability computed with March and April over the same 2014–2018 time period).

On average, the ASD data are about 3 cm lower than the AMSR2B data and 6 cm lower than the DuST data. ASD data present the stronger snow depth gradients (associated with a higher spatial standard deviation) with clear distinctions between low snow depths over thin sea ice regions (e.g. Laptev Sea region) and higher snow depths over regions of thicker sea ice (e.g. Queen Elizabeth Islands region). Except in the MYI zones, AMSR2B data are quite smooth and characterised by a weak variability compared to that of the two altimetric products. In MYI zones (towards Queen Elizabeth islands) AMSR2B has the highest snow depths. We assume that it might be related to the snow depth retrieval algorithm, inadequate over thick ice. Indeed, the statistical results which are presented in the Appendix (Table A1) show that the two Ka–Ku products (ASD and DuST) are in good agreement in terms of spatial variability and the magnitude of their annual cycle considering the entire set of data over the 6 winter months. Note that DuST presents the highest climatic “seasonal” mean solution, even higher than W99m (see also Fig. 2). In addition, DuST and ASD spatial distributions are comparable, showcasing that the difference between the two Ka–Ku altimetric products is likely a bias resulting from the re-calibration of the DuST data with OIB (this point will also be discussed in the next sections). This underlines the consistency of this methodology used to compute bi-frequency snow depth estimations in the absence of external data (ASD) since DuST calibration is limited by the spatial and temporal availability of OIB.

Figure 4 Comparisons between the ASD data and model snow depth products (PIOMAS, NESOSIM, MERCATOR) in the Arctic. The first column presents the “climatic mean” snow depth maps, the second column the “annual variability” and the third column the “interannual variability” as presented in Sect. 4.1. These diagnoses are calculated from March 2013 to April 2019 for ASD (first row), PIOMAS (second row) and MERCATOR (fourth row) and from March 2013 to April 2017 for NESOSIM (third row).

Figure 4 compares the ASD snow depth estimations with the outputs of three sea ice models in the Arctic (NESOSIM, MERCATOR and PIOMAS). The statistical diagnoses presented in Sect. 4.1 are used. They are summarised in Table A1. Amongst the three models, the MERCATOR model provides lower snow depth estimations throughout with a climatic mean of 7.4 cm and lower mean interannual and annual variability of about 3 cm. On the contrary, climatic means of PIOMAS and NESOSIM are similar both in terms of spatial distribution and pan-Arctic mean (with values of respectively 15.5 and 16.2 cm). Considering annual and interannual variability, PIOMAS presents high values that are more evenly spread. In contrast, NESOSIM is characterised by lower values everywhere except in the Greenland and the Kara seas.

In spite of some comparable large-scale spatial patterns, these results highlight the large deviations between the ASD data and the model solutions. One striking feature is the transport of MYI along the Canadian coast by the Beaufort Gyre (see also Fig. 1), which is nearly non-existent in model solutions. The patterns of deep snow over MYI are well represented in the models but do not extend along the Canadian coastline towards the Beaufort Sea, as is the case for ASD. Overall, the models tend to overestimate snow depths nearly everywhere compared to the ASD data. The ASD data climatic mean is lower (12.3 cm), and zones of thin snow layers (e.g. around 120^∘ E) are less pronounced in models (except by MERCATOR, where snow depth is low everywhere). Models also present higher variability. However, in all datasets the maximum variability occurs between 0 and 90^∘ E, which is relevant because of the proximity to the open ocean and the existence of strong meteorological events (e.g. Semenov et al., 2019; Dong et al., 2019). Although investigating what causes these discrepancies is beyond the scope of this article, the strong sensitivity of the models to the reanalysis snowfall forcing data (e.g. Boisvert et al., 2018; Petty et al., 2018) could play a predominant role here.

4.2.2 Results in the Antarctic

ASD is currently the only publicly available altimetric snow depth product in the Antarctic. One major advantage compared to the Arctic is that the SARAL orbit does not affect the data coverage since all the Antarctic sea ice is below 81.5^∘ S. Similarly to what has been done in the previous section for the Arctic, Fig. 5 compares ASD with the AMSR2-NSIDC data and the GIOMAS and the MERCATOR sea ice model. The statistics are summarised in Table A2.

In spite of a comparable climatic mean and some coherent spatial distributions such as the transport of snow in the Weddell Gyre, ASD and AMSR2-NSIDC data show significant discrepancies. Unlike what we have observed in the Arctic, ASD data are relatively smooth with a weak spatial variability. In comparison, AMSR2-NSIDC data are characterised by strong contrasts between large snow depth patterns in the Weddell Sea and to a lesser extent in the western part (from 150^∘ E to 90^∘ W) and areas with nearly no snow (eastern part, between 0 and 150^∘ E). As in the Arctic, AMSR2-NSIDC annual and interannual variability is lower and much more localised than in the ASD data. ASD and AMSR2-NSIDC reproduce well the thin snow pattern of the Ross Ice Shelf. One other relevant difference between these two datasets is the systematic decrease in snow depth in October in the AMSR2-NSIDC data (see Table A2).

The presence of thicker snow in East Antarctica for ASD is consistent with Worby et al. (2008b) in that comparisons with ARISE in situ data have shown radiometric measurement snow depth underestimations. Wet-snow conditions due to flooding might be responsible for radiometric brightness temperature contrasts. The low variabilities in the Weddell Sea for ASD were not expected since winter snow properties are extremely variable (Massom et al., 1997). In addition the large bias with AMSR2-NSIDC raises questions. Since the Ku band is supposed to better penetrate in cold and dry snow (Willatt et al., 2011), one hypothesis is that the presence of saline and warm moist basal snow layers, even in the absence of flooding (Massom et al., 2001; Perovich and Richter-Menge, 1994), lead to ASD underestimations. For AMSR2-NSIDC, the snow depth retrieval algorithm is probably not well adapted for rougher snow that can be compared to MYI in the Arctic.

With a climatic spatial mean of about 34.6 cm, GIOMAS simulates everywhere very high snow depths compared to the other solutions. Considering what is provided in the literature (e.g. Worby et al., 2011), these values clearly seem to be overestimated. With the exception of the Weddell Sea, the MERCATOR model and ASD present comparable snow patterns with substantial coastal values not seen in the AMSR2-NSIDC data. However, the MERCATOR model displays thicker snow everywhere but in the Weddell Sea, while the presence of deep snow is a well observed feature, which is supported by ASD (e.g. Massom et al., 1997; Eicken et al., 1994). Because of highly dynamic weather conditions, with persistent strong winds in the Antarctic, snow thickness distribution is not directly related to snowfalls (Massom et al., 2001). Conversion from precipitation to snow depth is then very different to the Arctic, and Antarctic snow pack is not a uniform slab. These features should partly explain model difficulties and differences between the two hemispheres.

Figure 5 Comparisons between ASD data and other snow depth estimations in the Antarctic (AMSR2-NSIDC, MERCATOR, GIOMAS). The first column presents the “climatic mean” snow depth maps, the second column the “annual variability” and the third column the “interannual variability” as presented in Sect. 4.1. These diagnoses are calculated from May 2013 to October 2019 for the ASD data (first row), the AMSR2-NSIDC (second row) and the MERCATOR model (third row) and from May 2013 to October 2016 for the GIOMAS model data (fourth row).

Figure 6 Along-track comparison between the snow depth products and the OIB snow radar data acquired on 29 March 2015 in the Arctic. The red envelope refers to the ASD uncertainties. The geographical map specifies the geographical location of this particular 1 d OIB track.

4.3 Comparison with in situ and airborne data

4.3.1 OIB

OIB data are compared to all the datasets presented in the previous sections. We only consider OIB observations below 81.5^∘ N from the 2014–2018 campaigns. Figure 6 shows an along-track example from the OIB flight on 29 March 2015, and Fig. 7 presents scatter plots comparing the various snow depth products with the OIB data over the entire 2014–2018 time period (including all spring OIB campaigns within this period). Statistical results are summarised in Table 1.

One striking feature is the good consistency between ASD and OIB in terms of magnitude and spatial variability. The OIB data are almost within the ASD envelope of uncertainties (in shaded red) at all times. These consistencies between OIB and ASD data are also demonstrated in several other OIB tracks provided in the Supplement (Figs. S14 to S17). The DuST Ka–Ku data present an along-track variability similar to ASD but tend to overestimate snow depth compared to OIB. This reinforces the hypothesis that the main difference between ASD and DuST is a bias. The AMSR2B data also reproduce the large-scale snow depth variability quite well but overestimate at a level comparable with DuST when compared to OIB. The DuST and AMSR2B overestimations are quite consistent with Kwok et al. (2017), who found that this OIB product tends to have thinner snow than the NOAA Wavelet Airborne Snow Radar dataset, which has been used for re-calibrations. In fact, it is important to be aware that various OIB datasets have been produced, leading to a variety of snow estimations that can reach 7 cm on FYI and 12 cm on MYI (Kwok et al., 2017). The variabilities in snow properties from OIB daily basis data to the monthly means of the other datasets certainly explain some of the deviations.

Models hardly represent the along-track variability and do not reproduce small-scale patterns. A resolution that is too coarse could partly explain this feature. As expected, the MERCATOR model snow thicknesses are far below the others, whereas the PIOMAS and NESOSIM solutions provide comparable results (see also Table 1). As expected, the W99m climatology strongly overestimates snow depths and is clearly a less optimal solution in this case.

Figure 7 Scatter plots between the snow depth products and OIB snow radar data considering the 2014 to 2018 missions. The red line is the linear fitting line. For NESOSIM, 2018 is not included.

Table 1 Statistics of the comparison between OIB snow radar data and the various snow depth products considering the OIB spring campaigns of 2014–2018.

4.3.2 IMB

The comparison between snow depth estimations and the IMB data for the winters of 2013–2014, 2014–2015 and 2015–2016 is presented in Fig. 8. Statistical results are summarised in Table 2. Except for AMSR2B, which is only available for the months of March and April, all comparisons are performed using the exact same spatial coverage considering only data below 81.5^∘ N. Note that only very limited IMB data are available south of 81.5^∘ N for the winter of 2015–2016.

Figure 8 and Table 2 show how the satellite-derived mean values are closer to the in situ observations than the model. However, the seasonal changes captured in the IMBs are not reproduced well by the satellite products or the models. This could reflect sampling differences as the IMB observations consist of hourly point data that are not representative of the wider averages of the gridded satellite or model products. DuST and ASD along-track variability is once again very similar, with DuST being overall higher, highlighting the very likely bias between these two products.

Figure 8 Along-track comparison between the snow depth products and the IMB in situ data. (a) IMB buoys compared with satellite data, both passive and active, and (b) IMB buoys compared with model data. The right column shows geographical maps from data for all three periods (2013–2014, 2014–2015, 2015–2016). The maps present the geographical location of the three IMB buoys.

Table 2 Statistics of the comparison with IMB data and the various snow depth products considering data acquired during 2013–2016.

4.3.3 CryoVEx

Figure 9 presents the snow depth products projected along the CryoVEx flight track measured on 31 March 2017. This figure gathers the CryoVEx snow depth estimations calculated from (a) the difference between the Ku-band ASIRAS radar freeboard and the ALS laser (ALS/ASR) total freeboard and from (b) the Ku-band ASIRAS radar freeboard and the Karen (Ka/ASR) total freeboard. Note that because of SARAL's spatial coverage limitation (latitude < 81.5^∘) and the absence of the SARIN mode in the CS-2 ESA Baseline-B GOP (used to construct the ASD data; see Sect. 2), only the track of the 31 March 2017 CryoVEx survey of Baffin Bay can be compared with the ASD data.

Figure 9 Along-track comparison, in the Arctic, between the snow depth products (models and satellites) and the CryoVEx snow depth data calculated from the Ka-band KAREN radar (blue) and the ALS laser (dark blue) with the Ku-band ASIRAS radar during the 31 March 2017 campaign. The map specifies the geographical location of the CryoVEx track.

Ka/ASR estimation exhibits very thin snow thickness (<10 cm). Although it might be expected in this area (Landy et al., 2017), this solution still contains unrealistic negative values due to the fact that ASIRAS and KAREN freeboards are nearly equivalent over FYI (without negative values). The main difference between snow depths calculated from the laser and the Ka-band-radar snow depth estimations seems to be a bias. Except around 71.3^∘ N, where a specific event may have occurred, the ASD snow product tracks the magnitude of the CryoVEx airborne data, with ALS/ASR and Ka/ASR acting respectively as an upper and lower bound and being within ASD's range of uncertainty (in red shading). However, the along-track variability is not well reproduced by any dataset. The MERCATOR and NESOSIM models also show similar magnitudes as the CryoVEx snow depth estimations (ALS/ASR and Ka/ASR), while PIOMAS and the W99m climatology clearly overestimate. AMSR2B and DuST overestimate snow depth, while the DuST along-track variability is comparable to ASD.

5 Towards an ASD snow depth climatology

In spite of the consistency of the ASD data highlighted in the previous sections from comparisons with airborne and in situ data, the limited temporal coverage (only available post-2013) remains an important limitation to the use of these data. Considering the low interannual snow precipitation variability compared to that of the spatial variability (e.g. Figs. 2 and 4), one solution would be to develop an ASD climatology. For that purpose, we have constructed a preliminary altimetric snow depth climatology in the Arctic by averaging all the ASD snow depth maps of each month during the 2013–2019 period (designated as ASD-clim). An equivalent climatology could also be constructed for the Antarctic, but the lack of validation prevents its validity to be demonstrated at present. To demonstrate the relevance of this climatology in the Arctic for the years prior to 2013, the ASD-clim data are projected on all the tracks of the four OIB missions occurring between 2009 and 2012 and are presented in Fig. 10.

Figure 10 Scatterplots comparing the ASD and the W99m climatologies with OIB by considering all the tracks of the 2009, 2010, 2011 and 2012 OIB missions.

Figure 10 shows that the ASD climatology would be a better solution than W99m, at least for the CS-2 time period. This is highlighted with the lower bias (−0.1 m compared with −0.14 m), the higher correlation coefficients (0.71 compared with 0.06) and the lower RMSE (0.06 m compared with 0.17 m). Furthermore, even as a climatology, we obtained a correlation between the ASD-clim and OIB of the same order of magnitude as obtained by direct correlation between ASD and OIB measurements (Fig. 7). However, as future work, a more refined computation should be investigated. In particular, an important point would be to provide a relevant snow product representing the Envisat era. Since the W99m can be considered as relatively consistent until the year 2000, a mixed ASD and W99m product (ASD-W99m), including, for instance, a time dependency, could be considered. Also, it could be relevant to take into account the sea ice type interannual variability as is done for W99m.

Another critical point of the ASD data (for instance for sea ice volume estimations) is the reduced spatial coverage due to the SARAL orbit. While during the time period of Envisat, this does not become an issue since its orbit is equivalent to that of SARAL; the absence of data north of 81.5^∘ N is a major limitation to SIT estimations. A combination of ASD with W99m and/or model data could be done to extrapolate the data. Potentially, one could also investigate a combination of ASD with W99m and/or models to extrapolate the data. In the Antarctic, we do not encounter the same issue since the orbit does not reduce the spatial coverage, and FYI is generally considered to be the only sea ice type present in the Southern Ocean. Then, such a climatology would be a solution to compute SIT estimations over the Envisat and CS-2 time period provided that the temporal variability in precipitation is investigated. Note that experimental SIT estimations (Hendricks et al., 2018) have already been done in the Antarctic in the framework of the Sea Ice Climate Change Initiative (SI-CCI). These data used a snow climatology elaborated from AMSR-E and AMSR2 data in a nearly equivalent manner to what has been presented here.

6 Impacts of snow depth on SIT estimation

Since there is still no validated freeboard product in the Antarctic, this analysis has only been performed in the Arctic. From the CS-2 radar freeboard product presented in Guerreiro et al. (2017), extended to 2019, we have computed SIT estimates using the various snow depth datasets presented in this study. They are calculated from the CS-2 L1b ESA Baseline-C SAR waveforms (Bouffard et al., 2018b) with the methodology presented in Sect. 2. Freeboard heights are converted into SIT assuming hydrostatic equilibrium between snow-covered sea ice and the ocean (Eq. 8).

$$\begin{array}{}\text{(8)}& \mathrm{SIT}={\displaystyle \frac{{\mathit{\rho }}_{{\mathrm{wfb}}_{\mathrm{ice}}}+{\mathit{\rho }}_s\mathrm{sd}}{{\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_i}}\end{array}$$

As in Ricker et al. (2014), we assume the seawater density (ρ_w) is set to 1024 kg/m³. Consistent with the approach of Laxon et al. (2013), sea ice density (ρ_i) is set to 882 kg/m³ for MYI and 917 kg/m³ for FYI (Alexandrov et al., 2010). Furthermore, the snow density (ρ_s) is set to 300 kg/m³. A distinction between FYI and MYI follows the EUMETSAT OSI SAF sea ice type classification. The sea ice freeboard fb_ice is calculated from the radar freeboards, taking into account the reduction in radar velocity within the snow layer, which depends on snow depth data following Eq. (3). SIT monthly maps are calculated for each snow depth product presented in previous sections. These monthly maps are provided for the winter of 2015 in the Supplement (Figs. S18 to S24). Figure 11 presents the annual mean time series of the spatial mean of all these SIT solutions (except for AMSR2B which is only calculated in March and April).

Figure 11 Arctic annual mean time series (2014–2019) of the averaged SIT solutions computed using various snow depth products (models and active/passive satellite-based). Spatial coverage is defined by the smallest sea ice extent. The dotted blue lines define the annual mean, minimum and maximum values between all solutions (except for MERCATOR). The AMSR2 data are not represented since they are not available over the period.

Except for the MERCATOR model (which has shown to be overall significantly low), using the ASD data provide the lowest pan-Arctic mean SIT estimations. The highest SIT estimations are obtained with the DuST product (and W99m), highlighting the impact of the bias previously identified between these two solutions (ASD and DuST). Considering the various snow products together (except for the MERCATOR model), the spatial mean SIT annual mean values are distributed between 1.2 and 1.6 m. This leads to a pan-Arctic SIT ensemble mean of 1.4 m (mean of the dotted blue line in Fig. 11) associated with a mean maximum deviation of 29 cm (mean of differences between the maximum and minimum dotted blue lines). This mean maximum deviation, which is approximately the mean difference between ASD and DuST, represents about 20 % of the mean SIT. Considering the good consistency between ASD and the validation snow depth data (demonstrated in Sect. 4.3), actual SIT products computed using W99m should tend to significantly overestimate sea ice volume in the Arctic. For instance, taking a deviation of about 0.3 m (≈ the mean deviation between ASD and W99m in 2015) and a mean Arctic sea ice extent of 13×10⁶ km² (approximately the mean sea ice extent of 2015), the resulting SIT overestimation would lead to a sea ice volume excess of about 4×10³ km³, which represents about 3.5×10¹⁵ L of freshwater. Such an uncertainty is not negligible compared to the 8000 ± 2000 km³ change in the western Arctic freshwater content between 1995–1996 and 2009–2010 relayed in Giles et al. (2012). To give an illustrating example, this uncertainty of 3.5×10¹⁵ L of freshwater would be roughly equivalent to 6 months of the Amazon water discharge (considering a debit of about 2×10⁸ L/s). Deviations of only a few centimetres of snow depth (≈ 6 cm of mean deviation between ASD and W99m snow depths in 2015; see Fig. 1) have strong impacts on the freshwater budget.

Although this provides a relevant overview, such large-scale pan-Arctic means do not show the spatial and temporal variability in the various SIT solutions. For instance, while DuST and W99m provide comparable SIT pan-Arctic mean estimations, the spatial and temporal information included in these two datasets are very different (e.g. Figs. 2 and 3). To better characterise the impacts of snow depth on SIT, we propose to consider the deviations between these SIT estimates for each month of each year and at each grid point. Such an ensemble of SIT solutions will allow us to characterise the level of uncertainty due to discrepancies between the different snow depths. For the month m of year y, the standard deviation (SD_y, m(SIT)) and the maximum deviation (maxdev_y, m(SIT)) maps are calculated with Eqs. (9) and (10). ${\mathrm{SIT}}_{\mathrm{y},\mathrm{m},\mathrm{p}}$ is the SIT solution provided by the snow product p, and $\stackrel{\mathrm{‾}}{{\mathrm{SIT}}_{\mathrm{y},\mathrm{m}}}$ is the mean map between all these SIT estimations. N_p is the number of snow depth products, and ${max}_{\mathrm{p}}\left({\mathrm{SIT}}_{\mathrm{y},\mathrm{m},\mathrm{p}}\right)$ and ${min}_{\mathrm{p}}\left({\mathrm{SIT}}_{\mathrm{y},\mathrm{m},\mathrm{p}}\right)$ are respectively the maximum and minimum monthly maps between the SIT solution of the snow depth products.

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathrm{SD}}_{\mathrm{y},\mathrm{m}}\left(\mathrm{SIT}\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{p}}}}\sum _{\mathrm{p}}({\mathrm{SIT}}_{\mathrm{y},\mathrm{m},\mathrm{p}}-\stackrel{\mathrm{‾}}{{\mathrm{SIT}}_{\mathrm{y},\mathrm{m}}}{)}^{\mathrm{2}}}\text{(10)}& {\displaystyle }{\mathrm{maxdev}}_{\mathrm{y},\mathrm{m}}\left(\mathrm{SIT}\right)=\underset{\mathrm{p}}{max}\left({\mathrm{SIT}}_{\mathrm{y},\mathrm{m},\mathrm{p}}\right)-\underset{\mathrm{p}}{min}\left({\mathrm{SIT}}_{\mathrm{y},\mathrm{m},\mathrm{p}}\right)\end{array}$$

For each winter month, interannual mean, minimum and maximum maps of these two statistics (standard deviation SD_y, m(SIT) and the maximum deviation maxdev_y, m(SIT)) can then be computed from Eqs. (11) to (16), where N_y is the number of years.

$$\begin{array}{}\text{(11)}& {\displaystyle }\stackrel{\mathrm{‾}}{{\mathrm{SD}}_{\mathrm{m}}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{y}}}}\sum _{\mathrm{y}}\left({\mathrm{SD}}_{\mathrm{y},\mathrm{m}}\right(\mathrm{SIT}\left)\right)\text{(12)}& {\displaystyle }\stackrel{\mathrm{‾}}{{\mathrm{maxdev}}_{\mathrm{m}}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{y}}}}\sum _{\mathrm{y}}\left({\mathrm{maxdev}}_{\mathrm{y},\mathrm{m}}\right(\mathrm{SIT}\left)\right)\text{(13)}& {\displaystyle }\max \left({\mathrm{SD}}_{\mathrm{m}}\right)=\underset{\mathrm{y}}{max}\left({\mathrm{SD}}_{\mathrm{y},\mathrm{m}}\right(\mathrm{SIT}\left)\right)\text{(14)}& {\displaystyle }\max \left({\mathrm{maxdev}}_{\mathrm{m}}\right)=\underset{\mathrm{y}}{max}\left({\mathrm{maxdev}}_{\mathrm{y},\mathrm{m}}\right(\mathrm{SIT}\left)\right)\text{(15)}& {\displaystyle }\min \left({\mathrm{SD}}_{\mathrm{m}}\right)=\underset{\mathrm{y}}{min}\left({\mathrm{SD}}_{\mathrm{y},\mathrm{m}}\right(\mathrm{SIT}\left)\right)\text{(16)}& {\displaystyle }\min \left({\mathrm{maxdev}}_{\mathrm{m}}\right)=\underset{\mathrm{y}}{min}\left({\mathrm{maxdev}}_{\mathrm{y},\mathrm{m}}\right(\mathrm{SIT}\left)\right)\end{array}$$

These statistics quantify the SIT uncertainty, caused only for snow depth, assuming that deviations between snow depth solutions refer to the snow depth level of uncertainty. We define two cases: (1) “obs snow products”, which considers ASD, DuST, W99m and AMSR2B, and (2) “all snow products”, which includes in addition the PIOMAS, MERCATOR and NESOSIM model solutions. Note that ASD-clim is not used in calculations in order to avoid the redundancy with ASD data. Considering “obs snow products”, Fig. 12 presents maps of the metrics presented in Eqs. (11) to (16) for April. In addition, the monthly minimum, mean and maximum maps over the period 2014–2019 are presented for April in the Supplement (Fig. S25). All the results are summarised in Table A3 which presents, for each winter month, the spatial mean value of these maps. The spatial mean of the monthly mean SIT maps ($\stackrel{\mathrm{‾}}{{\mathrm{SIT}}_{\mathrm{y},\mathrm{m},\mathrm{p}}}$) are also included.

Figure 12 Interannual minimum (min(SD_m), first column), mean ($\stackrel{\mathrm{‾}}{{\mathrm{SD}}_{\mathrm{m}}}$, second column) and maximum (max(SD_m), third column) maps of the standard deviation (first row) and the maximum deviation (second row) for the month of April (m=04). The snow products used to compute the maps are ASD, DuST, AMSR2B and W99m. These maps correspond to the case “obs snow products”.

The highest impacts on SIT are generally located near coastal zones and marginal ice zones. In these areas, maximum standard deviation can reach up to 50 cm and maximum deviation to 80 cm when considering “obs snow products”. Due to larger differences between model and observations, in particular MERCATOR, these values are higher with “all snow products” (see Table A3). Except for these specific locations that mostly account for maximum values, the standard deviation maps are rather uniform. This suggests that spatial mean values (presented in Table A3) are a relevant indicator.

The large deviations between minimum and maximum (standard deviation and/or maximum deviation) maps indicate significant variations in the SIT level uncertainty due to snow depth variability from one year to another. On the contrary, spatial mean standard deviations do not vary significantly in between months. Considering “obs snow products” and all winter months, we obtained a SIT mean standard deviation of ≈ 20 cm (14 %) ± 11 cm (31–9 cm interval; see Table A3) and mean maximum deviation of ≈49 cm (35 %) ± 25 cm. These values increase to ≈27 cm (20 %) ± 9 cm for the standard deviation and ≈76 cm (58 %) ± 25 cm for the maximum deviation when we take into account “all snow products”. Therefore, compared to the deviation (roughly 30 cm) that can be extracted from the difference between ASD (red line) and the SIT max line (dotted blue) in Fig. 11, we obtain higher maximum deviations using this methodology. It suggests that considering large-scale global means could tend to underestimate impacts of snow depth on SIT and therefore the SIT uncertainty. Note that in addition to a more refined methodology, taking into account AMSR2B, which can reach high values, might have a small impact here as well.

Besides providing a relevant quantification of the impacts of snow depth on SIT, such a “multi-observation” approach also provides a space and time variability in uncertainties. It also avoids the necessity of relying on the strong, and usual, assumption of Gaussian propagation of de-correlated errors. These uncertainty estimates are very important for data assimilation, in which the characterisation of observation error is of main importance to constrain models (e.g. Stewart et al., 2008; Bunzel et al., 2016; Janjić et al., 2018). Note that there are a growing number of research programmes and international initiatives which aim at coordinating efforts for better characterising uncertainties and provide traceable quality indicators to sea ice earth observation (EO) products (e.g. World Meteorological Organisation (WMO) workshop at AWI, European Union/European Commission (EU/EC) Polar Science Week).

7 Discussions and conclusions

In this study, we have presented a Ka–Ku snow depth product computed, in both hemispheres, from the data obtained by altimeters on board the satellites SARAL (Ka-band) and CS-2 (Ku-band). This product has been compared over the period 2013–2019 with the other Ka–Ku snow depth product (DuST, only in the Arctic), AMSR2 passive radiometer snow depth estimations, W99m (only in the Arctic) and snow depth solutions of the models PIOMAS (GIOMAS in the Antarctic), MERCATOR and NESOSIM (only in the Arctic). Thereafter, all these products have been evaluated, in the Arctic, against OIB airborne data, IMB snow buoys and CryoVEx airborne campaign data, including Ka–Ku, as well as ALS and Ku, snow depths. In the Antarctic, the lack of validation data remains a major obstacle to the assessment of sea ice products. It is currently a primary focus of the ESA CSAO+ project.

In the Arctic, we observed a good agreement between the two Ka–Ku products (ASD and DuST) in terms of spatial distributions and annual and interannual variability, but DuST is always high biased by about 6.5 ± 0.5 cm (see Table A1). Although this is likely due to the re-calibration of the DuST product with OIB data mainly measuring thick snow, this bias probably also indicates variabilities in Ku-band penetration related to snow properties. Only the surface roughness is considered in the ASD calculation, but the impact of other snow properties should be analysed and considered.

The AMSR2B (and also AMSR2-NSIDC over FYI) data also tend to overestimate snow depths compared to ASD in the Arctic. This overestimation is lower than DuST, on average of about 3 cm, and does not seem to only be a bias since AMSR2B spatial patterns present noticeable differences. We reiterate that comparisons with AMSR2B have only been conducted in March and April because the data are not available over MYI for the other months. This absence of variability except in the MYI zone where AMSR2B has been re-calibrated is a surprising feature that needs to be considered when using these data.

Comparisons with models in the Arctic have highlighted that the MERCATOR model is always biased low when compared to the other products and simulates limited snow depth variability. This is mainly due to the use of the ERA-Interim forcing fields that underestimate precipitation. The impacts of other parameters such as the albedo or the amount of precipitation falling into the ocean should be investigated. In contrast, NESOSIM and PIOMAS spatial patterns are consistent with ASD and have a bias on average of a few centimetres high. We found comparable representations in PIOMAS and NESOSIM, although NESOSIM presents higher snow depths over MYI and higher variability. In spite of consistencies, these two models hardly represent the MYI sea ice drift in the Beaufort Gyre that can be observed in satellite observations. Note that we have not considered in this study the recent SnowModel-LG model that should consistently capture the Arctic snow depth over sea ice spatial variability (Stroeve et al., 2020; Liston et al., 2020). We refer to Zhou et al. (2021) for a comparison between the SnowModel-LG and various snow depth observation products.

In evaluation and for comparison of the snow depth products, we have compared OIB, CryoVEx and in situ IMB snow depths over the Arctic. Comparisons with OIB have shown satisfying consistencies with the ASD estimations. We found a RMSE of about 6 cm with a correlation coefficient of 0.66. By comparison, the DuST product has a RMSE of 8.7 cm and AMSR2B of 7 cm. Note that this higher RMSE for DuST, in spite of a re-calibration, probably comes from the OIB version used for re-calibrating. Indeed, DuST is re-calibrated against the OIB NOAA Wavelet Airborne Snow Radar dataset which is expected to present higher snow depth values. In addition, further analysis could be done in order to link the space and time variability in these biases with snow properties The various snow depth products are less consistent when comparing with IMB, although ASD mean values present a consistent order of magnitude. Due to the data coverage, validation with CryoVEx has been limited to one single track (where both Ku-band radar (ASIRAS), Ka-band radar (KAREN) and lidar scanner (ALS) were available). We observed that KAREN and ASIRAS snow depth estimations are low, including a non-negligible amount of negative values. ASD data always remain in between the laser-Ku (ALS-ASR) and Ka–Ku (Ka/ASR) CryoVEx snow depth estimations.

Since the results of such comparisons with validation data can vary from one methodology to another (e.g. grid sizes, smoothing kernel, dataset versions), they do not aim to assess the best snow depth product. However, they demonstrate how ASD provides a relevant snow depth solution, in good agreement with several validation data, and that ASD allows for the characterisation of the deviations between the different snow depth products. A more refined comparison with the CryoVEx data (including various tracks) is mandatory to understand the relative impacts of roughness and penetration and their link with snow properties. For instance, Willatt et al. (2011) show that the Ku-band dominant scattering surface can significantly vary with snow temperature, with a reduced penetration when temperatures increase. This feature could enhance an underestimation of ASD data with the warming of temperatures due to climate change. It also points out the difficulty to retrieve snow depth from altimetry beyond the 6 months of winter as variations in snow properties are most important in the summer period.

In the Antarctic, for satellite-based observations, spatial mean snow depth values of AMSR2-NSIDC and ASD are quite close. However, we found that AMSR2-NSIDC estimations are very low in zones of thin snow, while it tends to strongly overestimate snow depths in areas of thick ice. Such patterns of very thin snow have not been observed in the Arctic. Also patterns of thick snow are very likely due to the retrieval algorithm that is not adapted for thick ice. Regarding models, over the Antarctic, the MERCATOR model simulates mean values in better agreement with ASD. However, thicker snow patterns and the clockwise drifts in the Weddell Sea are not captured by the model. The GIOMAS model simulates very high snow depths nearly everywhere, likely leading to unrealistic snow depth representations. Sea ice models are generally very dependent on atmospheric forcings, and differences in snow accumulation in models can lead to very different snow depths. Model biases are largely due to the lack of consistent observations, which has, so far, limited the tuning of snow parameters such as albedo or snow diffusion coefficients (Chevallier et al., 2017; Uotila et al., 2019). These comparisons highlight snow properties in the Antarctic that are very different than in the Arctic. Such opposite model behaviours in the two hemispheres argue that a single model tuning of parameters may not be adapted for a global configuration. In the context of data assimilation, the positive impacts of the integration of SIT (Fiedler et al., 2021; Massonnet et al., 2015) could also be enhanced using ASD either by improving the SIT observations or by enabling the simultaneous use of freeboard and snow depth data to constrain the models.

An important limitation of the ASD snow depth data is its temporal coverage, imposed by SARAL which has only been available since 2013. To circumvent this limitation, we have shown that a simple monthly climatology, constructed as the average of ASD snow depth maps over the 2013–2019 period, could be a relevant option. More specifically, we have demonstrated that such a climatology could provide a more reliable solution than the W99m climatology, which strongly overestimates snow depths and is useless outside the central Arctic due to extrapolation. Indeed, comparison between the ASD climatology and the 2009–2012 OIB missions shows good consistency (comparable results of similar magnitude as when comparing the actual ASD product with OIB campaigns from 2014–2018), while W99m is on average biased about 14 cm high with an associated RMSE of 17 cm (when comparing with OIB). Such a climatology should be refined in order to provide relevant snow depth representations for at least the Envisat and CS-2 time period. For instance, the possibility of constructing a climatology using both ASD, AMSR and W99m information should be investigated.

A better snow depth representation significantly improves SIT estimates. In this context, this study has investigated the impacts of different snow depth products on SIT. Meanwhile, we have also presented a methodology to characterise the SIT level of uncertainty due to the snow depth deviations between products. The approach constructed an ensemble of SIT solutions using the different snow depth products over the 2013–2019 period, thus providing an estimation of the SIT level of uncertainty due to snow depth discrepancies and its space and time variability. Note that equivalent methodologies are frequently used in multi-model approaches. Taking only into account the satellite products AMSR2B, ASD, DuST and the W99m climatology over the 2013–2019 period, we found a spatially averaged mean standard deviation of 20 cm (14 % of the pan-Arctic mean SIT) and mean maximum deviation of 49 cm (35 % of the pan-Arctic mean SIT) between the different SIT estimations. Deviations between SIT estimations reached up to 77 cm (55 % of pan-Arctic mean SIT) in marginal and coastal zones.

We have highlighted in this study that the DuST and ASD products, albeit computed from data acquired by the same sensors at the same time, present some significant differences due to the methodology applied. The main difference is that for ASD we compare LRM data (PLRM for CS-2), while DuST compares LRM data with SAR data. The advantage of comparing LRM data together is that they have comparable footprint sizes, by which the difference reduces the effects of surface roughness to highlight the effects of volume backscatter and penetration. DuST, on the other hand, relies on the calibration of the measurements on OIB without taking into account effects of the surface roughness. In fact, differences in surface roughness and the effect it has on radar penetration and retrieval of surface elevation is currently one of the other main contributors to the SIT uncertainty budget (e.g. Hendricks et al., 2010; Landy et al., 2020), which still needs to be further investigated. We have only considered the roughness as mainly depending on the size of the footprint, but the backscattering coefficient variability is not as trivial and depends on several snow physical properties such as density, grain size, salinity and/or moisture (e.g. Lundberg et al., 2006; Adodo et al., 2018; Nandan et al., 2020). That is, specular surfaces have stronger nadir scattering of the radar pulse, affecting the shape of the radar echo and, in turn, the retrieval of the surface elevation. Empirical re-trackers like TFMRA have been shown to be sensitive over rougher surfaces (e.g. Ricker et al., 2014) In particular, compared to physical retracking they have been shown to overestimate sea ice freeboards over rougher ice while also underestimating sea ice freeboard over smooth, thin ice (e.g. Laforge et al., 2020; Landy et al., 2020). In turn, when assuming that snow depth can be retrieved from the difference between the Ku-band and Ka-band freeboards, the effect that surface roughness, and particularly the Ku-band volume retrodiffusion/penetration, has on the retrieval of the surface elevation warrants further investigation. An analysis of these dependencies is out the scope of this study, but it is necessary to further use the ASD data to improve our understanding of snow properties.

The use of two satellites, operating in different modes, with different orbits, also leads to important uncertainties and differences that cannot be neglected. The high-priority candidate mission, CRISTAL, which shall require a dual-frequency Ka–Ku SAR altimeter, could address some of these issues. The fact of having coincident and collocated SAR measurements obtained from the same platform, with similar ground footprints for both frequencies, will enable direct comparisons and a better understanding of the relative impacts of the surface roughness, radar penetration and volume backscattering. In this context, ASD demonstrates the possibilities of dual-frequency snow depth estimations, which could be further improved using physical retracking (as opposed to empirical re-trackers which have their own limitations, e.g. Ricker et al., 2014). In addition, CRISTAL will also mitigate the issue of the spatial coverage limitation of 81.5^∘ N imposed by SARAL in the Arctic. It is also important to note that the CRISTAL mission is the only altimetric mission planned to succeed CS-2 and avoid a gap of observations. Note that this Copernicus mission is not yet officially selected. If it is not confirmed, sea ice thickness measurement by radar altimetry above 81.5^∘ N will end with the CS-2 mission. Prior to potential CRISTAL data, comparisons with the very recent ICESat-2 and CS-2 snow depth estimations (Kwok et al., 2020) computed from the difference between laser and Ku-band measurements certainly provide important first results. Moreover, the recent change in CS-2 orbit to align with ICESat-2 (CRYO2ICE), providing collocated measurements with ≈ 3 h temporal delay, also allows for the investigation of Ku-band penetration and the discrepancies that different radar footprints can introduce.

Appendix A: Appendices

Table A1 Statistical comparison of the different snow depth estimations in the Arctic over the 6 winters of the 2013–2019 period. AMSR2's data are not indicated since only the months of March and April were available.

Table A2 Statistical comparison of the different snow depth estimations in the Antarctic over the seven winters of the 2013–2019 period. MIV signifies mean interannual variability, and MAV signifies mean annual variability.

Table A3 Statistical impact of the various snow depths on SIT in the Arctic over the 2013–2019 time period. The upper part considers all the snow depth products while the upper part consider only the products based on observations (ASD, DuST, AMSR2B and W99m). Note that AMSR2B SIT estimations on MYI are only taken into account for the months of March and April.

Lawrence et al. (2018)Rostosky et al. (2018)Meier et al. (2018)Schweiger et al. (2011)Lellouche et al. (2018)Petty et al. (2018)

Table A4 Time period of snow depth products used in this article. Note that the current data availability may have changed since this article was written.

Appendix B: List of acronyms

Acronym	Long name
ALS	Airborne laser scanner
ASD	Altimetric snow depth
AMSR	Advanced Microwave Scanning Radiometer
ARISE	Antarctic Remote Ice Sensing Experiment
ASIRAS	ESA Airborne Ku‐band Interferometric Radar
ASPeCt	Antarctic Sea Ice Process and Climate
ASR	Arctic System Reanalysis
AVISO+	Archiving, Validation and Interpretation of Satellite Oceanographic data
AWI	Alfred Wegener Institute
BROMEX	Bromine, Ozone and Mercury Experiment
CMEMS	Copernicus Marine Environment Monitoring Service
CNES	Centre national d'études spatiales
CPOM	Centre for Polar Observation and Modelling
CRISTAL	Copernicus Polar Ice and Snow Topography Altimeter
CRREL	Cold Regions Research and Engineering Laboratory
CryoSeaNICE	CryoSat SciEnce-oriented data ANalysis over Sea-ICE areas
CryoVEx	CRYOsat Validation Experiment
CS-2	CryoSat-2
CSAO+	CryoSat+ Antarctic Ocean
DTU	Technical University of Denmark
DuST	Dual-altimeter Snow Thickness
EASE	Equal-Area Scalable Earth
ECMWF	European Centre for Medium-Range Weather Forecasts
ERA-Interim	Re-Analysis Interim
ESA	European Space Agency
EUMETSAT	European Organisation for the Exploitation of Meteorological Satellites
FMCW	Frequency modulated continuous wave
FYI	First year ice
GCOM-W	Global Change Observation Mission – Water
GOP	Geophysical ocean product
HPCM	High Priority Copernicus Candidate Mission
ICESat-2	Ice, Cloud, and Land Elevation Satellite-2
JRA	Japanese Meteorological Agency
JAXA	Japan Aerospace Exploration Agency
LEGOS	Laboratoire d'Etudes en Géophysique et Océanographie Spatiales
LIM	Louvain-la-Neuve Sea Ice Model
LRM	Low-resolution mode
MERRA	Modern-Era Retrospective analysis for Research and Applications
MYI	Multi-year ice
N-ICE	Nansen Basin during the Norwegian young sea ice
NASA	National Aeronautics and Space Administration
MOSAIC	Multidisciplinary drifting Observatory for the Study of Arctic Climate
NCAR	National Center for Atmospheric Research
NCEP	National Center for Environmental Prediction
NEMO	Nucleus for European Modelling of the Ocean
NSIDC	National Snow and Ice Data Center
SMMR	Scanning Multichannel Microwave Radiometer
SMOS	Soil Moisture Ocean Satellite
SLA	Sea level anomaly

SSM/I	Special Sensor Microwave/Imager
SPICES	Space-borne Observation for Detecting and Forecasting Sea Ice Cover Extremes project
IFS	Integrated Forecast System
IMB	Ice mass balance
ODATIS	Online Data Extraction Service
OIB	Operation Ice Bridge
OSISAF	Ocean and Sea Ice Satellite Application Facility
PLRM	Pseudo low-resolution mode
PP	Pulse peakiness
SAR	Synthetic aperture radar
sd	Snow depth
SGDR	Sensor Geophysical Data Record
SIT	Sea ice thickness
SHEBA	Surface HEat Budget of the Arctic Ocean
TFMRA	Threshold first maximum re-tracker algorithm
UCL	University College of London
W99m	Warren-99 modified climatology