A New Retracking Algorithm for Retrieving Sea Ice Freeboard from CryoSat-2 Radar Altimeter Data during Winter–Spring Transition

Abstract

A new method called Bézier curve fitting (BCF) for approximating CryoSat-2 (CS-2) SAR-mode waveform is developed to optimize the retrieval of surface elevation of both sea ice and leads for the period of late winter/early spring. We found that the best results are achieved when the retracking points are fixed on positions at which the rise of the fitted Bézier curve reaches 70% of its peak in case of leads, and 50% in case of sea ice. In order to evaluate the proposed retracking algorithm, we compare it to other empirically-based methods currently reported in the literature, namely the threshold first-maximum retracker algorithm (TFMRA) and the European Space Agency (ESA) CS-2 in-depth Level-2 algorithm (L2I). The results of the retracking procedure for the different algorithms are validated using data of the Operation Ice Bridge (OIB) airborne mission. For two OIB campaign periods in March 2015 and April 2016, the mean absolute differences between freeboard values retrieved from CS-2 and OIB data were 9.22 and 7.79 cm when using the BCF method, 10.41 cm and 8.16 cm for TFMRA, and 10.01 cm and 8.42 cm for L2I. This suggests that the sea ice freeboard data can be obtained with a higher accuracy when using the proposed BCF method instead of the TFMRA or the CS-2 L2I algorithm.

1. Introduction

Sea ice thickness is an important environmental parameter, which is related to the ice mass balance and which controls the local and regional heat budget. Thus, it affects also the global climate system. Due to limitations in spatial and temporal coverage of in situ ice thickness data, there is a demand on improving the retrieval of ice thickness from satellite data [1]. Various methods have been proposed, using data from passive microwave sensors (e.g., special Sensor Microwave/Imager [2,3], Advanced Very High Resolution Radiometer [4], and Advanced Microwave Scanning Radiometer-Earth Observing System [5]), optical sensors (e.g. Moderate Resolution Imaging Spectroradiometer [6]), and synthetic aperture radar [7]). In recent years, numerous satellite altimeter missions have been launched that have provided data for retrieving sea ice thickness [8,9,10,11,12,13,14]. Data from spaceborne altimeter instruments are important to cover the entire polar regions over longer time periods and to provide information on the spatial, annual, and interannual variability of sea ice thickness [15]. In addition, altimeters can be also used to support large-scale sea ice charting [16]. Sea ice thickness data were or are obtained from data of the European Remote Sensing (ERS) satellites, ENVISAT, ICESat, Sentinel-3, and CryoSat-2. CryoSat-2 (CS-2), which was launched in April 2010, is equipped with a Ku-band Synthetic Aperture Interferometric Radar Altimeter (SIRAL) and offers much higher spatial resolution (~360 m along track) [17] and larger global coverage (up to 88°N) than earlier space-borne altimetry missions. It is one of the most advanced satellite altimeters for sea ice observation, providing more detailed information about spatial and temporal ice thickness variations than hitherto available.

Acquisitions of sea ice freeboard data is the prerequisite for calculating sea ice thickness from altimeter data and hence has a significant effect on the accuracy of the ice thickness calculations. Ice freeboard is obtained from subtracting the local sea surface elevation from ice elevation. The elevation of a surface is determined from the so-called tracking point in the received radar pulse. In the following, the received radar pulse is denoted as the ‘waveform’. For sea ice and water, this is the point on the leading edge of the echo that corresponds to the mean level of scattering from the surface [18]. Ideally, the tracking point is set at the center of a predefined range window in which the waveform reflected from the surface is expected to appear. However, due to varying surface topographies and the variation of the penetration depths of radar signal into the snow/ice volume, the actual tracking point may reveal an offset relative to the predefined tracking point. In order to obtain accurate surface elevation data, it is necessary to consider this offset between the center of the range window and the real tracking point [8]. This process is called retracking correction. The performance of the retracking correction can introduce random noise into the retracking correction (i.e., affecting precision) and has a direct and essential impact on the accuracy of the ice freeboard retrieval. Imperfectly fitted waveforms can also introduce random noise into the retracking correction, i.e., affect precision.

The common retracking methods used for sea ice can be divided into three classes: the threshold first-maximum retracker algorithms (TFMRA), the physical model (PM) approach, and non-physical (i.e., purely mathematical) waveform fitting (WfF). For TFMRA, the retracking point is fixed to the position at which the received power reaches a certain level relative to the first peak of the waveform [9,11]. Commonly used thresholds are 40% or 50% [1,11]. Ricker et al. [19] oversampled the waveform and applied a boxcar averaging procedure to generate a smoothed CS-2 echo, with retracking points for ice and leads at the 40%, 50%, and 80% level of the first local power maximum, respectively, which are values that have been used in recent studies. They found that the 40% and 50% threshold tracked above the ice surface, whereas the 80% threshold tracked below the ice surface compared with airborne laser altimetry data. With this method, problems occur if snow and ice conditions (snow depth, density, brine content, sea ice salinity, and sea ice surface roughness) in the investigated area are highly variable. This means that the positions of thresholds would have to be adjusted accordingly on spatial scales corresponding to the scales of snow and ice properties variations, which also may change as a function of time. This step, however, cannot be carried out with TFMRA and WfF retrackers due to the lack of information required for such adjustments. In the PM approach, the retracking point position is adjusted based on the shape of the waveform, which in turn is related to sea ice properties [20].

In applications of the PM approach, the shape of the radar waveform is simulated considering technical specifications of the altimeter system and different sea ice parameters and then compared to the measured radar echo. This step includes to find the echo delay time for which the simulated waveform deviates least from the observed one, which in turn determines the retracking point position [12,21,22,23]. Kurtz et al. developed a physically-based approximation of the waveform which takes into account the radar backscatter coefficient, incidence angle, and ice surface roughness [12]. Relative to Operation Ice Bridge (OIB) data acquired on three airborne campaigns from March 2011 to March 2013, taken as reference, they found RMSE values ranging from 7.4 to 11.1 cm when applying their PM approach, and 14.1–19.8 cm for a TFMRA approach with a threshold of 50% for sea ice. Note that in [12], the PM approach is denoted as CS2 WfF. However, the difference between CS2 WfF and TFMRA is not valid in general as the results of [6] indicate.

Poisson et al. proposed a modified Brown model [24] coupled with a surface identification method, and used an adaptive maximum-likelihood retracker [25]. Similarly, a Brown-Hayne model that well matches the radar signal received from ocean surfaces was applied for adaptive retracking over leads by Passaro et al. [26]. Landy et al. [20] developed a numerical facet model of statistical surface topographies to simulate the heterogeneous sea ice radar echoes. The results show that this model can improve retrievals of key sea ice properties (including freeboard, surface roughness, and snow depth). Since the PM approach needs to consider explicitly ice properties and interactions between the radar signal and the ice, it is mathematically more complex.

For WfF, an appropriate mathematical function is selected to fit the echo shape, and then a threshold for the retracking point is determined empirically [10,14]. For example, Giles et al. used a Gaussian function for the leading edge of the waveform, followed by a linear function and a negative exponential for the tail as retracker to retrack leads [13]. The retracking point is positioned at the breakpoint between the Gaussian and the linear function. Yi et al. used three different retracking methods (75% threshold method, Gaussian peak method, and the centroid method) to estimate elevation [14]. The results indicate that the Gaussian and centroid methods yield more accurate surface elevation data than the application of a fixed threshold.

The WfF approaches have become popular in recent years [10,13,14]. One reason is that they are easier to implement compared to PM approaches. On the other hand, unlike TFMRA, WfF algorithms can be adapted to variations of echo waveform changes. Hence in this study, we focus on a new algorithm belonging to the group of WfF algorithms. Because the radar echo is sensitive to the properties of snow and ice, the shape of the waveform tends to be complicated and highly variable, which means that a predefined fitting function usually does not result in an optimal fit, and the estimated elevation may reveal a lower accuracy than potentially achievable. To tackle this problem, we propose a new retracking correction based on Bézier curve fitting (in the following sections denoted as BCF) applied to the CS-2 echo waveform. The Bézier curve has been widely used to estimate complex curves in computer graphics due to its good curve fitting performance [27]. In comparison with airborne radar altimeter data from the OIB mission, we show that the proposed method can deliver highly accurate ice freeboard data.

The paper is organized as follows: Section 2 describes the used data sets and Section 3 introduces the procedure to fit CS-2 waveform using Bézier curve. A procedure for the retrieval of ice freeboard and a retracking correction based on Bézier curve fitting are both introduced in Section 4. An evaluation and comparison of the proposed method is given in Section 5. The discussion and conclusions are provided in Section 6 and Section 7.

2. Data Sets

2.1. CryoSat-2 Data

The data used in this study are Baseline C CS-2 data from March 2014, March 2015 and April 2016. The spatial resolution of the CS-2 measurement is approximately 0.3 km by 1.5 km in along and across track direction, respectively [17]. CS-2 has a full 369-day cycle and a subcycle of 30 days, which ensure that the whole Arctic is covered in one month. CS-2 is operated in three modes: “Low Resolution” (LRM) for measurements over the open water, “SAR” for sea ice, and “SARIn” for the ice sheets of Greenland and Antarctica. Here only CS-2 SAR-mode data was used.

During the CS-2 commission phase, it became clear that the SAR mode data over specular surfaces, like sea ice leads, were aliased. The aliasing effect was minor over diffuse surfaces (sea ice floes). This has been corrected by oversampling the CS-2 SAR echoes in data series since 2012 [17]. The aliasing effect and its correction are discussed in [23] and [28]. The CS-2 Level-1b (L1b) data used here are based on the oversampled radar echoes, and a detailed processor can be seen in [19].

The L1b data provide the two-way travel time of the radar echoes, which is used to calculate the range between the satellite’s center of mass and the reflecting horizon on the ground. The processing of the L1b data includes instrumental corrections, multi-looking, corrections of position errors, noise suppression, and geophysical corrections. The latter considers the wet and dry tropospheric delay time; propagation through the ionosphere; atmospheric pressure loading; and ocean, solid earth, and pole tides [17]. We applied the correction values that are provided in the metadata file, to the L1b data according to [17]. The data quality flag in CS-2 L1b data is also used to filter poor-quality waveform data which ensure the used data are reliable [17]. In the SAR mode, the echoes are sampled in 256 range bins, each covering a time interval of 1.563 ns (0.234 m). Detailed specifications of CS-2 SIRAL are presented in

Table 1. Specifications of the CS-2 data.

	LRM	SAR	SARIn
Center frequency	13.575 GHz
Bandwidth	350 MHz
Pulse Repetition Frequency	1.97 kHz	17.8 kHz	17.8 kHz
Range window	~60 m	~60 m	~240 m
Samples per echo	128	256	1024
Sample interval	0.4684 m	0.2342 m	0.2342 m

.

For comparison, we used CS-2 In-depth Level-2 (L2I) data, which contain surface elevation and sea ice freeboard information, to evaluate the performance of the proposed method. The spatial and temporal resolutions of CS-2 L2I correspond to the L1b data (approximately 0.3 km by 1.5 km in along and across track direction). All necessary corrections, separation of waveforms in lead and ice floes, and sea surface height interpolations, are identically applied in all used method (BCF, CS-2 L2I and TFMRA, the latter two will be introduced in Section 5).

2.2. Operational IceBridge Data

In this study, we compare sea ice freeboard data acquired as part of NASA’s OIB airborne mission with the freeboard retrieved from nearly coincident CS-2 data acquired in the Arctic in March 2015 and April 2016. The OIB mission, initiated by NASA, utilizes a laser altimeter for large-scale surveys of Arctic sea ice on an annual basis. The primary goal of OIB, which was started in 2009, was to generate a laser altimetry time series bridging the gap between the ICESat missions [29]. Here, the OIB sea ice freeboard, snow depth, and thickness quick look data from the periods of 20–27 March 2014, 24–30 March 2015, and 20–21 April 2016 are used. The time differences between the OIB mission and CS-2 satellite acquisitions were less than 3 hours.

The OIB mission is based on three instruments: an airborne topographic mapper (ATM) laser altimeter [30], a digital mapping system (DMS) camera [31], and a snow radar [32]. The ATM is a 532-nm wavelength scanning laser altimeter for retrieving the surface height. The footprint size of each individual elevation measurement is determined by the laser beam divergence and flight altitude. The footprint of the measurements used in this study is about 1 m [29]. The position accuracy is better than 1 m and the overall elevation accuracy from the ATM is about 0.1 m [33,34]. The DMS and snow radar are used to detect ice types and snow depth, respectively. Snow depth is derived from reflections at the air–snow and snow–ice interfaces [29,35]. The vertical resolution of the snow radar is about 0.06 m, and the nominal spatial resolution is 40 m. Freeboard, thickness, and snow depth from all three instruments were resampled at a 40 m spatial resolution [29], and uncertainty estimates are also provided with the data products. In this study, we used ice freeboard for which the uncertainty is less than 0.1 m and snow depth data from the OIB Quick Look products. Since the ATM laser pulses are reflected at the air–snow interface, the snow depth must be subtracted from the retrieved ATM freeboard before comparing it to the CS-2 radar freeboard. Figure 1 shows the ice freeboard obtained from the OIB Quick Look products in March 2015 and April 2016 on a map of the Arctic with a 25 km grid.

2.3. Auxiliary Data

2.3.1. Sea Ice Type Information

We used daily data provided by the Ocean and Sea Ice Satellite Application Facility (OSI-SAF) to separate first-year ice (FYI) and multi-year ice (MYI). The sea ice classification is based on data from the Special Sensor Microwave/Imager (SSM/I) and Advanced SCATerometer (ASCAT) [36]. The temporal variation of MYI area coverage is used for a quality assessment of OSI SAF ice-type data as the area of multi-year ice is assumed to change only slightly over a month. Aaboe et al. showed that monthly standard deviations of the daily variation of the MYI extent are lower than 100 km² [36].

2.3.2. Snow Depth and Mean Sea Surface Height Data

To retrieve ice freeboard over the whole Arctic region, snow depth data from the University College London climatology model (UCL04) [17,37] are used. UCL04 is a hybrid mean sea surface model, which provides monthly snow depth data for the Arctic with spatial resolution of 0.0625 degrees. [17]. Comparing to OIB data, snow depths from UCL04 model are larger by about 0–0.1 m, and this difference can be particularly observed over areas of FYI [11], which we considered in our analysis (see discussion). Together with the ULC04 data, the regionally and temporally limited but highly accurate OIB snow depths serve for investigations of the uncertainties of ice freeboard retrievals caused by the snow cover. Furthermore, we employed DTU15 mean sea surface height (MSS) data [38,39]. An assessment of its accuracy can be found in [40,41].

3. CS-2 Waveform Fitting Based on Bézier Curves

To consider the observed shape variations of the echo waveforms from different surface types (ice and water), we apply a fitting approach based on composite Bézier curves. The segmentation strategy proposed below ensures the high-quality and high-efficiency of the fitting procedure to the CS-2 waveforms. The fitted waveforms can be used to calculate different waveform parameters but, in this paper, we focus on the retrieval of the retracking point.

3.1. Bézier Curves

A Bézier curve is a parametric curve frequently used to fit functions with complex shapes. It has been widely applied in computer graphics (e.g., animations), user interface design, and smoothing cursor trajectories [42]. An n-order Bézier curve fit is calculated from a set of control points p₀ to p_n [43]. The formula can be expressed as

$$B\left(t\right)={\displaystyle \sum _{i=0}^n{b}_{i,n}\left(t\right)p_i},0\le t\le 1$$

(1)

where $b_{i,n}\left(t\right)=C_n^i{(1-t)}^{n-i}{t}^i$, which are known as the Bernstein basis polynomials of degree n, $C_n^i$ are the binomial coefficients, p_i denotes the control points. The independent variable t is restricted to the interval between 0 and 1, t∈(0, 1).

Since the shape of the n-order Bézier curve is sensitive to the positions of the control points, their determination is critical to curve fitting [44]. In general, the positions of control points can be derived from the sample points in the original curve. For example, if an n-order Bézier curve is chosen to model the surface waveform, m+1 sample points q(t₀), q(t₁),…, q(t_m) from the measured waveform can be used to solve for the control points by applying the equation

$$\mathsf{\Phi }P=Q$$

(2)

Here, $\mathsf{\Phi }$ is a Bernstein matrix with (m+1) × (n+1) dimensions

$$\mathsf{\Phi }=\left(\begin{array}{cccc}{b}_{0,n,0}& b_{1,n,0}& \dots & b_{n,n,0}\\ b_{0,n,1}& b_{1,n,1}& \dots & b_{n,n,1}\\ ⋮& ⋮& \ddots & ⋮\\ b_{0,n,m}& b_{1,n,m}& \dots & b_{n,n,m}\end{array}\right)$$

(3)

$P={[p_0,p_1,\dots ,p_n]}^T$ is the n+1 dimensional control point vector and $Q={[q(t_0),q(t_1),\dots ,q(t_m)]}^T$ is the m+1 dimension sample point vector.

The number of sample points m+1 is usually larger than the order n of the Bézier curve, which will lead to an overdetermined system of equations. Hence, Equation (2) always has a solution. The least squares method described in [45] is applied to find an optimal solution

$$\min \left|\right|\mathsf{\Phi }P-Q\left|\right|$$

(4)

The solution to Equation (2) can be written using a normal equation

$$P={({\mathsf{\Phi }}^T\mathsf{\Phi })}^{-1}{\mathsf{\Phi }}^{T}Q$$

(5)

3.2. Fitting Strategy

The altimeter echo waveform can be closely approximated by a single Bézier curve with sufficiently high order n. However, Lang et al. (2016) pointed out that the application of composite Bézier curves has a higher computational efficiency and stability [44]. Therefore, we choose to fit composite cubic Bézier curves to the CS-2 altimeter echo waveforms. To this end each waveform is divided into several segments. The equation for the cubic Bézier curve equation can be written as

$$B(t)={(1-t)}^3{p}_0+3t{(1-t)}^2{p}_1+3t^2(1-t)p_2+t^3{p}_3$$

(6)

In our implementation, the Bézier curve is set to pass through its first and last control point in each segment, i.e., q(t₀) = p₀ and q(t_m) = p₃ using the notation of Equation (2). To simplify calculations, the optimal position of the control points p_i in each segment is determined by minimizing the squared differences between sample points and fitted data

$$\min {\displaystyle \sum _{k=0}^m{[q(t_k)-B(t_k)]}^2}$$

(7)

The performance of the composite cubic Bézier curve fitting is directly dependent on the division of the waveform into different segments. Three common segmentation strategies are analyzed: strategy a: segments with exponentially distributed breakpoints, the fixed breakpoints in terms of the variable t being (0,1/16,1/8,1/4,1/2,1); strategy b: segments with a set of uniformly distributed breakpoints (0, 1/5, 2/5, 3/5, 4/5, 1); and strategy c: fixed segments with a double sampling density compared to strategy b. Figure 2 shows examples of typical CS-2 waveforms for lead and sea ice using the three segmentation strategies. For the leading edge, strategy c results in much better fitting results than strategies a and b. We note that the absence of enough breakpoints over the leading edge of the waveform causes bad matches in case of all three strategies. Hence, we propose another strategy for segmentation: the waveform is divided into five segments with two breakpoints on the leading edge and two breakpoints on the trailing edge. The first and the last breakpoint are assigned to the positions where the signal rises and drops to 5% of the maximum of the first peak, and the other two breakpoints are located at the nearest bins in front and behind the first peak. With this strategy, we obtain good fitting results for the raising and trailing edge and around the waveform peak.

Figure 2 demonstrates the effectiveness of our proposed segmentation strategy (in the figure denoted as ‘strategy d’) when comparing it with the results of strategies a, b, and c. It shows that the new segmentation strategy matches the radar waveforms very well with less segment points, except for short-scale oscillations in the trailing edge of the ice waveform. The quality of the fit on the leading edges enables a good retracking correction.

To quantitatively evaluate the fitting performance, we applied the root mean square error (RMSE) between fitted and original waveform. The RMSE in terms of waveform power and range bin are separately shown in

Table 2. Performance of different segmentation strategies for waveform fitting (Parm: parameters for comparison; le: leading edge; w: whole waveform). The used CS-2 data were collected in March 2014. Further explanations are given in the text. The RMSE in terms of waveform power and range bin are separately shown here.

	Parm	(a)	(b)	(c)	(d)	Parm	(a)	(b)	(c)	(d)
Leads (le)	RMSE 1	0.0822	0.0902	0.0097	0.0072	RMSE 2	3.13	4.77	0.38	0.27
Leads (w)	RMSE 1	0.1298	0.1335	0.0437	0.0128	RMSE 2	4.37	4.47	1.31	0.39
	ACT 3 (s)	0.1712	0.1730	0.1851	0.1510	ACT(s)	0.1712	0.1730	0.1851	0.1510
Sea ice (le)	RMSE 1	0.0566	0.0447	0.0298	0.0078	RMSE2	2.57	4.12	2.55	1.62
Sea ice (w)	RMSE 1	0.0776	0.0589	0.0438	0.0335	RMSE2	2.19	4.02	2.26	1.51
	ACT(s)	0.1816	0.1855	0.1901	0.1715	ACT(s)	0.1816	0.1855	0.1901	0.1715

¹ RMSE in terms of waveform power, ² RMSE in terms of range bin, ³ ACT: average computation time.

based on all valid CS-2 waveforms (~1,400,000) obtained during measurements in March 2014. Both results demonstrate that the segmentation strategy d performs best in fitting composite Bézier curves to both lead and ice floe waveforms. This is valid both for considering only the leading edge and the whole waveform. The average computation time of the proposed method is about 0.1625 second per waveform, with MATLAB R2015b 64-bit software on a personal computer having an i7 Quad Core (3.40 GHz) processor with 4 GB of RAM.

4. Sea Ice Freeboard Retrieval

In this section, sea ice freeboard is retrieved using the following methodology: the first step is to distinguish open water, leads, and sea ice. In the next step, a composite Bézier curve is fitted to the CS-2 waveform to determine the retracking point (hereinafter, we refer to this step as ‘BCF’). Following this step, the surface elevation is calculated. At last, the freeboard is obtained from the difference between the sea ice elevation and the local sea surface level. In Figure 3a flowchart of the scheme for retrieving the ice freeboard is depicted.

4.1. Discrimination of Surface Types

Different waveform shapes can be assigned to different surface types. The detection of leads relies on the fact that the radar echo is dominated by specular reflections [8,46]. This is the case even if the width of the lead is smaller than the altimeter footprint [46]. Our way to detect leads is based on the pulse peakiness (PP) and the stack standard deviation (SSD). The latter uses stacks of different beam echoes and thus provides a measure of the variation in surface backscatter with incidence angle [37]. PP is calculated by

$$PP=\frac{\max (wp(i))}{{\displaystyle \sum _{i=1}^{256}wp(i)}}\cdot 256$$

(8)

where wp(i) represents the normalized echo waveform power in the range bin with index i. PP is sensitive to the shape of the echo waveform. It can discriminate between specular reflections and diffuse scattering with reasonable accuracy. The SSD provides a measure of the variation in surface backscatter [37]. Returns from leads are identified by their high PP (>40) and low SSD (<4) [19]. In addition to analyzing the altimeter waveform, we employed daily OSI SAF data to separate FYI, MYI, and open water. As the OSI SAF data are nearly coincident with the CS-2 data the uncertainty caused by the sea ice drift can be ignored.

4.2. Waveform Retracking

Here, we describe how the optimal power threshold for setting the retracking point in waveforms approximated by composite cubic Bézier curves was determined. For testing, we assigned retracking points to the range bins for which the power of the leading edge was at 50, 70, 90, and 100% relative to the first maximum of the Bézier curve, covering the commonly used thresholds in the current literature. Ice freeboard was retrieved from CS-2 waveforms using 16 different pairs of thresholds (in each pair the first for leads, the second for ice floes). Because the backscattering from leads is dominated by specular reflection and from ice the backscattering is more diffuse, it is unphysical to assume that the threshold for leads could be lower than for ice floes. Hence, all pairs with lead threshold lower than ice threshold are not used in this study. The test dataset is from March 2014. The results of the retrievals are compared with the corresponding OIB data. The latter are averaged into 25 × 25 km grid cells, resulting in totally 595 grid cells which contain both CS-2 and OIB freeboards (the latter is OIB snow freeboard – OIB snow depth). In the comparison, about 1,400,000 individual waveforms are included, and the covered area is about 371,875 km².

Table 3. Retrieved sea ice freeboard using BCF with different retracking thresholds, compared with the freeboard determined from OIB March 2014 data for overall, FYI, and MYI separately. Unit: cm.

Thresholds (Leads, Ice)	50%, 50%	70%, 50%	70%, 70%	90%, 50%	90%, 70%	90%, 90%	100%, 50%	100%, 70%	100%, 90%
MD 1	11.73	11.47	36.30	14.49	28.90	65.98	21.89	22.45	55.79
RMSE	15.41	14.05	43.50	17.16	37.22	73.65	25.16	30.72	64.83
Correlation	0.4223	0.4295	0.1596	0.4074	0.1507	0.0292	0.3540	0.1298	0.0156
FYI
MD 1	10.11	9.72	27.21	14.29	20.32	49.70	24.23	17.34	38.80
RMSE	14.74	11.58	33.04	16.71	27.05	56.35	27.03	22.69	48.06
Correlation	0.1247	0.0488	-0.0259	0.0433	−0.0266	−0.0239	0.0085	−0.0413	−0.0350
MYI
MD 1	12.29	11.82	38.23	14.47	30.78	69.57	21.27	23.61	59.62
RMSE	16.09	14.50	45.38	17.19	39.04	76.91	24.61	32.19	67.99
Correlation	0.4348	0.4384	0.2117	0.4203	0.2052	0.1099	0.3802	0.1925	0.1024

¹ MD (Mean freeboard difference): |CS-2−OIB|.

lists the mean absolute difference (MD) and RMSE of ice freeboard retrievals for the different threshold combinations. When calculating the freeboard, we considered the lower propagation speed of the radar signal in snow (the detailed processing is described in Section 4.3). Table 3 reveals that the optimal thresholds are 70% for leads and 50% for ice. For this combination we obtained 11.47 cm for the MD, 14.05 cm for the RMSE and the correlation coefficient is 0.43. The selected threshold combination performed also well for the data from March 2015 and April 2016, see Section 6.1. In addition, in order to investigate the effect of ice surface conditions on the retracking point estimation, we also compared the different retracking thresholds based on FYI and MYI separately. The results show that the same thresholds (70% for leads and 50% for ice) are optimal for both FYI and MYI (Table 3). Hence this combination is selected as threshold for BCF. We remark that in order to include as many data point as possible, the freeboard filter described in Section 4.3 was not applied here. Unreasonable thresholds (e.g., 100% for leads and 90% for ice) can cause a significant amount of unrealistic negative freeboard values.

4.3. Ice Freeboard Calculation

The elevation of each measurement point above the WGS84 ellipsoid is calculated by

$$H=H_S-R-H_R-H_G$$

(9)

where H represents the elevation of the main scattering interface and H_S the elevation of the satellite, both relative to the WGS84 ellipsoid, R is the range, and H_R and H_G represent retracking and geophysical corrections, respectively. R is calculated by multiplying the speed of light in a vacuum with the one-way travel time T of the radar echo. The values of T, H_S, and H_G are taken from the CS-2 L1b metadata file. H_R is calculated using BCF.

For the determination of the sea level we started with a linear interpolation connecting the retrieved spots of open-water or thin ice in leads. We then applied a low-pass filter by using a running mean with 25 km width to smooth jumps that occur in dense lead clusters due to signal noise [19]. This procedure was carried out for each CS-2 track, yielding the sea surface anomaly (SSA, given relative to the mean sea surface height) at the time of data acquisition. In the next step, the MSS and SSA were subtracted from retrieved surface elevation data H, which were identified as sea ice in the surface-type discrimination [19]. The radar freeboard (F_R), in which the lower wave propagation speed in the snow layer is not considered yet, is then obtained from

$$F_R=H-MSS-SSA$$

(10)

where the MSS values were received from the DTU15 model. In our analysis, only radar freeboard values within the interval −0.1 m < F_R < 2.1 m (2 m ± 0.1 m) were considered as realistic [19]. The value of 0.1 m represents the random uncertainty of the range measurement due to speckle [37]. In our data analysis, we have followed the approach used in [19]. Some studies allow larger negative freeboards [47]. For example, in March 2015 using BCF and TFMRA over 10% (for CS2 L2I less than half of it) of all the calculated freeboards were negative. The treatment of these freeboards affects the results.

As the wave propagation speed in the snow layer is lower than in the air, the ice elevation will be underestimated. Hence a correction term must be added to F_R that depends on the snow depth on sea ice [12]. The sea ice freeboard (F, unit: meter) is finally calculated by

$$F_{}=F_R+h_c$$

(11)

Here, h_c (m)is the correction factor given by [12]

$$h_c=h_s(1-\frac{c_{snow}}{c})$$

(12)

where h_s is snow depth (cm) which is usually acquired from an external data source (e.g., the meteorological model, remote sensing observation, or in situ data). In this study, the snow depth data are from OIB and UCL04. The speed of light in the snow, C_snow (m/s), is calculated according to [12]

$$c_{snow}=\frac{c}{\sqrt{1+1.7{\rho }_s+0.7{\rho }_s^2}}$$

(13)

where ρ_s is the density of snow, set to ρ_s = 320 g/cm³. This value agrees with Warren’s estimate for snow density in March and April [48] and has been widely used for freeboard calculation in other studies.

However, it is difficult to consider the effect of snow load on the ice freeboard retrieval. Usually, it is assumed that the radar signal can penetrate through the snowpack to the ice surface when snow is cold and dry, but in several studies, it was found that snow moisture content, brine content, and snow layers of large density reduce the radar penetration [49,50,51,52,53,54]. We will discuss this issue further in Section 6.2. Note that in our analysis we do not consider penetration into the ice volume, which is justified since the strongest reflection is from the ice surface or from vertically higher scattering horizon in the snow layer.

5. Results

In this section, we compare the freeboard results using BCF to those using TFMRA and CS-2 L2I products from March 2015 and April 2016, and also to the airborne measurements from NASA’s OIB campaign. In Section 4.2, we determined optimal thresholds for BCF through comparison with OIB data from March 2014. For further analysis, it is necessary to repeat this procedure for TFMRA, since its performance also depends on the choice of thresholds. We therefore re-programmed the TFMRA retracking algorithm according to [19], including waveform oversampling, noise reduction, and retracking correction. Final results are provided in

Table 4. Sea ice freeboard results from application of TFMRA with different retracking thresholds relative to freeboard data from OIB March 2014 data. Unit: cm.

Thresholds (Leads, Ice)	50%, 50%	70%, 50%	70%, 70%	90%, 50%	90%, 70%	90%, 90%	100%, 50%	100%, 70%	100%, 90%
MD 1	15.21	13.63	35.54	16.24	26.96	66.55	20.05	23.64	60.44
RMSE	23.75	21.06	46.71	21.73	39.90	79.39	24.89	36.67	74.50
Correlation	0.2879	0.2877	0.0853	0.2838	0.0821	–0.0304	0.2622	0.0710	−0.0373

¹ MD (Mean freeboard difference): |CS-2−OIB|.

. As for the results listed in Table 3, negative freeboards were not filtered. We found that a 70% threshold for leads and 50% threshold for ice are optimal, in agreement with [55].

For studying the differences between the results of different retrieval methods, the three methods (BCF, TFMRA, CS-2 L2I) were compared with data collected during the OIB mission in March 2015 and April 2016 (hence assuming that the thresholds determined with March 2014 data remain optimal for the data acquisitions in 2015 and 2016—this item is addressed in the discussion below). All OIB freeboard data are presented in the same 25 km polar stereographic grid, in 374 grid cells for March 2015 and 184 cells for April 2016. In

Table 5. CS-2 ice freeboard retrievals from the BCF, TFMRA, and CS-2 L2I products compared to OIB data in March 2015 and April 2016. Unit: cm.

	March 2015				April 2016
Number of grid points	374				184
Retracker	OIB	BCF	TFMRA	L2I	OIB	BCF	TFMRA	L2I
Mean ice freeboard	21.46 ± 8.59	27.62 ± 14.60	29.83 ± 14.62	30.24 ± 9.55	23.52 ± 10.28	27.15 ± 8.68	22.53 ± 11.07	24.07 ± 11.82
MD 1		9.22	10.41	10.01		7.79	8.16	8.42
RMSE		15.14	15.99	12.42		10.42	10.68	11.03
Correlation		0.3705	0.3917	0.5152		0.4598	0.5079	0.5059

¹ MD (Mean absolute freeboard difference): |CS-2−OIB|.

, we summarize the comparison between BCF, TFMRA, CS-2 L2I products and OIB data. The values listed in Table 5 reveal that the retrieval of sea ice freeboard based on the BCF retracking method provides the best match to the OIB elevation profiles for the used data sets.

Figure 4 shows the probability density function (PDF) of ice freeboard corresponding to the results of BCF, TFMRA, and L2I, without considering regional variations. Overall, the BCF has a minimum mean bias relative to OIB compared with the other methods. For March 2015 and April 2016, the OIB freeboard values are distributed mainly in the interval between 0 m and 0.6 m (Figure 4), and BCF, TFMRA, and L2I tend to overestimate the ice freeboard, which can also be seen in Table 5. The spatial resolution difference between OIB and CS2 data is one reason for deviations between freeboard values. The smallest freeboard estimates can be found for BCF and TFMRA. The slope of the leading edge of the waveform is a function of surface roughness. It decreases with increasing roughness. If the roughness is, e.g., Gaussian distributed around the mean surface elevation, and the 50% waveform power point in the leading edge is kept fixed at the same time of signal arrival independent of roughness, there will be a shift for all other thresholds. Since we use thresholds ≥50%, the arrival times we obtain are shifted towards larger values, which correspond to smaller freeboards. With a fixed retracker we should hence often get too small freeboards over rough ice. In Figure 4, the freeboard differences reveal some deviations between the three methods. The differences are small in April 2016. For March 2015, the mean freeboard difference with respect to OIB is 0.0598 ± 0.1393 m for BCF, 0.0819 ± 0.1375 m for TFMRA, and 0.0861 ± 0.0897 m for L2I, respectively. The corresponding mean freeboard difference values for April 2016 are 0.0012 ± 0.1106 m for BCF, −0.0142 ± 0.1061 m for TFMRA and 0.0320 ± 0.0995 m for L2I, respectively. According to Table 5, the MD (RMSE) values for BCF are 0.0922 m (0.1514 m) in March 2015 and 0.0779 m (0.1042 m) in April 2016. The corresponding numbers are 0.1041 m (0.1599 m) and 0.0816 m (0.1068 m) for TFMRA and 0.1001 m (0.1242 m) and 0.0842 m (0.1103 m) for L2I. Hence, compared to TFMRA and L2I, BCF gives the best results for this data set in both test months.

The sea ice freeboard retrieval errors are further investigated for the FYI and MYI regions respectively (

Table 6. The ice freeboard retrievals from the BCF, TFMRA, and CS-2 L2I products are compared with OIB data for first-year ice (FYI) and multi-year ice (MYI) in March 2015 and April 2016. Unit: cm. Snow data were obtained from the airborne OIB mission.

	March 2015				April 2016
Retracker	OIB	BCF	TFMRA	L2I	OIB	BCF	TFMRA	L2I
Number of grid points in FYI	59				45
Number of grid points in MYI	315				139
OIB Snow Depth in FYI	15.03 ± 7.01				8.39 ± 5.34
OIB Snow Depth in MYI	23.37 ± 9.67				24.87 ± 11.17
Mean FYI freeboard	19.39 ± 5.92	19.17 ± 8.21	21.43 ± 8.32	24.11 ± 9.45	17.81 ± 18.17	22.02 ± 4.94	14.36 ± 5.49	15.47 ± 5.86
MD 1 for FYI		6.70	6.88	7.89		13.42	12.02	12.60
RMSE for FYI		9.03	9.27	10.55		18.59	18.65	18.60
Correlation		0.2015	0.2135	0.2460		0.1063	0.0844	0.07539
Mean MYI freeboard	21.90 ± 9.05	28.78 ± 14.86	30.97 ± 14.87	31.14 ± 9.44	25.81 ± 10.27	28.68 ± 9.00	24.99 ± 11.21	26.67 ± 12.00
MD for MYI		9.67	11.00	10.48		7.74	8.59	8.87
RMSE for MYI		15.80	16.71	12.83		10.50	10.99	11.39
Correlation		0.3701	0.3926	0.5340		0.4522	0.4787	0.4851

¹ MD (Mean freeboard difference): |CS-2−OIB|.

). The ice types are separated based on the OSI SAF ice type data, see Section 2.3.1. Based on the procedure described in Section 4.2, we determined the optimal threshold for ice without making a difference between MYI and FYI. The results listed in Table 6 show that the mean ice freeboard difference based on BCF is lower than the results based on CS-2 L2I and TFMRA in 2015, and for MYI in 2016. For FYI in April 2016, BCF reveals a slightly lower accuracy than TFMRA and CS-2 L2I. The reason for that is not clear due to the lack of detailed information about snow loading and structure. For BCF, CS-2 L2I, and TFMRA, the accuracy of the estimated FYI freeboard is higher than for the MYI freeboard in March 2015, while in April 2016 the situation is opposite. Overall, the BCF freeboard shows an overestimation except for FYI in March 2015. TFMRA overestimates the ice freeboard in March 2015 and underestimates it in April 2016. CS-2 L2I overestimates the ice freeboard except for FYI in April 2016.

6. Discussion

The selected positions of retracking points and the snow loading on sea ice affect the accuracy of ice freeboard estimation. We examine these issues in the following sub-section.

6.1. Thresholds of the BCF Method

The undersampling of the waveform signal is common for radar altimeter echoes. Especially the peaky and narrow echoes of specular surfaces such as leads are undersampled, which influences the results of retracking [28]. Since undersampling only reveals average values for each range bin, the real maximum value of the echo is missed. Aliasing has to be considered also for CS-2 even though L1b data are now oversampled as noted in Section 2.1. However, the sampling frequency of CS-2 waveform is four times larger than for the ERS-2 altimeter echo and two times higher than for the echo of the ENVISAT altimeter. Therefore, the application of a threshold retracker is justifiable for elevation estimates from CS-2 data. In particular on the steep slope of the leading edge, the Bézier curve can provide a close nonlinear interpolation between adjacent sampling points. The quality of the interpolation and the resulting retracking position depends on the number and position of breakpoints, see Figure 2.

In order to quantitatively analyze how the location of breakpoint affects the calculated retracking point, we artificially changed the location of the breakpoints on the leading edge for strategy d (see Section 3.2). Afterwards, we recalculated the new retracking points for the optimal thresholds (70% for leads and 50% for ice). The results show that if the location of breakpoints is changed by +/− 1 bin, the new retracking point varies by 0.22/0.21 bins (RMSE: 0.27/0.27 bins) for leads and 0.30/0.41 bins (RMSE: 0.41/0.49 bins) for ice floes, respectively, compared with the original retracking point. Hence, the choice of the initial positioning of breakpoints procedure has a non-negligible influence on the tracking point location. However, this can, at least partly, be balanced by the redetermination of the optimal thresholds (which we did not change here).

In Section 4.2, we determined the optimal thresholds for the BCF retracking method by using OIB data from March 2014. To prove the robustness of our selected retracker thresholds (70% for leads and 50% for sea ice), we tested whether the same thresholds are still optimal for the data sets from March 2015 and April 2016.

Table 7. Retrieved sea ice freeboard using BCF with different retracking thresholds relative to the freeboard determined from OIB data acquired in March 2015 and April 2016.

Thresholds (Leads, Ice)	50%, 50%	70%, 50%	70%, 70%	90%, 50%	90%, 70%	90%, 90%	100%, 50%	100%, 70%	100%, 90%
2015
MD 1	14.58	12.75	40.46	14.13	33.46	69.86	19.37	26.82	60.91
RMSE	18.92	16.26	48.77	17.33	43.04	78.27	22.95	37.35	70.69
Correlation	0.3428	0.3486	0.2434	0.3296	0.2328	0.1885	0.2884	0.2152	0.1772
2016
MD 1	19.75	12.85	50.97	14.62	41.63	80.16	16.54	30.06	66.78
RMSE	24.09	17.17	58.34	17.82	50.22	89.75	19.34	39.70	77.82
Correlation	0.3403	0.3968	−0.0109	0.2669	−0.0085	−0.1568	0.3538	−0.0066	−0.1557

¹ MD (Mean freeboard difference): |CS-2−OIB|.

lists the retrieved sea ice freeboard values for different retracking thresholds. As mentioned above, negative freeboards were not filtered when comparing freeboards with different retracking thresholds. In Table 7, the mean difference and RMSE for the threshold combination of (70%, 50%) are smallest also for these two data sets. They are 12.75 cm and 16.26 cm, respectively, relative to the OIB data for March 2015 and 12.85 cm and 17.17 cm for April 2016. This finding indicates that the applied BCF threshold combination is optimal for late winter–early spring measurements (March–April) in the Arctic. Here we can safely assume that the OIB profile data were not affected by moisture in the snow layer during measurements since mean air temperatures of −21.15 °C in April 2016 and −27.15 °C in March 2015 were reported for this area (data from Danish Meteorological Institute). Nevertheless, the snow structure may reveal effects of snow metamorphism processes that have taken place earlier, e.g. in the previous year’s freeze-up period. This, in turn, may shift the main scattering horizon from the snow–ice interface towards the air–snow interface. Another factor that influences the optimal threshold is the roughness of the ice surface within the antenna footprint. In practical measurements, the properties of the snow layer and the ice surface remain unknown (but can to a certain degree be retrieved indirectly by PM approaches, dependent on how realistic the corresponding approach models the snow layer and ice surface roughness). When using fixed thresholds in the retracker, variations of the scattering horizon cannot be accounted for and must be treated as kind of inherent noise or error. Since OIB data are not available for other months, it was not possible to investigate the selection of optimal thresholds for other periods of the year.

6.2. Uncertainty Caused by Snow Loading

For an assessment of the potential impact of snow on the retrieval of ice freeboard on our results we have to consider the differences in the measurements. OIB data are available with spatial resolution of 40 m. The snow depth is determined separately using the snow radar. The data are carefully filtered to take into account different situations (e.g., no snow on leads, moist snow, loss of snow radar signal) and subsequently considered in the calculations of the ice freeboard [29]. In the retrieval of freeboard using the CS-2 data with spatial resolutions of 300 m along-track and 1500 m across-track, the thresholds are determined by the average effect of the snow cover within the footprint, and differences of snow cover are not taken into account (see Section 4.2). Although we have determined the optimal retracking thresholds separately for FYI and MYI, and the same thresholds are obtained, the detailed regional snow cover difference is still not considered. Dependent on snow properties the mean scattering horizon may vary considerably, so that even if the spatial variations of snow depth are known in detail (which they usually are not in the case of satellite data), the position of the scattering horizon may vary between the snow–ice and air–snow interface. Vertical variations of the scattering horizon are larger for thicker snow layers.

Willatt et al. found that the main scattering of Ku-band radar pulses occurs at the air–ice interface when there is no snow cover on the ice [50]. However, when snow covers the ice surface, the main scattering surface is shifted upwards to a horizon between the air–snow and snow–ice interfaces, especially on MYI on which a large snow thickness and complex snow layering structures are typical. Based on the 2008 CryoVEx field experiment, Willatt et al. reported that 80% of Ku band radar returns were closer to the snow–ice interface than to the air–snow interface when the snow was cold and dry [50]. During the CryoVEx 2006 experiment, surface temperatures were close to the melting point and the snow stratigraphy was complex. Therefore, only 25% of Ku-band radar returns were closer to the snow–ice interface. Based on in situ measurements of thermo-physical snow measurements, Nandan et al. theoretically derived the position of the scattering horizon in brine-wetted snow on FYI and found that its shift towards the snow–air interface cannot be neglected [54]. Kurtz et al., for example, reported that Ku-band radar waves did not fully penetrate the snow cover, and the freeboard error associated with the assumption of total penetration was < 0.04 m for FYI, and < 0.08 m for MYI [12]. By comparing year-round measurements of snow depth and ice freeboard data obtained from mass-balance buoys with CS-2 freeboard retrievals, Ricker et al. [53] found that the ice freeboard error due to snow load can be much higher than the errors theoretically estimated by [12].

The mean snow depths obtained from OIB are 15.03 cm and 23.37 cm for FYI and MYI, respectively, in March 2015, and 8.39 cm and 24.87 cm in April 2016. Hence, they are large enough for shifts of the scattering horizon on the order of the RMSE and MD given in Table 6 (2015: FYI 9.03/6.70, MYI 15.80/9.67 2016: FYI 18.59/13.42, MYI 10.50/7.74, all in cm). The correlation coefficients (Table 6) indicate that the spatial variations of the freeboard compare only moderately between the OIB profiles and our data, which can at least be partly attributed to the fact that local snow depth and structure variations cannot be directly incorporated in the freeboard retrieval from satellite data.

Figure 5 shows the Arctic-wide ice freeboard retrieved with BCF, together with the PDF of freeboard differences between BCF and TFMRA and between BCF and L2I, for the data sets from March 2015 and April 2016. Here we used the snow depths from UCL04, since for OIB, we only have single profiles. Snow density are from Warren’s estimation in March and April [48]. From the area around 85° N, 120° W of Arctic Ocean, which is called the “Wingham Box” [56], reliable SAR-mode data cannot be obtained and the SIRAL sensor is mostly switched to the SARIn mode [17]. Hence, all CS-2 data used in this study are located outside of the “Wingham Box”. The major spatial patterns are common to all three retracking methods: highest freeboard values are found north of Greenland (MYI region) and low values in Baffin Bay. At the east coast of Greenland and in the central Arctic the freeboard lies between the values observed in Baffin Bay and north of Greenland. The ice freeboards in the CS-2 L2I products are in general larger than the freeboards by BCF in March 2015 and lower than BCF in April 2016, and BCF ice freeboard is similar to TFMRA in both March 2015 and April 2016 (see Figure 5b,d).

The different results for the sea ice freeboard are due to the use of different retrackers over assumed or modeled fixed snow depths and given mean sea surface heights. Because of varying snow depths, snow structures and ice surface properties, one retracker may be perform better over a region but be worse over another region. For example, Xia and Xie (2017) found that the performance of different retrackers varies for different ice types [55]. Consistent higher or lower ice freeboard for the same ice type may be due to the rather similar CS-2 echo waveforms in the area of interest. The sea ice freeboard we obtained by L2I is the highest in March 2015 (BCF: 0.21 ± 0.10 m, L2I: 0.24 ± 0.07 m, TFMRA: 0.18 ± 0.09 m.), whereas in April 2016, BCF has the highest freeboard (0.19 ± 0.09 m) followed by TFMRA (0.17±0.08 m). Note that Table 5 shows the difference of freeboards along the OIB tracks, whereas here we consider the entire Arctic sea ice cover. The influence of snow on the spatial patterns of freeboard differences between BCF and TFMRA/L2I can only be analyzed in detail when detailed snow property data are available.

6.3. Comparison between BCF and TFMRA

Although BCF and TFMRA revealed rather similar freeboard differences relative to the OIB data (Table 5), the mean difference for BCF is closer to the OIB data than TFMRA (Figure 4). In addition, the MD and RMSE, which are derived from local comparisons of data gridded in cells of 25 km length, indicate a better statistical accuracy for BCF than for TFMRA although the difference is not large (Table 5). This difference shows that with the local increase/decrease of the ice freeboard, BCF freeboard is in better agreement with the OIB data than with the TFMRA data, which in turn seems to indicate that the BCF retracker follows the local elevation changes present in the OIB data more closely than the TFMRA retracker. For correlation, slightly lower values can be found in BCF freeboard than in TFMTA freeboard data. However, due to the small dynamic range of the freeboard values, the correlation coefficient is not a robust measure for the goodness of fit in this case.

In Section 6.1, we found that the applied threshold combination in BCF is optimal both in late winter and early spring measurements in three different years (Table 3). We have used the same thresholds (70%, 50%) also for TFMRA every year since we found that this is the best possible choice (Table 4). From Table 3 and Table 4, we can infer that the difference in the ice freeboard retrieval results between the optimal and the near-optimal thresholds is more sensitive for BCF than for TFMRA.

7. Conclusions

In the study we presented here, a new method is introduced for the retracking correction of the Cryosat-2 radar waveform, which we denote ‘BCF’ (Bézier Curve Fitting). It is based on using composite cubic Bézier curves to match the altimeter echo shape. Included in our investigation is the finding of an optimal scheme for the position of break points between single Bezier curves in the composite. Based on comparison with the OIB sea ice freeboards, snow depth, and thickness quick look data over Arctic sea ice in March 2014, we found that the retracking points should be fixed on positions at which the leading edge reaches 70% of the matching Bézier curve peak for leads and 50% for sea ice, respectively. In addition, we investigated BCF, TFMRA, and L2I products generated from the CS-2 data acquired during two OIB campaign periods in March 2015 and April 2016. The mean absolute differences between retrieved freeboard values and OIB measurement are 9.22 and 7.79 cm for BCF, 10.41 cm and 8.16 cm for TFMRA, and 10.01 cm and 8.42 cm for L2I for profiles that were mainly located in multi-year ice. This suggests that the proposed algorithm has a good potential to improve sea ice freeboard retrieval accuracy. Our method obtained better fitting results than TFMRA and CS-2 L2I methods for the analyzed data set. We plan to study whether this method is compatible with other altimeter (HY-2, Sentinel-3, etc.) records to produce a more complete time series of sea ice thickness records in the future.