Arctic sea ice radar freeboard retrieval from the European Remote-Sensing Satellite (ERS-2) using altimetry: toward sea ice thickness observation from 1995 to 2021

Sea ice volume's significant interannual variability requires long-term series of observations to identify trends in its evolution. Despite improvements in sea ice thickness estimations from altimetry during the past few years thanks to CryoSat-2 and ICESat-2, former ESA radar altimetry missions such as the Environmental Satellite (Envisat) and especially the European Remote-Sensing Satellite (ERS-1 and ERS-2) have remained under-exploited so far. Although solutions have already been proposed to ensure continuity of measurements between CryoSat-2 and Envisat, there is no time series integrating ERS. The purpose of this study is to extend the Arctic radar freeboard time series back to 1995. The difficulty in handling ERS measurements comes from a technical issue known as the pulse blurring effect, altering the radar echoes over sea ice and the resulting surface height estimates. Here we present and apply a correction for this pulse blurring effect. To ensure consistency of the CryoSat-2, Envisat and ERS-2 time series, a multiparameter neural-network-based method to calibrate Envisat against CryoSat-2 and ERS-2 against Envisat is presented. The calibration is trained on the discrepancies observed between the altimeter measurements during the mission-overlap periods and a set of parameters characterizing the sea ice state. Monthly radar freeboards are provided with uncertainty estimations based on a Monte Carlo approach to propagate the uncertainties all along the processing chain, including the neural network. Comparisons of corrected radar freeboards during overlap periods reveal good agreement between the missions, with a mean bias of 0.30 cm and a standard deviation of 9.7 cm for Envisat and CryoSat-2 and a 0.20 cm bias and a standard deviation of 3.8 cm for ERS-2 and Envisat. The monthly corrected radar freeboards obtained from Envisat and ERS-2 are then validated by comparison with several independent datasets such as airborne, mooring, direct-measurement and other altimeter products. Except for two datasets, comparisons lead to correlations ranging from 0.41 to 0.94 for Envisat and from 0.60 to 0.74 for ERS-2. The study finally provides radar freeboard estimation for winters from 1995 to 2021 (from the ERS-2 mission to CryoSat-2).

1 Introduction

Several indicators illustrate the evolution of sea ice in response to climate change. Arctic sea ice extent has strongly decreased since the beginning of the satellite observation era by radiometry (Stroeve et al., 2012; Meier et al., 2014; Stroeve and Notz, 2018). The proportion of perennial ice has decreased significantly since 1984: the amount halved in April between 1984 and 2018 (Stroeve and Notz, 2018). The end of summer 2021 had the second-lowest amount of multiyear ice since 1985 (Meier et al., 2021). To improve our knowledge and forecast its evolution, an additional dimension becomes crucial: the thickness. Thick and old ice is disappearing and being replaced by younger, thin ice that has a higher mechanical sensitivity. Thin ice is more prone to deformation (Stroeve and Notz, 2018) that induces area changes and is more sensitive to climate hazards such as cyclones or strong winds (Rheinlænder et al., 2022). Thickness is a key parameter for sea ice study: it varies a lot according to the regions and modulates the sea ice volume evolution in the Arctic Ocean (Landy et al., 2022). Various campaigns have been carried out in the Arctic since the middle of the 20th century to measure sea ice thickness (Lindsay and Schweiger, 2013; Krishfield et al., 2014). However, these space- and time-limited measurements do not allow conclusions to be drawn from basin-scale sea ice volume variations. A quasi-global approach is possible through satellite altimetry, especially with radar altimetry, which is not impacted by the cloud cover and whose missions have been continuous since 1991.

Sea ice thickness estimation by spatial altimetry, in its modern form, was introduced by Laxon (1994) and Peacock and Laxon (2004) based on the freeboard methodology. The radar freeboard is obtained by taking the difference between the height measured above the floes and the height over leads interpolated below the floes. The radar freeboard has to be corrected for the radar signal slowdown within the snow layer to retrieve the ice freeboard, which is the thickness of the emerged part of the floe. Given the snow depth and the density of water, ice and snow, it is possible to derive the thickness of the ice assuming hydrostatic equilibrium. Lead and floe heights can be estimated using a heuristic retracker or, more recently developed, physical retrackers (Kurtz et al., 2014; Landy et al., 2019; Laforge et al., 2020). The implementation of this method requires the assumption that the Ku-band radar wave completely penetrates the snow layer, which is still widely discussed and is not the subject of a definitive consensus (Ricker et al., 2014, 2015; Nandan et al., 2017).

The launch of the CryoSat-2 (CS-2) mission featuring a high-resolution synthetic aperture radar (SAR) mode has enabled important advances in the estimation of sea ice thickness. The benefits are many: especially when compared with altimeters from past missions (European Remote-Sensing Satellite ERS-1 and ERS-2 and the Environmental Satellite – Envisat) operating with older technology, the low-resolution mode (LRM) has a larger surface footprint size, making thickness estimation more difficult. To reconstruct Envisat sea ice thickness estimation, Guerreiro et al. (2017), Paul et al. (2018) and Tilling et al. (2019) relied on the differences between Envisat and CS-2 during their common flight period to be able to calibrate the Envisat freeboard. Getting back to ERS missions, an additional problem coming from the instrument appears: the “pulse blurring” described in Peacock (1998) and in Peacock and Laxon (2004). Laxon et al. (2003) and Giles et al. (2008) published thereafter the first and last ERS thickness estimations for the Arctic and Antarctic sea ice so far (as a map averaging all the estimates of the different winters over the whole flight period of ERS-1 and ERS-2).

This study presents a method to recover a homogenous time series of the Arctic sea ice radar freeboard back to ERS-2 that is aimed for use in describing sea ice thickness changes over the last few decades. To minimize inter-mission bias along with the series, ERS-2 freeboard estimates are adjusted on Envisat radar freeboard estimates, which in turn were previously adjusted on CryoSat-2, taking advantage of the respective common flight periods. Consistency between missions is ensured by using the same processing chain regardless of the mission (before calibration), starting with the chosen retracking algorithm: the empirical threshold first-maximum retracker algorithm (Helm et al., 2014) with a threshold of 50 % (TFMRA50). Note that the term “radar freeboard” refers to the TFMRA50 radar freeboard in the KU band for both the LRM and SAR modes (depending on the mission). Since this LRM–TFMRA50 radar freeboard will be corrected to be consistent with CryoSat-2 and not conventionally obtained by calculating the difference in the height over floes and height over leads, it will be specified as neural network (NN) FBr, which stands for radar freeboard adjusted using the neural network. This study does not present any new understanding concerning LRM waveform retracking on sea ice and does not draw any conclusions about links between surface properties and TFMRA50 FBr from LRM but will focus only on recovering a long and homogeneous time series between different altimeter technologies. We also present the method used to correct the ERS-2 measurements for the effect of pulse blurring, which is a prerequisite for using ERS measurements over sea ice. The adjustment of LRM measurements on CryoSat-2 is performed using machine learning based on the surface state of the ice in Sect. 3. The associated uncertainties are derived using a Monte Carlo approach. Section 4 compares the monthly Envisat and ERS-2 radar freeboard data with various in situ, spaceborne datasets or other altimetry products available during this period. The time series is finally presented as a radar volume time series with trend estimation, providing a first overview of Arctic sea ice changes over the last 27 years.

2 Data

2.1 Satellite altimetry data

2.1.1 CryoSat-2

CryoSat-2 is an ESA altimeter mission launched in 2010. With a nearly polar and geodetic orbit, it enables observations of up to 88^∘ N, which makes it particularly adapted for cryosphere observations. Additionally, CryoSat-2 incorporates nadir SAR and synthetic aperture radar interferometric (SARin) technologies (Wingham et al., 2006). The two altimetry approaches exploit the Doppler capabilities of the instrument to reduce the along-track footprint from several kilometers to approximately 300 m compared to LRM (corresponding to a reduction in the footprint area from 5 to 180 km², Stammer, 2018). Increasing the along-track resolution of the aperture radar has led to considerable advances in estimating sea ice thickness. For this study, we use SAR and SARin data at 20 Hz from the ESA baseline-D L1b product. We derive the radar freeboard for the 7 coldest Arctic months (from October to April) of each year from November 2010 to the present.

2.1.2 Envisat

The ESA's Envisat mission was launched in 2002, reaching latitudes of 81.5^∘ N and 81.5^∘ S. The satellite carried the radar altimeter RA-2 (operating in LRM), with a high pulse-repetition frequency (PRF) of 1795 Hz allowing a large number of measurements per second to be performed, resulting in a better accuracy. The return pulses are averaged in batches of 100 to constitute each waveform (Roca et al., 2009). RA-2 L1b version-3 products including waveforms provided by the ESA are used in this analysis. Sea ice freeboard is computed for the coldest 7 months of each year from October 2002 to March 2012.

The surface illuminated by the satellite is significantly larger in LRM than in synthetic aperture radar mode (SARM) (by a factor of about 30). In addition, surface roughness determines the size of the illuminated footprint in LRM; the greater the roughness, the larger the footprint, whereas it is constant in SARM (Chelton et al., 1989; Raney, 1995). Therefore, surface roughness will have a greater effect on LRM range retrieval than the more nadir-focused SAR techniques (see Sect. 3.4 for more information). There is no waveform model for sea ice to account for the effect of roughness, and conventional retracking methods do not allow relevant radar freeboard estimation using LRM information alone. Our approach is to exploit these processing-mode differences to derive an LRM-corrected freeboard. To this end, we compare Envisat and CryoSat-2 datasets during the mission overlap period from November 2010 to March 2012 (see Sect. 3.4).

2.1.3 European Remote-Sensing Satellite (ERS-2)

In the 1990s, the ESA launched two ERS satellites, ERS-1 in July 1991 and ERS-2 in April 1995. ERS-1 was able to perform nominally until June 1996 and ERS-2 until November 2003. To extend the time series while ensuring continuity with the Envisat mission, ERS-2 products from the ESA Reaper project (Brockley et al., 2017) were used until July 2003. The ERS RA operated at a lower PRF than the Envisat RA-2 (1020 Hz against 1795 Hz, respectively). Thus, the 20 Hz waveforms are made up of 50 elementary echoes instead of 100 for RA-2. This leads to a higher speckle noise for ERS missions than for Envisat (see Sect. 3.5 for further details).

Another significant difference between RA and RA-2 comes from the tracker-board control loop that aims at centering the expected echo in the altimeter acquisition window. The delay between the transmission and the reception of the radar waveform depends on the vertical distance between the altimeter and the Earth's surface. This distance varies along with the satellite orbit and the ground topography. The time between the transmission of the radar wave and the opening of the acquisition window must therefore be constantly adapted (called the window delay (s) or the tracker range (m) if this is within a time or a distance). The distance or range between the altimeter and the measured surface is then equal to the sum of the tracker range and the epoch, i.e., the position of the waveform in the window. Since a waveform is an average of 50 individual pulses, it is important that each pulse is correctly centered in the window (to be aligned with the others). Otherwise, the resulting averaged 20 Hz waveforms will be blurred. This is unfortunately what happened to ERS altimeters over sea-ice-covered surfaces (Peacock and Laxon, 2004).

2.2 Ancillary data

Whether it is for the calculation of the radar freeboard itself, the LRM calibration or the comparison of our results to in situ data, we use various additional datasets. We present here additional datasets that have been used for this purpose.

The sea ice concentration field is needed to restrain the freeboard computation over a sea-ice-covered area. The product used is the NSIDC 0051 product based on Nimbus-7 Scanning Multi-channel Microwave Radiometer (SMMR) and Defense Meteorological Satellite Program (DMSP) Special Sensor Microwave/Imagers (SSM/I) Special Sensor Microwave Imager/Sounder (SSMIS) passive microwave data (Cavalieri et al., 1996). This product is also used to compute radar freeboard volume in Sect. 4.3 and for the LRM and SARM calibration. The study also requires a sea ice type product; this information is derived from the NSIDC 0611 sea ice age product (Tschudi et al., 2019) that is aggregated into two classes (multiyear ice, MYI, and first-year ice, FYI) according to the oldest age of the ice within the grid cell (FYI: ice age between 0 and 1 year; MYI: ice age of at least 1 year). Data are available as daily and weekly maps, respectively, with a 12.5 km grid resolution. The proportion of the MYI of a given grid cell refers, in this study, to the mean ice type observed by all the tracks (for each month of each mission) that pass within a 25 km radius of this grid cell. This value is computed during the gridding step. The proportion would consequently be overestimated compared to what can be estimated with ice-age-tracking algorithms.

SnowModel-LG (Liston et al., 2020a; Stroeve et al., 2020) is a snow depth product from a snow evolution model forced by different reanalyses: we use the version forced by ERA5 (Liston et al., 2020b). The dataset is available from 1 August 1980 and 30 July 2018 at a 25 km resolution. Although the most commonly used product is still the Warren snow depth climatology (W99) (Warren et al., 1999), it is no longer consistent for the recent period, and an altimetry-based product such as the altimetric snow depth (ASD) climatology (Garnier et al., 2021) would not be a relevant choice before the 2010s; we therefore justify the use of SnowModel-LG to ensure the relevance of the comparisons presented in Sect. 4.

2.3 Validation data

The results obtained in this study are compared with the different independent datasets presented in this section. Comparisons are detailed in Sect. 4. Most of the following datasets are included in the Lindsay and Schweiger (2013) dataset and are further described in the corresponding publication Lindsay and Schweiger (2015). Data availability is summarized in Fig. 1.

Figure 1 Summary of various available datasets for Envisat and ERS validation. Colors distinguish the different types of data: dark blue for satellite products, light blue for airborne data, yellow for submarines, orange for anchored moorings, green for buoys and red for direct measurements. (a) Temporal availability. (b) Spatial availability and extent of mission data gaps. The blue disks represent areas not covered by altimeters due to their orbit's inclination (no data above 81.5^∘ N for Envisat, 86^∘ N for ICESat-1 and 88^∘ N for CryoSat-2).

2.3.1 Airborne

Operation Ice Bridge (OIB) was a mission led by NASA. It consisted of airborne measurement campaigns using scanning lidar altimeter and snow radar to measure both snow depth and ice thickness (Kurtz et al., 2013). The data we use are from the Unified Sea Ice Thickness Climate Data Record of Lindsay and Schweiger (2013), Operation Ice Bridge Version 2 being processed by Kurtz et al. (2013). These measurements were carried out between 2009 and 2013 during each early spring or early fall near the coasts of the Canadian Arctic Archipelago and Alaska.

Airborne electromagnetic induction (AirEM) can measure total thickness (snow plus sea ice): the methodology is described in Haas et al. (2009). AirEM data that are used in this study are provided by Lindsay and Schweiger (2013) and are available from 2001 to 2013 for 22 campaigns in the Arctic Ocean and Fram Strait.

2.3.2 Moorings and submarines

The following data are all measured with upward-looking instruments that are installed either on anchored moorings or on board submarines. These instruments measure the sea ice draft, i.e., the height of the immersed part, from which the sea ice thickness can be derived for comparison with altimetry data.

The Beaufort Gyre Expedition Project (BGEP) is composed of a network of four moorings located in the Beaufort Sea (Krishfield et al., 2014). The moorings, equipped with upward-looking sonar (ULS), record drafts every 2 s with a precision evaluated to ±0.3 cm. Data are currently available from August 2003 to September 2018. The data were collected and made available by the BGEP based at the Woods Hole Oceanographic Institution.

Belter et al. (2020) performed and diffused a daily sea ice draft dataset based on upward-looking acoustic Doppler current profilers (ADCPs). These data are located in the Laptev Sea and are available from August 2003 to September 2016.

The Institute of Ocean Sciences (IOS) provides two ULS draft measurement datasets named IOS-Eastern Beaufort Sea (IOS-EBS) and IOS-Chukchi Sea (IOS-CHK) (Melling, 2008). For IOS-EBS, data are available from April 1990 to September 2003 in a network of nine sites in the Beaufort Sea. IOS-CHK is composed of data from a single site located in the Chukchi Sea between August 2003 and August 2005. The sea ice draft product for IOS-EBS comes from Melling (2008), and for IOS-CHK it comes from Lindsay and Schweiger (2013). Draft can be measured with a precision of about 0.05 m for young ice and can be overestimated up to 0.3 m for older and rougher ice.

The Davis Strait sea ice draft (Davis St) product from anchored moorings was detailed in Drucker et al. (2003). Data used in the study come from the Unified Sea Ice Thickness Climate Data Record (Lindsay and Schweiger, 2013) and are available from 2005 to 2008.

The Alfred Wegener Institut (AWI) mooring sea ice draft dataset is composed of 11 moorings in the Greenland Sea and Fram Strait processed by Witte and Fahrbach (2005). The data span the period from 1991 to 2002, and the draft is recorded with a 5 min frequency with an accuracy of ±0.2 m.

The last sea ice draft dataset presented in this section is derived from data collected by both U.S. Navy and Royal Navy submarines in the Arctic Ocean from 1975 to 2005 (National Snow and Ice Data Center, 2006; Wadhams and Horne, 1980; Wadhams, 1984; Wensnahan, 2005). It gathers data from 39 cruises. According to Rothrock and Wensnahan (2007), sea ice drafts are estimated to have an overall bias of 29 cm and a standard deviation of 25 cm from the actual draft.

2.3.3 Coastal stations

Environment and Climate Change Canada compiled weekly measurements from 27 monitoring stations along the coasts of the Canadian Arctic Archipelago in one product named CanCoast. Measurement methods can vary from one station to another (boreholes, hot-wire thickness gauges, etc.), but all the stations provide at least sea ice thickness and snow depth estimation with an accuracy of less than 1 cm.

2.3.4 Satellite altimetry products

Three satellite altimeter sea ice freeboard products have also been used for comparisons.

Guerreiro et al. (2017) presents the first Envisat radar freeboard dataset consistent with the CryoSat-2 mission. Radar freeboards are available as monthly maps from November to April between 2002 and 2012 between 65^∘ N and 81.5^∘ N. This dataset will be referred to as Envisat LEGOS-PP because LRM TFMRA50 FBr was corrected using the pulse peakiness (PP) and processed by LEGOS (Laboratory of Space Geophysical and Oceanographic Studies).

The product generated by the ESA Climate Change Initiative program for sea ice thickness estimation (SI-CCI) (Hendricks et al., 2018) includes the whole Arctic sea-ice-covered region for all winters (October–April) of the Envisat mission (2002–2012). It provides monthly grids of sea ice thickness, radar and ice freeboard combined with the related uncertainties. The corresponding methodology is described in Paul et al. (2018). In this study, this product will be named Envisat CCI.

Finally, ICESat-1 was a mission operated by NASA, launched in 2003 and ceasing operations in 2009. It was composed of a laser altimeter that allowed retrieval of the total sea ice freeboard (snow depth plus sea ice freeboard). The ICESat-1 product provides estimations for 15 periods of about 30 d between February 2002 and November 2008 at a 25 km grid resolution. The version used is the NASA Goddard one processed by Zwally et al. (2008) and Yi and Zwally (2009). The ICESat-1 total sea ice freeboard measurement accuracy is estimated to be about 0.05 m.

3 Methods

3.1 ERS pulse blurring correction

The values of the onboard tracker heights are h_trk, while ERS-2 flights over sea ice reveal instabilities of several meters that cannot be explained by the sea surface topography. This phenomenon is better observed by computing the surface height anomaly (h_rtrk) over all types of surfaces. This height anomaly is calculated according to Eq. (1), where alt is the altitude of the satellite, range is the onboard tracker range, epoch is obtained after the waveform retracking using TFMRA50, MSS is the DTU15 mean sea surface (Andersen et al., 2016), and geophysical_corrs is the sum of all the geophysical corrections. Figure 2 shows the TFMRA50-retracked height estimation for Envisat and ERS-2 along the collocated pass 25 for cycles 12 and 80, respectively (beginning of January 2003). Regarding the behavior of Envisat measurements, ERS-2 height anomalies show instabilities of about 1 m that make the measurements unusable.

$$\begin{array}{}\text{(1)}& h_{\mathrm{rtrk}}=\text{alt}-\text{range}-\text{epoch}-\text{MSS}-\text{geophysical\_corrs}\end{array}$$

Figure 2 Profiles of surface height anomaly over sea ice and ocean for pass 25 between 78 and 81^∘ N for Envisat in blue green (cycle 12), ERS-2 in blue and ERS-2 blurring corrected in orange (cycle 80). The red line represents the limit of the 50 % concentration of sea ice, as with the limit between the open ocean and an ice-covered area. The dark-blue line shows the location of the pass between Svalbard and Greenland. (a) The surface height along the latitude and (b) the probability density function of surface height for the three passes with the associated statistics, the average and the standard deviation (SD). The color legend is identical for both subfigures.

The instabilities of the height anomalies (mainly over sea ice) are known as pulse blurring and are a consequence of the onboard tracker settings. This phenomenon occurs for both the ERS-1 and ERS-2 missions. A simplified version of the ocean-mode tracking system is represented in Peacock (1998, p. 71). The tracking system is composed of three tracking loops to maintain echoes within the radar acquisition window: the height-tracking loop (HTL), the slope-tracking loop (STL) and the automatic gain control (AGC). The role of the HTL is to maintain the successive waveforms in the middle of the acquisition window. For this purpose, the tracking system is able to estimate the position of the tracker for each individual echo that comprised the average sequence. The tracker position is therefore adapted at 1020 Hz with a low-pass αβ filter described by Eqs. (2) and (3) according to the HTL error ε. T is the interval for the low-pass filter to be updated, T is ${\displaystyle \frac{\mathrm{1}}{\mathrm{PRF}}}\mathrm{s}$ for the HTL, and h_n is the tracker height for the nth echo, with n between 0 and 49. In the ocean mode, the algorithm used to estimate ε is a suboptimal maximum likelihood estimator (SMLE). Nevertheless, the SMLE has been developed for Brown-like waveforms that can be found over the open ocean, but it is not suitable for specular waveforms found over sea ice.

$$\begin{array}{}\text{(2)}& {\displaystyle }{h}_n=h_{n-\mathrm{1}}+\mathit{\alpha }\mathit{\epsilon }+T{\dot{h}}_n\text{(3)}& {\displaystyle }{\dot{h}}_n={\dot{h}}_{n-\mathrm{1}}+{\displaystyle \frac{\mathit{\beta }\mathit{\epsilon }}{T}}\end{array}$$

Both Eqs. (2) and (3) can be combined to give Eq. (4), which shows that the range window correction increases with a n² factor.

$$\begin{array}{}\text{(4)}& h_n=h_{\mathrm{0}}+nT{\dot{h}}_{\mathrm{0}}+n\mathit{\alpha }\mathit{\epsilon }+{\displaystyle \frac{n}{\mathrm{2}}}(n+\mathrm{1})\mathit{\beta }\mathit{\epsilon }\end{array}$$

A high ε (due to an inappropriate sea ice height error estimation algorithm, Roca et al., 2009) coupled with this low-pass filter could drive the large variation of range windows inside the same averaging sequence, which will “blur” the final averaged waveform echo due “the bad overlay between the individual echoes” (Peacock and Laxon, 2004). Because the height error is estimated at the end of each averaging sequence, a sudden change in the range window is then possible, especially since the error estimate is also affected by the pulse blurring and can explain tracker height oscillations. The problem mainly comes from the choice of the SMLE to estimate the HTL error. Indeed, it has been elaborated for ocean-like waveforms with a long trailing edge, in contrast to peaky waveforms whose power decreases suddenly after the maximum power peak.

A methodology was developed by Peacock and Laxon (2004) to deal with this issue. This method consists in finding a relation between the height error parameter ε and the difference between the measured surface (over an area covered by ice) and the “same area if it was not covered by ice” (Peacock, 1998). We interpreted this as the difference between the raw surface measurements and the interpolated ocean-level measurements Δh.

Figure 3 Difference between the raw surface measurements and the interpolated ocean-level measurements as a function of the height error parameter ε. Panel (a) is taken from Peacock and Laxon (2004), and panel (b) is a reproduction of cycle 83 of ERS-2. In panel (b), squares and circles, respectively, are the median and mean of Δh for each value of ε (on the x axis).

Figure 3 illustrates that our interpretation of the Peacock (1998) theory (on the right) fits with his results (on the left). Our results reproduce the linear relation between Δh and ε found by Peacock (1998) and Peacock and Laxon (2004) with a slope equal to ${\displaystyle \frac{-\mathrm{1}}{\mathrm{5}}}$ when ε is negative:

$$\begin{array}{}\text{(5)}& \left\{\begin{array}{cccc}{h}_{\mathrm{corr}}=& h-\frac{\mathit{\epsilon }}{\mathrm{5}},& \mathit{\epsilon }& \le \mathrm{0},\\ h_{\mathrm{corr}}=& h,& \mathit{\epsilon }& >\mathrm{0}.\end{array}\right.\end{array}$$

This correction is applied and presented in Fig. 2. The correction of the pulse blurring effectively reduces the instabilities of measurements. The correction is similarly asymmetric, so that the variations toward the positive height anomaly are more corrected than the others. The corrected surface height anomaly of ERS-2 now appears more similar to Envisat in terms of the noise and amplitude of variation. For this particular pass, the standard deviation has been reduced by 16 % and gets closer to Envisat's one. Figure A1 shows more results on the impact of blurring correction on ASA (all-surface anomaly, e.g., floes and leads) noise reduction compared to Envisat during a whole cycle (80 for ERS and 12 for Envisat).

3.2 Along-track radar freeboard retrievals

This section aims to describe the FBr processing chain. This procedure is common to all missions to preserve homogeneity and continuity.

The FBr is the difference between the sea ice surface height measured over floes and the sea surface height measured over leads. The use of satellite altimetry for these estimates was introduced by Laxon (1994) and Laxon et al. (2003). Sea ice freeboard (FBi) can be derived from FBr after a correction of the lower wave propagation speed into the snow layer (Kwok, 2014). In this study, we focus on the radar freeboard (without the speed propagation correction) to avoid the introduction of errors related to snow depth estimation.

It is therefore necessary to discriminate leads and sea ice floes. This essential step is based on the peakiness of the radar waveform that quantifies the specularity of the surface (Laxon et al., 2003; Peacock and Laxon, 2004). We use the definition of Guerreiro et al. (2017) (Eq. 6). Pulse peakiness thresholds can depend on the radar altimeter. To ensure the continuity between Envisat and ERS-2, these thresholds have been adapted to keep the same lead–floe proportion during their common flight period; see Table A1.

$$\begin{array}{}\text{(6)}& \mathrm{PP}={\displaystyle \frac{max\left(\mathrm{WF}\right)}{{\sum }_{i=\mathrm{0}}^{{\mathrm{Nb}}_{{\mathrm{WF}}_{\mathrm{bins}}}}{\mathrm{WF}}_i}}\end{array}$$

Ranges are estimated using the TFMRA50 retracker for all surfaces and all missions mentioned here to maintain inter-mission continuity. Height anomaly measurements are then expressed relative to the DTU15 mean sea surface (Andersen et al., 2016) and corrected for the common geophysical corrections (oceanic, polar, solid earth and load tides as well as tropospheric and ionospheric corrections) according to Eq. (1). The tide model used is the FES14 from Carrere et al. (2015). The ERS-2 pulse blurring correction (detailed in Sect. 3.1) must be applied at this stage of the freeboard computation.

Surface height anomaly is then split into two variables using the pulse peakiness classification explained above, the ice-level anomaly (ILA) over floes and the sea-level anomaly (SLA) over leads. ILA and SLA outliers are removed by filtering data that are outside the interval: a rolling mean ±3 rolling standard deviations with a 60 km large sliding window. After filtering, ILA and SLA are smoothed using a rolling mean at 12.5 km, and then SLA and ILA are linearly interpolated (including below floes for SLA and above leads for ILA) and are again smoothed using a rolling mean at 12.5 km. No limit of distance is used to discard the radar freeboard, but the interpolation, smoothing and filtering are not done between values separated by land. Indeed, the processing is done within ocean segments separated by land in order to isolate statistics between segments.

In this study, we will only use the FBr measurements that are made over floes: indeed, the LRM data correction, explained in Sect. 3.4, is based on floe characteristics.

3.3 Data gridding

We performed the LRM FBr correction using monthly maps in EASE2 (Brodzik et al., 2012) with a 12.5 km resolution. The FBr gridding is done by averaging the values within a 25 km radius from each pixel weighted by the inverse of the uncertainty; for other variables, weighting is not applied. Only sea ice concentrations above 50 % are considered. To limit the outliers of LRM FBr at the sea ice–ocean boundary, data with ice concentrations lower than 85 % are removed when the waveforms have a leading-edge width (LEW) of higher than 2.5 gates.

3.4 Correction of LRM radar freeboards against SARM freeboards using neural networks

The radar freeboard maps obtained from the process presented in Sect. 3.2 for Envisat and CryoSat-2 are shown in Fig. 4.

Figure 4 Pan-Arctic radar freeboard maps for January 2012 for (a) Envisat uncorrected and (b) CryoSat-2.

Important differences between Envisat and CryoSat-2 can be noticed in terms of both patterns and mean values. Negative radar freeboards are mainly due to the retracker choice. Indeed, TFMRA50 is used to retrack heights on both leads and floes; this introduces a bias in the height over leads. The TFMRA threshold to retrack heights over leads should be closer to 80 %, and the use of a 50 % threshold corresponds to the position of the retrack point for ocean surfaces, not specular ones (Poisson et al., 2018). The surface over leads is measured higher than it is and even higher than the surface over floes. The SLA bias (over leads) is evaluated as constant for the SARM altimeter in the study of Laforge et al. (2020); this conclusion is also relevant for LRM altimeters as waveforms over leads are peaky and similar from one lead to another. This positive constant bias over leads results in a negative bias in the radar freeboard. To avoid this bias, the retracker threshold could be adapted for leads or the SLA could be corrected on the CryoSat-2 one. Nevertheless, a threshold of 50 % ensures the stability of the range (Poisson et al., 2018, Fig. 9), in contrast to higher thresholds (80 %–95 %) that could lead to up to 47 cm random error in the SLA. A TFMRA at 50 % for both leads and floes is preferred in this study as a constant bias is easier to correct than an undetermined random error. Nevertheless, beyond this bias, a lack of representative sea ice patterns can be observed. For instance, thick ice regions do not appear for the Envisat mission (see Fig. 4). This phenomenon also appears in the ERS-2 FBr since it is also in LRM with the TFMRA50 retracker. This inconsistency comes from the LRM itself (plus the retracker choice), which has a larger footprint than the resulting one in SARM. The larger footprint leads to a large impact of the surface roughness on the reflected echo. It is therefore impossible to distinguish between the contributions of heights and roughness without a physical model of the effect of sea ice roughness on the reflected echo.

Paul et al. (2018) and Guerreiro et al. (2017) showed that there is a significant correlation between the patterns of the parameters characterizing the surface roughness and the differences between the FBr of CS-2 and Envisat. These observations led to the development of correction methods of Envisat relative to CS-2, taking advantage of the mission-overlap period. In the same way, we propose correcting ERS-2 against Envisat and, once corrected, using the common flight period between Envisat and ERS-2.

So far, several empirical methods of Envisat freeboard correction using CryoSat-2 have been developed (Guerreiro et al., 2017; Paul et al., 2018; Tilling et al., 2019). In Guerreiro et al. (2017), the correction consists in finding a link between Envisat and CS-2 freeboard differences and the PP of Envisat's waveforms. The correction proposed by Paul et al. (2018) adapts the TFMRA threshold of Envisat waveform retracking according to the LEW of the waveforms and the surface backscatter. Tilling et al. (2019) suggested correcting the Envisat sea ice thickness (SIT) by utilizing the distance between leads and floes. All the methods are based on the comparison of monthly gridded data. The first two studies share the common approach of trying to correct the Envisat FBr using the surface roughness (characterized by one or more parameters as proxies). They both propose a third-degree polynomial function to link Envisat and CS-2 FBr differences to surface properties. Paul et al. (2018) were the first to propose using two distinct parameters that characterize two roughness scales to correct LRM measurements that can impact the waveform shape differently. Our method follows the same approach (pulse peakiness, leading-edge slope) with other additional parameters that define the sea ice state, such as ice concentration, ice type or season. The leading-edge slope refers to the leading-edge height divided by the LEW computed between 30 % and 70 % of the maximum power of the first waveform peak. The LRM-correction model procedure is based on a neural network in order to manage strong nonlinearities. The procedure is illustrated in Fig. 5.

Figure 5 Diagram illustrating the principle of freeboard correction by a neural network with the two main steps: in panel (a), the neural network training phase, and in panel (b), the prediction (correction) phase. Ω corresponds to the inputs and Γ to the output of the neural network. FB: radar freeboard; Conc: sea ice concentration; LES: leading-edge slope; PP: pulse peakiness, f_MYI: MYI proportion.

For each mission to be corrected, the NN is trained on the common flight period between the mission considered a reference and the one to be corrected. It takes as inputs monthly grids of the following parameters: LRM FBr (to be corrected), pulse peakiness, leading-edge slope, ice concentration, MYI proportion, the period and, as a target, the reference FBr (SARM FBr or LRM-corrected FBr). Note that the period of the year is taken to capture the seasonal variability better, as snow on sea ice and sea ice physical properties change along the seasons. Inputs and targets are standardized before the training step.

The neural network used is a multilayer perceptron (MLP). LRM FBr correction was performed with scikit-learn (Pedregosa et al., 2011). The MLP is composed of four hidden layers, each composed of 100 neurons. The choice of hyperparameters (number of neurons, the learning rate, the regularization term, the batch size, activation functions, solver for the weight optimization) was made using a grid search method. The evaluation criterion, called the score, is chosen as the determination coefficient. Models are trained on 90 % of the dataset and tested on the remaining 10 %, the splitting being random. During the tuning step, models are cross-validated, which means that they are each trained five times with the same combination of hyperparameters but without the same training or test dataset. The five scores are then analyzed to determine the best combination. Cross-validation gives a better idea of the model performance as the dependence on the training dataset is limited. The activation function for the hidden-layer neurons is a sigmoid related by possible negative radar freeboard values and the optimizer Adam (Kingma and Ba, 2014). Moreover, in order to avoid overfitting, an early stopping criterion is used to stop the model training as soon as the score is not improved during 10 consecutive iterations, with a defined tolerance (the validation fraction is set to 10 % for the early stopping).

Finally, once the hyperparameter combination is set (see Table A3 for the hyperparameter selection), the MLP is trained on the whole dataset to provide the correction function. The trained model is then applied to the LRM monthly grids to obtain a monthly LRM-corrected radar freeboard.

3.5 Radar freeboard uncertainty quantification

This section aims to estimate the uncertainties for Envisat and ERS-2 FBr estimations. The uncertainty budget is split into two steps corresponding to the two main parts of the freeboard processing chain. The first step covers the along-track processing up to the gridding, which is common to all missions, and the second step concerns the correction of the LRM freeboard, which consists in predicting the corrected FBr with the neural network (NN FBr).

Figure 6 Summary diagram of the uncertainty budget during the along-track, gridding and training correction steps. Panel (a) corresponds to the along-track to grid uncertainty budget. Panel (b) defines the notations for the Monte Carlo procedure: Ω for the neural network input parameter and Γ for the neural network output parameter (radar freeboard) with σ_Ω and σ_Γ the corresponding uncertainties. Panel (c) corresponds to the training of M models with noisy inputs and outputs. Panel (d) shows the predictions of the N noisy input with the M trained neural network. γ is the predicted radar freeboard estimation for one pixel of the M×N predictions. M=100, and N=200.

The uncertainty budget methodology concerning the first part is taken from Landy et al. (2020) and Ricker et al. (2014). We assume that for this step there are three sources of uncertainty. Two of them are random uncertainties: the speckle noise, largely discussed in Wingham et al. (2006), and the accuracy of the SLA measurement. The last one is linked to both retracker choices, surface roughness and snow radar signal penetration (Ricker et al., 2014).

According to Wingham et al. (2006), the speckle noise generates an “error” (σ_l1) from 7 to 14 cm (depending on the acquisition mode) in the range measurement for the CryoSat-2 altimeter. Estimations of the speckle noise error in the range for other missions can be relied on for the individual echo number used to compute the averaged waveforms (for LRM sensors). Indeed, the speckle noise generates an error in the range as a function of $\sqrt{n}$, where n is the ratio between the number of individual pulses used for each averaging sequence for CS-2 and for the mission for which we want to estimate the error (Calafat et al., 2017) (see Table A2 for LRM RA characteristics). All final σ_l1 values are summarized in Table 1. Note that Wingham et al. (2006) use the terminology “error”, but this refers to uncertainty.

Table 1 Summary of the “range” error for various missions governed by the speckle noise (Wingham et al., 2006).

^* Wingham et al. (2006).

SLA uncertainty (σ_SLA) estimation is taken from Ricker et al. (2014). The uncertainty in the SLA is estimated as the standard deviation of the SLA within a sliding window of 25 km if there are some leads within this window. If not, the SLA uncertainty is taken as the difference between the interpolated and smoothed SLAs and the mean SLA computed as the mean of raw SLA measurements at leads within a segment of the ocean (if the pass is over land, statistics are made segment of ocean by segment of ocean). Finally, we consider that σ_SLA and σ_l1 are not correlated and can be combined to give the random part of radar freeboard uncertainties (σ_R) following Eq. (7).

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{R}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{SLA}}^{\mathrm{2}}+{\mathit{\sigma }}_{l\mathrm{1}}^{\mathrm{2}}\end{array}$$

The radar freeboard (including uncertainties) gridding methodology is taken from Ricker et al. (2014, Sect. 2.4) in order to take into account the random uncertainties in the radar freeboard gridding process.

In Landy et al. (2020), the FBr systematic uncertainty budget is decomposed into two parts: on the one hand, the uncertainties due to the penetration of the signal in the snow (depending on its salinity or whether it is composed of metamorphic snow according to the type of ice) and, on the other hand, the surface roughness. We assume, as in Ricker et al. (2014), that the comparison of the freeboard from different retrackers does not enable us to separate the contribution of the roughness from the signal partial penetration. We therefore consider both sources to be one mixed contribution, estimated as about 20 % and 30 %, respectively, of the sea ice thickness for FYI and MYI (Landy et al., 2020). The systematic uncertainties can be underestimated, as the penetration of the radar waves in the snow uncertainty may be poorly handled. Note that this systematic uncertainty budget only concerns the CS-2 mission, which is afterward propagated to Envisat and ERS-2. Indeed, other missions will be corrected based on CS-2 estimates.

Therefore, we combined random and systematic uncertainties with a quadratic sum to get the total radar freeboard uncertainty in a grid and for the related radar freeboard estimation. The uncertainty of the other inputs (LES, PP, sea ice concentration, MYI proportion) is considered to be, for each grid cell, 2 times the standard deviation of the measurements used to calculate the average value (grid cell value) divided by the number of tracks passing through the corresponding grid cell.

As explained in Sect. 3.4, the LRM radar freeboard correction is predicted by a neural network. The uncertainty propagation through the neural network is not straightforward, since the inputs and outputs of the model are not linked by a mathematical relation. We have chosen to estimate the propagation of uncertainties through a Monte Carlo approach. This method allows us to propagate uncertainties through non-analytically represented systems. Nevertheless, the Monte Carlo method requires a representation of the distribution of the input and target parameters. As discussed above, there is no thorough knowledge of the distribution of the input data for the NN. The method consists in training a number M of NNs with noisy inputs. The noise has been added to all inputs (for each grid cell and each month) according to a Gaussian distribution centered on the estimated value and the corresponding uncertainty as the standard deviation. Training and predictions are done for all the noisy input and output, and then the distribution of M×N radar freeboard predictions (from the M noisy NN models applied to N noisy inputs) is analyzed for each grid cell and each month. The whole uncertainty budget process is summarized in Fig. 6.

Unfortunately, we have not been able to identify known distributions such as normal, log-normal or gamma for predicted FBr distributions. However, we can still derive various statistics, such as the quantiles at 2.5 % and 97.5 %, which represent 95 % of the values or the standard deviation to describe the FBr distribution for each pixel of all monthly grids.

4 Results

This section first presents the correction performance when applied to Envisat with respect to CS-2 and thereafter to ERS-2 with respect to Envisat. After this we present the Envisat and ERS-2 freeboard comparisons against independent validation datasets.

4.1 Correction performances

The LRM correction methodology is presented in Sect. 3.4. It was successively applied to Envisat (Env) with respect to CS-2 and then to ERS-2 with respect to Envisat NN FBr. Figure 7 compares Envisat NN FBr and CS-2 FBr during December 2010 and April 2011. Figure 8 presents the same feature for Envisat and ERS-2 NN FBr during December 2002 and April 2003.

Figure 7 Comparison of Envisat NN FBr against CryoSat-2 FBr for December 2010 in the upper half and April 2011 in the lower half. Maps (a) and (g) refer to Envisat with the corresponding CryoSat-2 radar freeboards (b) and (h). Maps below panels (d), (e), (j) and (k) are the related uncertainties. The right column presents freeboard difference maps (Env − CS-2) (c, i). Panels (f) and (l) show the distribution of Envisat FBr in red, CryoSat-2 FBr in blue and ΔFBr in grey. Histograms only include common data between Envisat and CryoSat-2, and data north or 81.5 ^∘ N are excluded. μ refers to the mean difference and SD to the standard deviation of the difference.

Figure 8 Same as Fig. 7 but for ERS-2 and Envisat during December 2002 and April 2003. Histograms only include data for the coinciding region between ERS-2 and Envisat. μ refers to the mean difference and SD to the standard deviation of the difference.

As can be noticed on these maps, Envisat correction allows recovery of typical patterns of the Arctic sea ice with thick ice near the coasts of the Canadian Arctic Archipelago and Greenland and thinner ice in the eastern part of the basin.

Nevertheless, compared to CS-2, the radar freeboard of thick ice is slightly underestimated and thin ice is slightly overestimated. However, the mean difference between the two products is close to zero (1 or 2 mm) for both months. We also note that the correction results in an asymmetric distribution with a tendency toward a log-normal distribution at the basin scale, whereas the CS-2 distribution for the same mask is centered and appears Gaussian-like.

Radar freeboards of ERS-2 and Envisat (see Fig. 8) are even more similar. Again, the bias is negligible, and the standard deviation (SD) of their difference is 2 to 3 times lower than for Env and CS-2 with a SD of 0.03 and 0.05 m.

The better performance of the ERS-2 correction can be explained because ERS-2 and Envisat carried similar LRM altimeters, and they flew on the same orbit with a 28 min delay between the two missions that allowed them to observe nearly the same surface. On the other hand, CryoSat-2 operates in SARM and flies on a quasi-polar orbit with cycles and subcycles very different from Envisat ones. These comparisons show that the neural network correction gives satisfying results, at least during the common flight periods. The correction function is then applied to the 10 years of the Envisat mission and the 8 years of the ERS-2 mission.

Tables A4 and A5, respectively, represent statistics for the Envisat–CryoSat-2 and ERS-2–Envisat radar freeboards for each month of the mission-overlap periods. For both corrections, the averaged radar freeboards are close. The highest mean difference reaches 7 mm in February 2011 for Envisat, i.e., 9.5 % of the mean Envisat NN FBr. Concerning ERS-2, the mean difference between ERS-2 and the Envisat NN FBr does not exceed 3 mm, or 3.3 % of the ERS-2 mean NN FBr. Concerning all the overlap periods, the mean FBr difference is 3 mm for Env and CS-2, 4.1 % of the Envisat mean NN FBr, and −2 mm for the ERS-2 and Env one, about 2.2 % of the ERS-2 mean NN FBr.

In both Figs. 7 and 8, the uncertainties presented for Envisat and ERS-2 are in fact 2 times the standard deviation of the respective NN FBr distributions as output of Monte Carlo simulations for each pixel grid. Detailed statistics for uncertainties are also provided in Tables A4 and A5.

The median uncertainty in the radar freeboard for the period 2010–2012 is 6.3 m for Envisat and 2 cm for CS-2. Regardless of the month, the mean and median uncertainties of Envisat are always larger than those of CS-2. Concerning the period 2002–2003, the median uncertainty is 8 cm for the ERS-2 radar freeboard and 7.3 cm for Envisat. Similarly, statistics on uncertainties are globally higher for ERS-2 estimates (see Tables A4 and A5 for detailed statistics).

4.2 Validation

In this section, Envisat and ERS-2 NN FBr are evaluated against a large set of independent data. These data are presented in Sect. 2.3 and include in situ, airborne and space-based measurements providing sea ice freeboard, draft or thickness. Data have been converted to sea ice thickness to make it comparable, except for the SI-CCI and LEGOS-PP Envisat products, which also provide FBr. The conversions are based on the hydrostatic balance assumption of sea ice covered by snow in seawater that is described in Eq. (8).

$$\begin{array}{}\text{(8)}& \mathrm{SIT}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{i}}}}\cdot {\mathrm{FB}}_i+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{i}}}}\cdot h_{\mathrm{s}}\end{array}$$

The snow depth (h_s) is taken from the dataset itself if given (e.g., OIB and CanCoast) with a constant snow density ${\mathit{\rho }}_{\mathrm{s}}=\mathrm{300}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. Otherwise, the SnowModel-LG and ERA5 version is used for the snow load. The water density values used are ${\mathit{\rho }}_{\mathrm{w}}=\mathrm{1024}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, ${\mathit{\rho }}_{\mathrm{i}}^{\mathrm{FYI}}=\mathrm{917}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ and ${\mathit{\rho }}_{\mathrm{i}}^{\mathrm{MYI}}=\mathrm{882}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ (Alexandrov et al., 2010). The sea ice density depends on the MYI proportion within the grid cell, such as described by Eq. 9).

$$\begin{array}{}\text{(9)}& {\mathit{\rho }}_{\mathrm{i}}=(\mathrm{1}-f_{\mathrm{MYI}}){\mathit{\rho }}_{\mathrm{i}}^{\mathrm{FYI}}+f_{\mathrm{MYI}}{\mathit{\rho }}_{\mathrm{i}}^{\mathrm{MYI}}\end{array}$$

For Ku-band measurements, the speed reduction in the wave in the snow layer is taken into account to obtain the ice freeboard ${\mathrm{FB}}_{\mathrm{i}}={\mathrm{FB}}_{\mathrm{r}}+h_{\mathrm{s}}(\frac{c}{c_{\mathrm{s}}}-\mathrm{1})$ (Mallett et al., 2020), with c the speed of light in a vacuum and c_s the speed of light in the snow estimated as $c_{\mathrm{s}}=c(\mathrm{1}+\mathrm{0.00051}{\mathit{\rho }}_{\mathrm{s}}{)}^{-\mathrm{1.5}}$ and determined by Ulaby et al. (1986).

The different datasets are then gridded into monthly EASE2 grids of a 12.5 km resolution so as to facilitate the comparison to Envisat and ERS-2 monthly grids. Concerning mooring or coastal measurement stations, data are averaged to get one value per month.

4.2.1 Envisat

Although numerous datasets are available during the Envisat flight period, the spatiotemporal coverage of the Arctic basin remains very patchy. The following comparisons are presented with several types of sensors to reinforce the relevance of the validation. The consistency and discrepancies are discussed in the following.

Figure 9 Comparative scatterplots between Envisat NN sea ice thickness or radar freeboard estimations and other datasets. The x axis indicates the sea ice thickness from the (a) OIB total ice freeboard, (b) AirEM snow plus ice thickness, (c) CanCoast ice thickness, (d) UK and US submarine draft and (e) ICESat-1 total freeboard. Panel (f) compares our Envisat radar freeboard with the SI-CCI Envisat solution and panel (e) with the Envisat PP solution from Guerreiro et al. (2017). Color bars represent the normalized density. log ₁₀ was applied before the normalization for panels (e) and (f) due to the large number of data. N is the number of the couple of values that are compared, “Med” refers to the median, SD refers to the standard deviation, RMSE refers to the root mean square error, and r refers to the correlation coefficient.

Figure 9 gathers comparisons between Envisat and different datasets coming from airborne ones, spaceborne ones, submarines, drifting buoys and fixed coastal stations, and Fig. 10 presents comparisons with fixed moorings.

Figure 10 Comparative scatterplots between Envisat NN sea ice thickness estimations and anchored mooring datasets. Each dot corresponds to a monthly averaged value. The x axis indicates the sea ice thickness from (a) BGEP, (b) BGEP vs. Envisat CCI, (c) Davis Strait, (d) IOS CHK and EBS and (e) Transdrift Laptev Sea ice draft. The color bar shows the MYI proportion. N is the number of the couple of values that are compared, “Med” refers to the median, SD refers to the standard deviation, RMSE refers to the root mean square error, and r refers to the correlation coefficient.

Comparisons presented in Fig. 9 provide satisfying statistics, with correlation values between 0.41 (CanCoast) and 0.71 (Envisat SI-CCI). Correlation between Envisat and OIB is in good agreement with Kurtz et al. (2014), who showed that the autocorrelation of OIB varies from 0.46 to 0.60. Nevertheless, AirEM data are poorly correlated with Envisat, and statistics reveal a high bias and root mean square error (RMSE) (up to 1 m). Disregarding CanCoast and spaceborne estimations, comparisons reveal a negative bias with Envisat from −3.9 to −2.9 cm, which could suggest an underestimation of Envisat sea ice thickness. The relevant statistics with CanCoast, with a bias of 8.9 cm and a RMSE of 6.4 cm, suggest that this underestimation could be attributed not only to Envisat, but maybe also to snow depth or other parameters. However, bias remains within the estimated range of uncertainty. The bias between OIB and Envisat estimation could also be attributed to the OIB snow depth whose estimation is sensitive to the algorithm used (Kwok and Haas, 2015; Kwok et al., 2017).

The comparison with ICESat-1 data reveals a strong dispersion with a low bias of 20 cm and a correlation of 0.50. The Envisat radar freeboard product established in the framework of the CCI and the version presented in this study are coherent, with a bias close to zero, a standard deviation from the difference of 5.2 cm and a fairly high correlation of 0.71. The Envisat CCI and Envisat NN datasets are consistent with a low mean bias, whereas the Envisat PP radar freeboard from Guerreiro et al. (2017) presents a thinner mean FBr than Envisat NN estimations.

The Envisat NN solution is also compared to several mooring (BGEP, IOS, Davis Strait and Transdrift Laptev Sea) datasets (Fig. 10). The four campaigns yield fairly high correlations, with Envisat NN estimates greater than or equal to 0.62 and reasonable standard deviations of 50–60 cm down to 32 cm for the Laptev Sea estimates from Transdrift. While Envisat has a low negative bias of −9.2 cm with respect to BGEP, other campaign biases are largely negative, with values from −53 to −76 cm.

Figure 10b also provides a comparison between the Envisat SI-CCI version and BGEP, similarly to Fig. 10a. With respect to BGEP estimations, the statistics are slightly better for the product presented in this study than for SI-CCIs. Nevertheless, the two products seem to be relatively consistent with each other.

4.2.2 ERS-2

The ERS-2 freeboards are also compared with measurements from AirEM, US and UK submarines, and CanCoast in Fig. 11.

Figure 11 Comparative scatterplots between ERS-2 NN sea ice thickness estimations and three in situ datasets. The x axis indicates the sea ice thickness from (a) AirEM total thickness, (b) UK and US submarine draft and (c) CanCoast sea ice thickness. The color bar indicates the normalized density. N is the number of the couple of values that are compared, “Med” refers to the median, SD refers to the standard deviation, RMSE refers to the root mean square error, and r refers to the correlation coefficient.

Figure 12 Comparative scatterplots between ERS-2 NN sea ice thickness estimations and two anchored mooring datasets. The x axis shows sea ice thickness estimations from (a) IOS Beaufort Sea and (b) AWI mooring sea ice drafts. The color bar indicates the respective MYI proportion. N is the number of the couple of values that are compared, “Med” refers to the median, SD refers to the standard deviation, RMSE refers to the root mean square error, and r refers to the correlation coefficient.

The results of ERS-2 validation are close to those of Envisat. Figure 11 shows similarly strong discrepancies with AirEM but even better results for CanCoast and submarine comparison in terms of correlations (0.60 and 0.74, respectively, instead of 0.41 and 0.55 for Envisat), and standard deviations of the differences are equivalent to those found for Envisat (62 and 57 cm compared to 63 and 53 cm). The comparisons between Envisat and the mooring illustrated in Fig. 12 are also relevant, with correlations of 0.67 and 0.66. Similarly to Envisat, comparisons with ULS data reveal non-negligible negative biases. As a comparison, the bias between CryoSat-2 and OIB between 2010 and 2019 is about 16 cm, and the RMSE is about 77 cm. Concerning CS-2–BGEP comparisons, the bias is 21 cm with the same overestimation of FYI thickness for CS-2, and with Transdrift Laptev Sea comparisons this shows a negative bias of −38 cm.

Although draft measurements by moorings are among the most accurate measurements (they measure 90 % of the total thickness), they will tend to overestimate the true thickness when the ice bottom surface is rough, which is inherent to the method. Indeed, the chaotic aspect of the lower surface of sea ice can impact the ULS-returned echoes. This strong deformation concerns mainly the thick and rough ice, which can explain the tendency of ULS measurements to overestimate thick ice relative to Envisat and ERS-2. Eventually, the methods used for the above comparisons can be questioned insofar as monthly averages are compared with punctual measurements (spatially and temporally), which may indeed induce biases (e.g., OIB, AirEM). Conversions from radar freeboard to sea ice thickness can also be suspected in bias comparisons, as this depends on the snow depth product. SnowModel-LG has been chosen because it is the only continuous dataset of the Arctic basin-scale product that covers more than the 27 years required in this study.

4.3 Radar freeboard volume time series from 1995 to 2021

The time series represented in Fig. 13 is derived from monthly maps processed as developed in Sect. 3 and validated in Sect. 4.2. Nevertheless, even after the gridding process, missing values can occur, especially where track density becomes low (e.g., close to the ice–ocean boundary). To ensure comparable volumes from one month to another, missing data have been replaced by interpolated values from the Gauss–Seidel relaxation method implemented in pangeo-Pyinterp (https://github.com/CNES/pangeo-pyinterp, last access: 15 January 2022), a Python library developed by CNES (the French National Space Agency). It is important to note that the volumes presented in Fig. 13 only consider values up to 81.5 ^∘ N (ERS-2 and Envisat orbit limitation).

Figure 13 Time series representing radar freeboard volume up to 81.5^∘ N for each winter month for ERS-2 in dark blue, Envisat in teal and CS-2 in dark red. Blue triangles are winter mean radar freeboard volumes. Red lines are linear regressions of winter mean volumes from 2002/2003 for the dashed line and 1995/1996 for the solid line, and estimated trends are $-\mathrm{4.57}\pm \mathrm{8.73}$ and $-\mathrm{9.11}\pm \mathrm{5.16}{\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$, respectively. Green triangles represent winter mean radar freeboard volumes computed with a climatology of radar freeboards between 1995 and 2021, dashed–dot green line is the regression for FBr volume with FBr climatology, and the estimated trend is $-\mathrm{1.80}\pm \mathrm{0.42}{\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$. ${}^{*}(\mathrm{1}-p)<\mathrm{0.5}$, ${}^{**}(\mathrm{1}-p)>\mathrm{0.99}$ and ${}^{***}(\mathrm{1}-p)>\mathrm{0.999999}$ are the probability values of the Mann–Kendall test.

The monthly radar freeboard volume has been finally estimated with Eq. (10) by summing radar freeboard volumes within the P cells, with S_i the grid cell area (12 500² m²) and FB_ri and SIC_i the individual radar freeboard and sea ice concentrations from NSIDC (introduced in Sect. 2.2) for each grid cell. We have decided not to convert estimations to sea ice volume to limit bias coming from snow depth estimates.

The evolution of the snow load is not taken into account in Fig. 13, which means that the evolution of the volume is not fully represented, in the same way as if the total volume were derived with a snow depth climatology. Indeed, a decrease in FBr volume may merely indicate that the snow depth is greater and the ice thickness unchanged.

$$\begin{array}{}\text{(10)}& V_{\mathrm{FBr}}=\sum _{i=\mathrm{1}}^P{S}_{\mathrm{i}}\cdot {\mathrm{FB}}_{\mathrm{ri}}\cdot {\mathrm{SIC}}_{\mathrm{i}}\end{array}$$

Three trend lines fitted on winter mean volumes are represented in Fig. 13: one is computed considering Envisat and CS-2 estimates only (dashed line), a second one from all mission estimations (solid line) and a third one with a FBr climatology from 1995 to 2021 (dashed–dot line). Trends are performed using the Theil–Sen estimator (Theil, 1950; Sen, 1968) and SciPy (Virtanen et al., 2020), providing uncertainties in the regression given with a 95 % confidence interval along with a statistical significance Mann–Kendall test (Mann, 1945; Kendall, 1990). All the trends are decreasing: $-\mathrm{4.57}\pm \mathrm{8.73}$, $-\mathrm{9.11}\pm \mathrm{5.16}$ and $-\mathrm{1.80}\pm \mathrm{0.42}{\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$, respectively. The negative trend is strengthened with the integration of ERS-2 estimates and becomes statistically significant. The uncertainties from regression remain high but are induced by the large interannual variability. The trend in mean winter FBr volume obtained from a FBr climatology between 1995 and 2021 is lower but more reliable, suggesting that FBr variation over the last 30 years dominates FBr volume variations, in contrast to sea ice concentration. In comparison, the Ob River mean discharge is about 400 km³ yr⁻¹ if the radar freeboard volume is roughly converted to total volume with a factor of 10, and the Arctic sea ice decline rate (up to 81.5^∘ N) is about one-quarter of the Ob River mean annual discharge.

5 Conclusions

This study presents a methodology to recover the radar freeboard by first correcting ERS-2 heights from the pulse blurring effect and then adjusting radar freeboards over the Envisat and CryoSat-2 missions. The pulse blurring effect correction is based on the Peacock (1998) and Peacock and Laxon (2004) approaches. The adjustment function developed in this study relies on a multilayer perceptron, which is trained during the common flight period between ERS-2 and Envisat missions using the Envisat radar freeboard as a reference. To ensure consistency along the three altimeters, the Envisat radar freeboard has been preliminarily corrected against CryoSat-2 using the same neural network. The choice of a fixed-threshold retracker (TFMRA) to process the waveforms was motivated by continuity purposes, as it can be used for all radar altimetry missions. The methodology to estimate ERS-2 and Envisat radar freeboards is based on the CryoSat-2 TFRMA50 radar freeboard, which is assumed to be the reference in this study. This hypothesis can be balanced regarding the progress in physical retrackers. The final NN FBr does not conventionally result from a difference in two retracked heights but corresponds to a TFMRA50 SAR-like radar freeboard corrected by a neural network.

The uncertainty estimation is initially tackled, referring to the previous studies of Ricker et al. (2014) and Landy et al. (2020). These uncertainties are then propagated through the neural network thanks to a Monte Carlo approach. Uncertainties in LRM-calibrated radar freeboards range from a few millimeters up to about 15 cm depending on the ice type and the density of the along-track measurements in the grid cell. One of the limitations of the uncertainty budget is the potential underestimation of the impact of radar penetration within the snow layer, which could lead to an underestimation of the radar freeboard uncertainty.

Envisat-corrected radar freeboards show good consistency with CryoSat-2 estimation, with a mean bias of 3 mm for both common winters and a SD of 9.8 cm. The Envisat radar freeboard was then compared to a large sample of validation data: 9 of the 10 datasets give consistent results, especially the strong correlation with the mooring (0.63 to 0.94) and CanCoast stations. Apart from CanCoast, these results are nearly systematically negatively biased, suggesting an underestimation of the radar freeboard. A part of the latter bias is probably due to the draft measurement's method, suggested by the increase in the bias depending on the thickness of the ice. In any case, these biases remain within the range of estimated uncertainties. The result is also consistent with the solution proposed by SI-CCI. Even if the main purpose of this study is to extend the radar freeboard time series to ERS-2, it is nevertheless fundamental to ensure the reliability of Envisat as a reference.

The ERS-2-corrected radar freeboard is close to Envisat-corrected ones, with a correlation of 0.88, a bias of 2 mm and a standard deviation of 3.8 cm for the difference. Comparisons with the few sets of in situ data reveal the same positive bias with CanCoast as for Envisat and high negative biases for submerged draft measurements. Except for AirEM, the comparisons provide consistent correlation values between 0.55 and 0.74. Indeed, these statistics are those expected for comparisons between measurements from different technologies (airborne lasers, ULS moorings, etc.), recording different physical quantities (draft, radar freeboard) with different spatial and temporal availabilities (point, monthly averages). Unless the comparison methodology is reconsidered, it seems difficult to obtain better correlations.

This work finally allows reconstruction of 27 years of Arctic radar freeboards up to 81.5^∘ N and suggests a decline in the sea ice radar freeboard volume of $\mathrm{9.11}\pm \mathrm{5.16}{\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ (about 1 order of magnitude more for the total volume). This decline is significantly greater than considering only the Envisat and CryoSat-2 period. Radar freeboard variations have a predominant influence on volume variations but also on the trend, in contrast to sea ice concentration, which seems to have a moderate impact. In the near future, the methodology will be extended to the ERS-1 mission and for austral sea ice to recover 30 years of sea ice volume variation for both hemispheres. These extended data will also be freely available to the community at large. This radar freeboard time series product based on CryoSat-2 estimations intends to provide a record of monthly sea ice changes over the last 3 decades and for climate studies.

Appendix A: Additional technical content and statistics

Figure A1 Comparison of the standard deviation of ASA between Envisat (cycle 12) and ERS (cycle 80) before (a) and after blurring correction (b) within each grid cell of a 12.5 km resolution grid. The median of the standard deviation difference for panel (a) is 7.5, and for panel (b) it is 0.80 cm.

Table A1 Pulse peakiness thresholds for lead and floe classification.

^* Guerreiro et al. (2017).

Table A2 Radar altimeter characteristics with “Number of echoes” the approximate number of individual echoes that have been summed up to deliver the 20 Hz or 18 Hz waveforms. σ_1b is the estimated error in the range from speckle noise.

^* Wingham et al. (2006).

Table A3 Neural-network-selected hyperparameters; other values were set by default in the MLP regression function from scikit-learn (Pedregosa et al., 2011).

Figure A2 An outline of the link established by the NN between some of the inputs (standardized LES, PP, sea ice concentration and TFMRA50 FBr) and Envisat NN FBr. WPU means the waveform power unit for Envisat correction.

Figure A3 An outline of the link established by the NN between some of the inputs (standardized LES, PP, sea ice concentration and TFMRA50 FBr) and ERS-2 NN FBr. WPU means the waveform power unit for ERS-2 correction.

Table A4 Monthly statistics (average and standard deviation) on the radar freeboard for each winter month from the overlap period between Envisat and CryoSat-2 and averaged for both winters (2010–2012). SD corresponds to the standard deviation, and “unc” stands for uncertainty.

Table A5 Statistics (average and standard deviation) on the radar freeboard for each winter month from the overlap period between ERS-2 and Envisat and averaged for the whole winter (2002–2003). SD corresponds to the standard deviation, and “unc” stands for uncertainty.

Table A6 Statistics of the sea ice thickness difference between Envisat and each validation dataset.

Table A7 Statistics of the sea ice thickness difference between ERS-2 and each validation dataset.

Appendix B: List of abbreviations

Abbreviations	Description
ASA	All-surface anomaly
AGC	Automatic gain control
ASD	Altimetric snow depth
CHK	Chukchi
CS-2	CryoSat-2
EBS	Eastern Beaufort Sea
Env	Envisat
ERS	European Remote-Sensing Satellite
FB	Sea ice freeboard
FBr	Radar freeboard
FBt	Total freeboard
FYI	First-year ice
HTL	Height-tracking loop
ILA	Ice-level anomaly
LRM	Low-resolution mode
MSS	Mean sea surface
MYI	Multiyear ice
MLP	Multilayer perceptron
NN	Neural network
PRF	Pulse-repetition frequency
RA	Radar altimeter
RMSE	Root mean squared error
SAR	Synthetic aperture radar
SARM	Synthetic aperture radar mode
SARin	Synthetic aperture radar interferometric
SD	Standard deviation
SIT	Sea ice thickness
SLA	Sea-level anomaly
STL	Slope-tracking loop
TFMRA	Threshold first-maximum retracker algorithm
ULS	Upward-looking sonar